Source code for refl1d.dist

Inhomogeneous samples

In the presence of samples with short range order on scale of the coherence
length of the probe in the plane, but long range disorder following some
distribution of parameter values, the reflectivity can be computed from
a weighted incoherent sum of the reflectivities for different values of
the parameter.

DistristributionExperiment allows the model to be computed for a single
varying parameter.  Multi-parameter dispersion models are not available.

import numpy as np
from bumps.parameter import Parameter, to_dict

from .experiment import ExperimentBase

[docs]class Weights(object): """ Parameterized distribution for use in DistributionExperiment. To support non-uniform experiments, we must bin the possible values for the parameter and compute the theory function for one parameter value per bin. The weighted sum of the resulting theory functions is the value that we compare to the data. Performing this analysis requires a cumulative density function which can return the integrated value of the probability density from -inf to x. The total density in each bin is then the difference between the cumulative densities at the edges. If the distribution is wider than the range, then the tails need to be truncated and the bins reweighted to a total density of 1, or the tail density can be added to the first and last bins. Weights of zero are not returned. Note that if the tails are truncated, this may result in no weights being returned. The vector *edges* contains the bin edges for the distribution. The function *cdf* returns the cumulative density function at the edges. The *cdf* function must implement the scipy.stats interface, with function signature f(x, a1, a2, ..., loc=0, scale=1). The list *args* defines the arguments a1, a2, etc. The underlying parameters are available as args[i]. Similarly, *loc* and *scale* define the distribution center and width. Use *truncated=False* if you want the distribution tails to be included in the weights. SciPy distribution D is used by specifying cdf=scipy.stats.D.cdf. Useful distributions include:: norm Gaussian distribution. halfnorm Right half of a gaussian. triang Triangle distribution from loc up to loc+args[0]*scale and down to loc+scale. Use loc=edges[0], scale=edges[-1] and args=[0.5] to define a symmetric triangle in the range of parameter P. uniform Flat from loc to loc+scale. Use loc=edges[0], scale=edges[-1] to define P as uniform over the range. """ def __init__(self, edges=None, cdf=None, args=(), loc=None, scale=None, truncated=True): self.edges = np.asarray(edges) self.cdf = cdf self.truncated = truncated self.loc = Parameter.default(loc) self.scale = Parameter.default(scale) self.args = [Parameter.default(a) for a in args]
[docs] def parameters(self): return {'args': self.args, 'loc': self.loc, 'scale': self.scale}
[docs] def to_dict(self): return to_dict({ 'type': type(self).__name__, 'edges': self.edges.tolist(), 'cdf': self.cdf.__name__, # TODO: can't lookup name 'args': self.args, 'loc': self.loc, 'scale': self.scale, 'truncated': self.truncated, })
def __iter__(self): # Find weights and normalize the sum to 1 centers = (self.edges[:-1]+self.edges[1:])/2 loc = self.loc.value scale = self.scale.value args = [p.value for p in self.args] cumulative_weights = self.cdf(self.edges, *args, loc=loc, scale=scale) if not self.truncated: cumulative_weights[0], cumulative_weights[-1] = 0, 1 relative_weights = cumulative_weights[1:] - cumulative_weights[:-1] total_weight = np.sum(relative_weights) if total_weight == 0: return iter([]) else: weights = relative_weights / total_weight idx = weights > 0 return iter(zip(centers[idx], weights[idx]))
[docs]class DistributionExperiment(ExperimentBase): """ Compute reflectivity from a non-uniform sample. *P* is the target parameter for the model, which takes on the values from *distribution* in the context of the *experiment*. The result is the weighted sum of the theory curves after setting *P.value* to each distribution value. Clearly, *P* should not be a fitted parameter, but the remaining experiment parameters can be fitted, as can the parameters of the distribution. If *coherent* is true, then the reflectivity of the mixture is computed from the coherent sum rather than the incoherent sum. See :class:`Weights` for a description of how to set up the distribution. """ def __init__(self, experiment=None, P=None, distribution=None, coherent=False): self.P = P self.distribution = distribution self.experiment = experiment self.probe = self.experiment.probe self.coherent = coherent self._substrate = self.experiment.sample[0].material self._surface = self.experiment.sample[-1].material self._cache = {} # Cache calculated profiles/reflectivities self._name = None
[docs] def parameters(self): return { 'distribution': self.distribution.parameters(), 'experiment': self.experiment.parameters(), }
[docs] def to_dict(self): return to_dict({ 'type': type(self).__name__, 'P': self.P, # Note: use parameter id to restore 'distribution': self.distribution, 'experiment': self.experiment, # Don't need self.probe since it is the experiment probe. 'coherent': self.coherent, })
[docs] def reflectivity(self, resolution=True, interpolation=0): key = ("reflectivity", resolution, interpolation) if key not in self._cache: calc_R = 0 for x, w in self.distribution: if w > 0: self.P.value = x self.experiment.update() Qx, Rx = self.experiment._reflamp() if self.coherent: calc_R += w*Rx else: calc_R += w*abs(Rx)**2 if self.coherent: calc_R = abs(calc_R)**2 Q, R = self.probe.apply_beam(Qx, calc_R, resolution=resolution, interpolation=interpolation) self._cache[key] = Q, R return self._cache[key]
def _max_P(self): x, w = zip(*self.distribution) idx = np.argmax(w) return x[idx]
[docs] def smooth_profile(self, dz=1): """ Compute a density profile for the material """ key = 'smooth_profile', dz if key not in self._cache: P = self._max_P() if self.P.value != P: self.P.value = P self.experiment.update() self._cache[key] = self.experiment.smooth_profile(dz=dz) return self._cache[key]
[docs] def step_profile(self): """ Compute a scattering length density profile """ key = 'step_profile' if key not in self._cache: P = self._max_P() if self.P.value != P: self.P.value = P self.experiment.update() self._cache[key] = self.experiment.step_profile() return self._cache[key]
[docs] def plot_profile(self, plot_shift=0.): import matplotlib.pyplot as plt from bumps.plotutil import auto_shift trans = auto_shift(plot_shift) z, rho, irho = self.step_profile() plt.plot(z, rho, '-g', z, irho, '-b', transform=trans) z, rho, irho = self.smooth_profile() plt.plot(z, rho, ':g', z, irho, ':b', transform=trans) plt.legend(['rho', 'irho'])
[docs] def plot_weights(self): import matplotlib.pyplot as plt x, w = zip(*self.distribution) plt.stem(x, 100*np.array(w)) plt.title('Weight distribution') plt.xlabel( plt.ylabel('Percentage')