Source code for refnx.reflect.spline

import numpy as np

from scipy.interpolate import PchipInterpolator as Pchip

from refnx.reflect import Structure, Component
from refnx.analysis import Parameter, Parameters, possibly_create_parameter

EPS = np.finfo(float).eps

[docs]class Spline(Component): """ Freeform modelling of the real part of an SLD profile using spline interpolation. Parameters ---------- extent : float or Parameter Total extent of spline region vs : Sequence of float/Parameter the real part of the SLD values of each of the knots. dz : Sequence of float/Parameter the lateral offset between successive knots. name : str Name of component interpolator : scipy.interpolate Univariate Interpolator, optional Which scipy.interpolate Univariate Interpolator to use. zgrad : bool, optional If true then extra control knots are placed outside this spline with the same SLD as the materials on the left and right. With a monotonic interpolator this guarantees that the gradient is zero at either end of the interval. microslab_max_thickness : float Maximum size of the microslabs approximating the spline. Notes ----- This spline component only generates the real part of the SLD (thereby assuming that the imaginary part is negligible). The sequence dz are the lateral offsets of the knots normalised to a unit interval [0, 1]. The reason for using lateral offsets is so that the knots are monotonically increasing in location. When each dz offset is turned into a Parameter it is given bounds in [0, 1]. Thus with an extent of 500, and dz = [0.1, 0.2, 0.2], the knots will be at [0, 50, 150, 250, 500]. Notice that there are two extra knots for the start and end of the interval (disregarding the `zgrad` control knots). If ``np.sum(dz) > 1``, then the knot spacings are normalised to 1. e.g. dz of [0.1, 0.2, 0.9] would result in knots (in the normalised interval) of [0, 0.0833, 0.25, 1, 1]. If `vs` is monotonic then the output spline will be monotonic. If `vs` is not monotonic then there may be regions of the spline larger or smaller than `left` or `right`. The slab representation of this component are approximated using a 'microslab' representation of spline. The max thickness of each microslab is `microslab_max_thickness`. A Spline component should not be used more than once in a given Structure. """ def __init__( self, extent, vs, dz, name="", interpolator=Pchip, zgrad=True, microslab_max_thickness=1, ): super().__init__() = name self.microslab_max_thickness = microslab_max_thickness self.extent = possibly_create_parameter( extent, name="%s - spline extent" % name, units="Å" ) = Parameters(name="dz - spline") for i, z in enumerate(dz): p = possibly_create_parameter(z, name=f"{name} - spline dz[{i}]") p.range(0.0000001, 1) self.vs = Parameters(name="vs - spline") for i, v in enumerate(vs): p = possibly_create_parameter( v, name=f"{name} - spline vs[{i}]", units="10**-6 Å**-2" ) self.vs.append(p) if len(self.vs) != len( raise ValueError("dz and vs must have same number of entries") self.zgrad = zgrad self.interpolator = interpolator self.__cached_interpolator = { "zeds": np.array([]), "vs": np.array([]), "interp": None, "extent": -1, } def __repr__(self): s = ( f"Spline({self.extent!r}, {self.vs!r}, {!r}," f" name={!r}, zgrad={self.zgrad}," f" microslab_max_thickness={self.microslab_max_thickness})" ) return s def _interpolator(self, structure): dz = np.array( zeds = np.cumsum(dz) # if dz's sum to more than 1, then normalise to unit interval. if len(zeds) and zeds[-1] > 1: # there may be no knots zeds /= zeds[-1] zeds = np.clip(zeds, 0, 1) # note - this means you shouldn't use the same Spline more than once in # a Component, because only the first use will be detected. try: loc = structure.index(self) # figure out SLDs for the bracketing Components. # note the use of the modulus operator. This means that if the # Spline is at the end, then the right most Component will be # assumed to be the first Component. This is to aid the use of # Spline in a Stack. left_component = structure[loc - 1] right_component = structure[(loc + 1) % len(structure)] except ValueError: raise ValueError( "Spline didn't appear to be part of a super Structure" ) if isinstance(left_component, Spline) or isinstance( right_component, Spline ): raise ValueError( "Spline must be bracketed by Components that" " aren't Splines." ) vs = np.array(self.vs) left_sld = structure.overall_sld( np.atleast_2d(left_component.slabs(structure)[-1]), structure.solvent, )[..., 1] right_sld = structure.overall_sld( np.atleast_2d(right_component.slabs(structure)[0]), structure.solvent, )[..., 1] if self.zgrad: zeds = np.concatenate([[-1.1, 0 - EPS], zeds, [1 + EPS, 2.1]]) vs = np.concatenate([left_sld, left_sld, vs, right_sld, right_sld]) else: zeds = np.concatenate([[0 - EPS], zeds, [1 + EPS]]) vs = np.concatenate([left_sld, vs, right_sld]) # cache the interpolator cache_zeds = self.__cached_interpolator["zeds"] cache_vs = self.__cached_interpolator["vs"] cache_extent = self.__cached_interpolator["extent"] # you don't need to recreate the interpolator if ( np.array_equal(zeds, cache_zeds) and np.array_equal(vs, cache_vs) and np.equal(self.extent, cache_extent) ): return self.__cached_interpolator["interp"] else: self.__cached_interpolator["zeds"] = zeds self.__cached_interpolator["vs"] = vs self.__cached_interpolator["extent"] = float(self.extent) # TODO make vfp zero for z > self.extent interpolator = self.interpolator(zeds, vs) self.__cached_interpolator["interp"] = interpolator return interpolator
[docs] def __call__(self, z, structure): """ Calculates the spline value at z Parameters ---------- z : float Distance along spline structure: refnx.reflect.Structure Structure hosting this Component Returns ------- sld : float Real part of SLD """ interpolator = self._interpolator(structure) vs = interpolator(z / float(self.extent)) return vs
@property def parameters(self): p = Parameters( p.extend([self.extent,, self.vs]) return p
[docs] def logp(self): return 0
[docs] def slabs(self, structure=None): """ Slab representation of the spline, as an array Parameters ---------- structure : refnx.reflect.Structure The Structure hosting this Component """ if structure is None: raise ValueError("Spline.slabs() requires a valid Structure") num_slabs = np.ceil(float(self.extent) / self.microslab_max_thickness) slab_thick = float(self.extent / num_slabs) slabs = np.zeros((int(num_slabs), 5)) slabs[:, 0] = slab_thick # give last slab a miniscule roughness so it doesn't get contracted slabs[-1:, 3] = 0.5 dist = np.cumsum(slabs[..., 0]) - 0.5 * slab_thick slabs[:, 1] = self(dist, structure) return slabs