"""
Define the resolution functions for the data.
This defines classes for 1D and 2D resolution calculations.
"""
from __future__ import division
import unittest
from scipy.special import erf # type: ignore
from numpy import sqrt, log, log10, exp, pi # type: ignore
import numpy as np # type: ignore
__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D",
"apply_resolution_matrix", "pinhole_resolution", "slit_resolution",
"pinhole_extend_q", "slit_extend_q", "bin_edges",
"interpolate", "linear_extrapolation", "geometric_extrapolation",
]
MINIMUM_RESOLUTION = 1e-8
MINIMUM_ABSOLUTE_Q = 0.02 # relative to the minimum q in the data
# According to (Barker & Pedersen 1995 JAC), 2.5 sigma is a good limit.
# According to simulations with github.com:scattering/sansresolution.git
# it is better to use asymmetric bounds (2.5, 3.0)
PINHOLE_N_SIGMA = (2.5, 3.0)
[docs]class Resolution(object):
"""
Abstract base class defining a 1D resolution function.
*q* is the set of q values at which the data is measured.
*q_calc* is the set of q values at which the theory needs to be evaluated.
This may extend and interpolate the q values.
*apply* is the method to call with I(q_calc) to compute the resolution
smeared theory I(q).
"""
q = None # type: np.ndarray
q_calc = None # type: np.ndarray
[docs] def apply(self, theory):
"""
Smear *theory* by the resolution function, returning *Iq*.
"""
raise NotImplementedError("Subclass does not define the apply function")
[docs]class Perfect1D(Resolution):
"""
Resolution function to use when there is no actual resolution smearing
to be applied. It has the same interface as the other resolution
functions, but returns the identity function.
"""
def __init__(self, q):
self.q_calc = self.q = q
[docs] def apply(self, theory):
return theory
[docs]class Pinhole1D(Resolution):
r"""
Pinhole aperture with q-dependent gaussian resolution.
*q* points at which the data is measured.
*q_width* gaussian 1-sigma resolution at each data point.
*q_calc* is the list of points to calculate, or None if this should
be estimated from the *q* and *q_width*.
*nsigma* is the width of the resolution function. Should be 2.5.
See :func:`pinhole_resolution` for details.
"""
def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA):
#*min_step* is the minimum point spacing to use when computing the
#underlying model. It should be on the order of
#$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes
#are computed with sufficient density to avoid aliasing effects.
# Protect against calls with q_width=0. The extend_q function will
# not extend the q if q_width is 0, but q_width must be non-zero when
# constructing the weight matrix to avoid division by zero errors.
# In practice this should never be needed, since resolution should
# default to Perfect1D if the pinhole geometry is not defined.
self.q, self.q_width = q, q_width
self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
if q_calc is None else np.sort(q_calc))
# Protect against models which are not defined for very low q. Limit
# the smallest q value evaluated (in absolute) to 0.02*min
cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
# Build weight matrix from calculated q values
self.weight_matrix = pinhole_resolution(
self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION),
nsigma=nsigma)
# Force positive q, even for events measured on the opposite side of
# the beam stop.
self.q_calc = abs(self.q_calc)
[docs] def apply(self, theory):
return apply_resolution_matrix(self.weight_matrix, theory)
[docs]class Slit1D(Resolution):
"""
Slit aperture with resolution function.
*q* points at which the data is measured.
*qx_width* slit width in qx
*qy_width* slit height in qy
*q_calc* is the list of points to calculate, or None if this should
be estimated from the *q* and *q_width*.
The *weight_matrix* is computed by :func:`slit_resolution`
"""
def __init__(self, q, qx_width, qy_width=0., q_calc=None):
# Remember what width/dqy was used even though we won't need them
# after the weight matrix is constructed
self.qx_width, self.qy_width = qx_width, qy_width
# Allow independent resolution on each point even though it is not
# needed in practice.
if np.isscalar(qx_width):
qx_width = np.ones(len(q))*qx_width
else:
qx_width = np.asarray(qx_width)
if np.isscalar(qy_width):
qy_width = np.ones(len(q))*qy_width
else:
qy_width = np.asarray(qy_width)
self.q = q.flatten()
self.q_calc = slit_extend_q(q, qx_width, qy_width) \
if q_calc is None else np.sort(q_calc)
# Protect against models which are not defined for very low q. Limit
# the smallest q value evaluated (in absolute) to 0.02*min
cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
# Build weight matrix from calculated q values
self.weight_matrix = \
slit_resolution(self.q_calc, self.q, qx_width, qy_width)
self.q_calc = abs(self.q_calc)
[docs] def apply(self, theory):
return apply_resolution_matrix(self.weight_matrix, theory)
[docs]def apply_resolution_matrix(weight_matrix, theory):
"""
Apply the resolution weight matrix to the computed theory function.
"""
#print("apply shapes", theory.shape, weight_matrix.shape)
Iq = np.dot(theory[None, :], weight_matrix)
#print("result shape",Iq.shape)
return Iq.flatten()
[docs]def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA):
r"""
Compute the convolution matrix *W* for pinhole resolution 1-D data.
Each row *W[i]* determines the normalized weight that the corresponding
points *q_calc* contribute to the resolution smeared point *q[i]*. Given
*W*, the resolution smearing can be computed using *dot(W,q)*.
Note that resolution is limited to $\pm 2.5 \sigma$.[1] The true resolution
function is a broadened triangle, and does not extend over the entire
range $(-\infty, +\infty)$. It is important to impose this limitation
since some models fall so steeply that the weighted value in gaussian
tails would otherwise dominate the integral.
*q_calc* must be increasing. *q_width* must be greater than zero.
[1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects
in Radially Symmetric Small-Angle Neutron Scattering by Numerical and
Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14.
https://doi.org/10.1107/S0021889894010095.
"""
# The current algorithm is a midpoint rectangle rule. In the test case,
# neither trapezoid nor Simpson's rule improved the accuracy.
edges = bin_edges(q_calc)
#edges[edges < 0.0] = 0.0 # clip edges below zero
cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
weights = cdf[1:] - cdf[:-1]
# Limit q range to (-2.5,+3) sigma
try:
nsigma_low, nsigma_high = nsigma
except TypeError:
nsigma_low = nsigma_high = nsigma
qhigh = q + nsigma_high*q_width
qlow = q - nsigma_low*q_width # linear limits
##qlow = q*q/qhigh # log limits
weights[q_calc[:, None] < qlow[None, :]] = 0.
weights[q_calc[:, None] > qhigh[None, :]] = 0.
weights /= np.sum(weights, axis=0)[None, :]
return weights
[docs]def slit_resolution(q_calc, q, width, height, n_height=30):
r"""
Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
$q_\perp$ = *width* and $q_\parallel$ = *height*. *n_height* is
is the number of steps to use in the integration over $q_\parallel$
when both $q_\perp$ and $q_\parallel$ are non-zero.
Each $q$ can have an independent width and height value even though
current instruments use the same slit setting for all measured points.
If slit height is large relative to width, use:
.. math::
I_s(q_i) = \frac{1}{\Delta q_\perp}
\int_0^{\Delta q_\perp}
I\left(\sqrt{q_i^2 + q_\perp^2}\right) \,dq_\perp
If slit width is large relative to height, use:
.. math::
I_s(q_i) = \frac{1}{2 \Delta q_\parallel}
\int_{-\Delta q_\parallel}^{\Delta q_\parallel}
I\left(|q_i + q_\parallel|\right) \,dq_\parallel
For a mixture of slit width and height use:
.. math::
I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp}
\int_{-\Delta q_\parallel}^{\Delta q_\parallel}
\int_0^{\Delta q_\perp}
I\left(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}\right)
\,dq_\perp dq_\parallel
**Definition**
We are using the mid-point integration rule to assign weights to each
element of a weight matrix $W$ so that
.. math::
I_s(q) = W\,I(q_\text{calc})
If *q_calc* is at the mid-point, we can infer the bin edges from the
pairwise averages of *q_calc*, adding the missing edges before
*q_calc[0]* and after *q_calc[-1]*.
For $q_\parallel = 0$, the smeared value can be computed numerically
using the $u$ substitution
.. math::
u_j = \sqrt{q_j^2 - q^2}
This gives
.. math::
I_s(q) \approx \sum_j I(u_j) \Delta u_j
where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the
difference between consecutive edges which have been first converted
to $u$. Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds
to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so
.. math::
W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j
= \frac{1}{\Delta q_\perp} \left(
\sqrt{q_{j+1}^2 - q_i^2} - \sqrt{q_j^2 - q_i^2} \right)
\ \text{if}\ q_j \in \left[q_i, \sqrt{q_i^2 + q_\perp^2}\right]
where $I_s(q_i)$ is the theory function being computed and $q_j$ are the
mid-points between the calculated values in *q_calc*. We tweak the
edges of the initial and final intervals so that they lie on integration
limits.
(To be precise, the transformed midpoint $u(q_j)$ is not necessarily the
midpoint of the edges $u((q_{j-1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$,
but it is at least in the interval, so the approximation is going to be
a little better than the left or right Riemann sum, and should be
good enough for our purposes.)
For $q_\perp = 0$, the $u$ substitution is simpler:
.. math::
u_j = \left|q_j - q\right|
so
.. math::
W_{ij} = \frac{1}{2 \Delta q_\parallel} \Delta u_j
= \frac{1}{2 \Delta q_\parallel} (q_{j+1} - q_j)
\ \text{if}\ q_j \in
\left[q-\Delta q_\parallel, q+\Delta q_\parallel\right]
However, we need to support cases were $u_j < 0$, which means using
$2 (q_{j+1} - q_j)$ when $q_j \in \left[0, q_\parallel-q_i\right]$.
This is not an issue for $q_i > q_\parallel$.
For both $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional
integration with
.. math::
u_{jk} = \sqrt{q_j^2 - (q + (k\Delta q_\parallel/L))^2}
\ \text{for}\ k = -L \ldots L
for $L$ = *n_height*. This gives
.. math::
W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel}
\sum_{k=-L}^L \Delta u_{jk}
\left(\frac{\Delta q_\parallel}{2 L + 1}\right)
"""
#np.set_printoptions(precision=6, linewidth=10000)
# The current algorithm is a midpoint rectangle rule.
q_edges = bin_edges(q_calc) # Note: requires q > 0
#q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
weights = np.zeros((len(q), len(q_calc)), 'd')
#print(q_calc)
for i, (qi, w, h) in enumerate(zip(q, width, height)):
if w == 0. and h == 0.:
# Perfect resolution, so return the theory value directly.
# Note: assumes that q is a subset of q_calc. If qi need not be
# in q_calc, then we can do a weighted interpolation by looking
# up qi in q_calc, then weighting the result by the relative
# distance to the neighbouring points.
weights[i, :] = (q_calc == qi)
elif h == 0:
weights[i, :] = _q_perp_weights(q_edges, qi, w)
elif w == 0:
in_x = 1.0 * ((q_calc >= qi-h) & (q_calc <= qi+h))
abs_x = 1.0*(q_calc < abs(qi - h)) if qi < h else 0.
#print(qi - h, qi + h)
#print(in_x + abs_x)
weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
else:
for k in range(-n_height, n_height+1):
weights[i, :] += _q_perp_weights(q_edges, qi+k*h/n_height, w)
weights[i, :] /= 2*n_height + 1
return weights.T
def _q_perp_weights(q_edges, qi, w):
# Convert bin edges from q to u
u_limit = np.sqrt(qi**2 + w**2)
u_edges = q_edges**2 - qi**2
u_edges[q_edges < abs(qi)] = 0.
u_edges[q_edges > u_limit] = u_limit**2 - qi**2
weights = np.diff(np.sqrt(u_edges))/w
#print("i, qi",i,qi,qi+width)
#print(q_calc)
#print(weights)
return weights
[docs]def pinhole_extend_q(q, q_width, nsigma=PINHOLE_N_SIGMA):
"""
Given *q* and *q_width*, find a set of sampling points *q_calc* so
that each point $I(q)$ has sufficient support from the underlying
function.
"""
try:
nsigma_low, nsigma_high = nsigma
except TypeError:
nsigma_low = nsigma_high = nsigma
q_min, q_max = np.min(q - nsigma_low*q_width), np.max(q + nsigma_high*q_width)
return linear_extrapolation(q, q_min, q_max)
[docs]def slit_extend_q(q, width, height):
"""
Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
that each point I(q) has sufficient support from the underlying
function.
"""
q_min, q_max = np.min(q-height), np.max(np.sqrt((q+height)**2 + width**2))
return geometric_extrapolation(q, q_min, q_max)
[docs]def bin_edges(x):
"""
Determine bin edges from bin centers, assuming that edges are centered
between the bins.
Note: this uses the arithmetic mean, which may not be appropriate for
log-scaled data.
"""
if len(x) < 2 or (np.diff(x) < 0).any():
raise ValueError("Expected bins to be an increasing set")
edges = np.hstack([
x[0] - 0.5*(x[1] - x[0]), # first point minus half first interval
0.5*(x[1:] + x[:-1]), # mid points of all central intervals
x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
])
return edges
[docs]def interpolate(q, max_step):
"""
Returns *q_calc* with points spaced at most max_step apart.
"""
step = np.diff(q)
index = step > max_step
if np.any(index):
inserts = []
for q_i, step_i in zip(q[:-1][index], step[index]):
n = np.ceil(step_i/max_step)
inserts.extend(q_i + np.arange(1, n)*(step_i/n))
# Extend a couple of fringes beyond the end of the data
inserts.extend(q[-1] + np.arange(1, 8)*max_step)
q_calc = np.sort(np.hstack((q, inserts)))
else:
q_calc = q
return q_calc
############################################################################
# unit tests
############################################################################
def eval_form(q, form, pars):
"""
Return the SAS model evaluated at *q*.
*form* is the SAS model returned from :fun:`core.load_model`.
*pars* are the parameter values to use when evaluating.
"""
from sasmodels import direct_model
kernel = form.make_kernel([q])
theory = direct_model.call_kernel(kernel, pars)
kernel.release()
return theory
def gaussian(q, q0, dq, nsigma=2.5):
"""
Return the truncated Gaussian resolution function.
*q0* is the center, *dq* is the width and *q* are the points to evaluate.
"""
# Calculate the density of the tails; the resulting gaussian needs to be
# scaled by this amount in ordere to integrate to 1.0
two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2
return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density)
def romberg_slit_1d(q, width, height, form, pars):
"""
Romberg integration for slit resolution.
This is an adaptive integration technique. It is called with settings
that make it slow to evaluate but give it good accuracy.
"""
from scipy.integrate import romberg # type: ignore
par_set = {p.name for p in form.info.parameters.call_parameters}
if any(k not in par_set for k in pars.keys()):
extra = set(pars.keys()) - par_set
raise ValueError("bad parameters: [%s] not in [%s]"
% (", ".join(sorted(extra)),
", ".join(sorted(pars.keys()))))
if np.isscalar(width):
width = [width]*len(q)
if np.isscalar(height):
height = [height]*len(q)
_int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
_int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
# If both width and height are defined, then it is too slow to use dblquad.
# Instead use trapz on a fixed grid, interpolated into the I(Q) for
# the extended Q range.
#_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height))
Iq = eval_form(q_calc, form, pars)
result = np.empty(len(q))
for i, (qi, w, h) in enumerate(zip(q, width, height)):
if h == 0.:
total = romberg(_int_w, 0, w, args=(qi,),
divmax=100, vec_func=True, tol=0, rtol=1e-8)
result[i] = total/w
elif w == 0.:
total = romberg(_int_h, -h, h, args=(qi,),
divmax=100, vec_func=True, tol=0, rtol=1e-8)
result[i] = total/(2*h)
else:
w_grid = np.linspace(0, w, 21)[None, :]
h_grid = np.linspace(-h, h, 23)[:, None]
u_sub = sqrt((qi+h_grid)**2 + w_grid**2)
f_at_u = np.interp(u_sub, q_calc, Iq)
#print(np.trapz(Iu, w_grid, axis=1))
total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), h_grid[:, 0])
result[i] = total / (2*h*w)
# from scipy.integrate import dblquad
# r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
# args=(qi,))
# result[i] = r/(w*2*h)
# r should be [float, ...], but it is [array([float]), array([float]),...]
return result
def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5):
"""
Romberg integration for pinhole resolution.
This is an adaptive integration technique. It is called with settings
that make it slow to evaluate but give it good accuracy.
"""
from scipy.integrate import romberg # type: ignore
par_set = {p.name for p in form.info.parameters.call_parameters}
if any(k not in par_set for k in pars.keys()):
extra = set(pars.keys()) - par_set
raise ValueError("bad parameters: [%s] not in [%s]"
% (", ".join(sorted(extra)),
", ".join(sorted(pars.keys()))))
func = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
total = [romberg(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
args=(qi, dqi), divmax=100, vec_func=True,
tol=0, rtol=1e-8)
for qi, dqi in zip(q, q_width)]
return np.asarray(total).flatten()
class ResolutionTest(unittest.TestCase):
"""
Test the resolution calculations.
"""
def setUp(self):
self.x = 0.001*np.arange(1, 11)
self.y = self.Iq(self.x)
def Iq(self, q):
"Linear function for resolution unit test"
return 12.0 - 1000.0*q
def test_perfect(self):
"""
Perfect resolution and no smearing.
"""
resolution = Perfect1D(self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
np.testing.assert_equal(output, self.y)
def test_slit_zero(self):
"""
Slit smearing with perfect resolution.
"""
resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
np.testing.assert_equal(output, self.y)
@unittest.skip("not yet supported")
def test_slit_high(self):
"""
Slit smearing with height 0.005
"""
resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
answer = [
9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
]
np.testing.assert_allclose(output, answer, atol=1e-4)
@unittest.skip("not yet supported")
def test_slit_both_high(self):
"""
Slit smearing with width < 100*height.
"""
q = np.logspace(-4, -1, 10)
resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf)
theory = 1000*self.Iq(resolution.q_calc**4)
output = resolution.apply(theory)
answer = [
8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
]
np.testing.assert_allclose(output, answer, atol=1e-4)
@unittest.skip("not yet supported")
def test_slit_wide(self):
"""
Slit smearing with width 0.0002
"""
resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
answer = [
11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
]
np.testing.assert_allclose(output, answer, atol=1e-4)
@unittest.skip("not yet supported")
def test_slit_both_wide(self):
"""
Slit smearing with width > 100*height.
"""
resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001,
q_calc=self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
answer = [
11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
]
np.testing.assert_allclose(output, answer, atol=1e-4)
def test_pinhole_zero(self):
"""
Pinhole smearing with perfect resolution
"""
resolution = Pinhole1D(self.x, 0.0*self.x)
theory = self.Iq(resolution.q_calc)
output = resolution.apply(theory)
np.testing.assert_equal(output, self.y)
# TODO: turn pinhole/slit demos into tests
@unittest.skip("suppress comparison with old version; pinhole calc changed")
def test_pinhole(self):
"""
Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001]
"""
resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x),
q_calc=self.x)
theory = 12.0-1000.0*resolution.q_calc
output = resolution.apply(theory)
# Note: answer came from output of previous run. Non-integer
# values at ends come from the fact that q_calc does not
# extend beyond q, and so the weights don't balance.
answer = [
10.47037734, 9.86925860,
9., 8., 7., 6., 5., 4.,
3.13074140, 2.52962266,
]
np.testing.assert_allclose(output, answer, atol=1e-8)
class IgorComparisonTest(unittest.TestCase):
"""
Test resolution calculations against those returned by Igor.
"""
def setUp(self):
self.pars = TEST_PARS_PINHOLE_SPHERE
from sasmodels import core
self.model = core.load_model("sphere", dtype='double')
def _eval_sphere(self, pars, resolution):
from sasmodels import direct_model
kernel = self.model.make_kernel([resolution.q_calc])
theory = direct_model.call_kernel(kernel, pars)
result = resolution.apply(theory)
kernel.release()
return result
def _compare(self, q, output, answer, tolerance):
#err = (output - answer)/answer
#idx = abs(err) >= tolerance
#problem = zip(q[idx], output[idx], answer[idx], err[idx])
#print("\n".join(str(v) for v in problem))
np.testing.assert_allclose(output, answer, rtol=tolerance)
def test_perfect(self):
"""
Compare sphere model with NIST Igor SANS, no resolution smearing.
"""
pars = TEST_PARS_SLIT_SPHERE
data_string = TEST_DATA_SLIT_SPHERE
data = np.loadtxt(data_string.split('\n')).T
q, _, answer, _ = data
resolution = Perfect1D(q)
output = self._eval_sphere(pars, resolution)
self._compare(q, output, answer, 1e-6)
@unittest.skip("suppress comparison with old version; pinhole calc changed")
def test_pinhole(self):
"""
Compare pinhole resolution smearing with NIST Igor SANS
"""
pars = TEST_PARS_PINHOLE_SPHERE
data_string = TEST_DATA_PINHOLE_SPHERE
data = np.loadtxt(data_string.split('\n')).T
q, q_width, answer = data
resolution = Pinhole1D(q, q_width)
output = self._eval_sphere(pars, resolution)
# TODO: relative error should be lower
self._compare(q, output, answer, 3e-4)
@unittest.skip("suppress comparison with old version; pinhole calc changed")
def test_pinhole_romberg(self):
"""
Compare pinhole resolution smearing with romberg integration result.
"""
pars = TEST_PARS_PINHOLE_SPHERE
data_string = TEST_DATA_PINHOLE_SPHERE
pars['radius'] *= 5
data = np.loadtxt(data_string.split('\n')).T
q, q_width, answer = data
answer = romberg_pinhole_1d(q, q_width, self.model, pars)
## Getting 0.1% requires 5 sigma and 200 points per fringe
#q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
# 2*np.pi/pars['radius']/200)
#tol = 0.001
## The default 2.5 sigma and no extra points gets 1%
q_calc = None # type: np.ndarray
tol = 0.01
resolution = Pinhole1D(q, q_width, q_calc=q_calc)
output = self._eval_sphere(pars, resolution)
if 0: # debug plot
import matplotlib.pyplot as plt # type: ignore
resolution = Perfect1D(q)
source = self._eval_sphere(pars, resolution)
plt.loglog(q, source, '.')
plt.loglog(q, answer, '-', hold=True)
plt.loglog(q, output, '-', hold=True)
plt.show()
self._compare(q, output, answer, tol)
def test_slit(self):
"""
Compare slit resolution smearing with NIST Igor SANS
"""
pars = TEST_PARS_SLIT_SPHERE
data_string = TEST_DATA_SLIT_SPHERE
data = np.loadtxt(data_string.split('\n')).T
q, delta_qv, _, answer = data
resolution = Slit1D(q, qx_width=delta_qv, qy_width=0)
output = self._eval_sphere(pars, resolution)
# TODO: eliminate Igor test since it is too inaccurate to be useful.
# This means we can eliminate the test data as well, and instead
# use a generated q vector.
self._compare(q, output, answer, 0.5)
def test_slit_romberg(self):
"""
Compare slit resolution smearing with romberg integration result.
"""
pars = TEST_PARS_SLIT_SPHERE
data_string = TEST_DATA_SLIT_SPHERE
data = np.loadtxt(data_string.split('\n')).T
q, delta_qv, _, answer = data
answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars)
q_calc = slit_extend_q(interpolate(q, 2*np.pi/pars['radius']/20),
delta_qv[0], 0.)
resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc)
output = self._eval_sphere(pars, resolution)
# TODO: relative error should be lower
self._compare(q, output, answer, 0.025)
def test_ellipsoid(self):
"""
Compare romberg integration for ellipsoid model.
"""
from .core import load_model
pars = {
'scale':0.05,
'radius_polar':500, 'radius_equatorial':15000,
'sld':6, 'sld_solvent': 1,
}
form = load_model('ellipsoid', dtype='double')
q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
width, height = 0.117, 0.
resolution = Slit1D(q, qx_width=width, qy_width=height)
answer = romberg_slit_1d(q, width, height, form, pars)
output = resolution.apply(eval_form(resolution.q_calc, form, pars))
# TODO: 10% is too much error; use better algorithm
#print(np.max(abs(answer-output)/answer))
self._compare(q, output, answer, 0.1)
#TODO: can sas q spacing be too sparse for the resolution calculation?
@unittest.skip("suppress sparse data test; not supported by current code")
def test_pinhole_sparse(self):
"""
Compare pinhole resolution smearing with NIST Igor SANS on sparse data
"""
pars = TEST_PARS_PINHOLE_SPHERE
data_string = TEST_DATA_PINHOLE_SPHERE
data = np.loadtxt(data_string.split('\n')).T
q, q_width, answer = data[:, ::20] # Take every nth point
resolution = Pinhole1D(q, q_width)
output = self._eval_sphere(pars, resolution)
self._compare(q, output, answer, 1e-6)
# pinhole sphere parameters
TEST_PARS_PINHOLE_SPHERE = {
'scale': 1.0, 'background': 0.01,
'radius': 60.0, 'sld': 1, 'sld_solvent': 6.3,
}
# Q, dQ, I(Q) calculated by NIST Igor SANS package
TEST_DATA_PINHOLE_SPHERE = """\
0.001278 0.0002847 2538.41176383
0.001562 0.0002905 2536.91820405
0.001846 0.0002956 2535.13182479
0.002130 0.0003017 2533.06217813
0.002414 0.0003087 2530.70378586
0.002698 0.0003165 2528.05024192
0.002982 0.0003249 2525.10408349
0.003266 0.0003340 2521.86667499
0.003550 0.0003437 2518.33907750
0.003834 0.0003539 2514.52246995
0.004118 0.0003646 2510.41798319
0.004402 0.0003757 2506.02690988
0.004686 0.0003872 2501.35067884
0.004970 0.0003990 2496.38678318
0.005253 0.0004112 2491.16237596
0.005537 0.0004237 2485.63911673
0.005821 0.0004365 2479.83657083
0.006105 0.0004495 2473.75676948
0.006389 0.0004628 2467.40145990
0.006673 0.0004762 2460.77293372
0.006957 0.0004899 2453.86724627
0.007241 0.0005037 2446.69623838
0.007525 0.0005177 2439.25775219
0.007809 0.0005318 2431.55421398
0.008093 0.0005461 2423.58785521
0.008377 0.0005605 2415.36158137
0.008661 0.0005750 2406.87009473
0.008945 0.0005896 2398.12841186
0.009229 0.0006044 2389.13360806
0.009513 0.0006192 2379.88958042
0.009797 0.0006341 2370.39776774
0.010080 0.0006491 2360.69528793
0.010360 0.0006641 2350.85169027
0.010650 0.0006793 2340.42023633
0.010930 0.0006945 2330.11206013
0.011220 0.0007097 2319.20109972
0.011500 0.0007251 2308.43503981
0.011780 0.0007404 2297.44820179
0.012070 0.0007558 2285.83853677
0.012350 0.0007713 2274.41290746
0.012640 0.0007868 2262.36219581
0.012920 0.0008024 2250.51169731
0.013200 0.0008180 2238.45596231
0.013490 0.0008336 2225.76495666
0.013770 0.0008493 2213.29618391
0.014060 0.0008650 2200.19110751
0.014340 0.0008807 2187.34050325
0.014620 0.0008965 2174.30529864
0.014910 0.0009123 2160.61632548
0.015190 0.0009281 2147.21038112
0.015470 0.0009440 2133.62023580
0.015760 0.0009598 2119.37907426
0.016040 0.0009757 2105.45234903
0.016330 0.0009916 2090.86319102
0.016610 0.0010080 2076.60576032
0.016890 0.0010240 2062.19214565
0.017180 0.0010390 2047.10550219
0.017460 0.0010550 2032.38715621
0.017740 0.0010710 2017.52560123
0.018030 0.0010880 2001.99124318
0.018310 0.0011040 1986.84662060
0.018600 0.0011200 1971.03389745
0.018880 0.0011360 1955.61395119
0.019160 0.0011520 1940.08291563
0.019450 0.0011680 1923.87672225
0.019730 0.0011840 1908.10656374
0.020020 0.0012000 1891.66297192
0.020300 0.0012160 1875.66789021
0.020580 0.0012320 1859.56357196
0.020870 0.0012490 1842.79468290
0.021150 0.0012650 1826.50064489
0.021430 0.0012810 1810.11533702
0.021720 0.0012970 1793.06840882
0.022000 0.0013130 1776.51153580
0.022280 0.0013290 1759.87201249
0.022570 0.0013460 1742.57354412
0.022850 0.0013620 1725.79397319
0.023140 0.0013780 1708.35831550
0.023420 0.0013940 1691.45256069
0.023700 0.0014110 1674.48561783
0.023990 0.0014270 1656.86525366
0.024270 0.0014430 1639.79847285
0.024550 0.0014590 1622.68887088
0.024840 0.0014760 1604.96421100
0.025120 0.0014920 1587.85768129
0.025410 0.0015080 1569.99297335
0.025690 0.0015240 1552.84580279
0.025970 0.0015410 1535.54074115
0.026260 0.0015570 1517.75249337
0.026540 0.0015730 1500.40115023
0.026820 0.0015900 1483.03632237
0.027110 0.0016060 1465.05942429
0.027390 0.0016220 1447.67682181
0.027670 0.0016390 1430.46495191
0.027960 0.0016550 1412.49232282
0.028240 0.0016710 1395.13182318
0.028520 0.0016880 1377.93439837
0.028810 0.0017040 1359.99528971
0.029090 0.0017200 1342.67274512
0.029370 0.0017370 1325.55375609
"""
# Slit sphere parameters
TEST_PARS_SLIT_SPHERE = {
'scale': 0.01, 'background': 0.01,
'radius': 60000, 'sld': 1, 'sld_solvent': 4,
}
# Q dQ I(Q) I_smeared(Q)
TEST_DATA_SLIT_SPHERE = """\
2.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
2.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
2.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
3.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
3.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
3.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
4.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
5.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
5.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
6.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
6.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
7.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
7.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
8.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
8.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
9.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
1.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
1.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
1.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
1.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
1.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
1.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
1.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
1.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
1.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
1.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
2.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
2.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
2.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
2.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
2.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
2.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
2.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
2.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
3.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
3.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
3.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
3.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
4.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
4.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
4.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
4.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
5.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
5.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
6.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
6.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
7.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
7.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
8.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
8.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
9.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
9.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
1.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
1.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
1.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
1.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
1.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
1.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
1.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
1.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
1.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
1.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
2.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
2.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
2.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
2.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
2.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
2.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
"""
def main():
"""
Run tests given is sys.argv.
Returns 0 if success or 1 if any tests fail.
"""
import sys
import xmlrunner # type: ignore
suite = unittest.TestSuite()
suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__]))
runner = xmlrunner.XMLTestRunner(output='logs')
result = runner.run(suite)
return 1 if result.failures or result.errors else 0
############################################################################
# usage demo
############################################################################
def _eval_demo_1d(resolution, title):
import sys
from sasmodels import core
from sasmodels import direct_model
name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder'
if name == 'cylinder':
pars = {'length':210, 'radius':500, 'background': 0}
elif name == 'teubner_strey':
pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643}
elif name in ('sphere', 'spherepy'):
pars = TEST_PARS_SLIT_SPHERE
elif name == 'ellipsoid':
pars = {
'scale':0.05, 'background': 0,
'r_polar':500, 'r_equatorial':15000,
'sld':6, 'sld_solvent': 1,
}
else:
pars = {}
model_info = core.load_model_info(name)
model = core.build_model(model_info)
kernel = model.make_kernel([resolution.q_calc])
theory = direct_model.call_kernel(kernel, pars)
Iq = resolution.apply(theory)
if isinstance(resolution, Slit1D):
width, height = resolution.qx_width, resolution.qy_width
Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
else:
dq = resolution.q_width
Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars)
import matplotlib.pyplot as plt # type: ignore
plt.loglog(resolution.q_calc, theory, label='unsmeared')
plt.loglog(resolution.q, Iq, label='smeared', hold=True)
plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True)
plt.legend()
plt.title(title)
plt.xlabel("Q (1/Ang)")
plt.ylabel("I(Q) (1/cm)")
def demo_pinhole_1d():
"""
Show example of pinhole smearing.
"""
q = np.logspace(-4, np.log10(0.2), 400)
q_width = 0.1*q
resolution = Pinhole1D(q, q_width)
_eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
def demo_slit_1d():
"""
Show example of slit smearing.
"""
q = np.logspace(-4, np.log10(0.2), 100)
w = h = 0.
#w = 0.000000277790
w = 0.0277790
#h = 0.00277790
#h = 0.0277790
resolution = Slit1D(q, w, h)
_eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
def demo():
"""
Run the resolution demos.
"""
import matplotlib.pyplot as plt # type: ignore
plt.subplot(121)
demo_pinhole_1d()
#plt.yscale('linear')
plt.subplot(122)
demo_slit_1d()
#plt.yscale('linear')
plt.show()
if __name__ == "__main__":
#demo()
main()