# source:sasmodels/sasmodels/resolution.py@0b9c6df

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0b9c6df was 0b9c6df, checked in by Paul Kienzle <pkienzle@…>, 2 years ago

suppress pinhole test for now

• Property mode set to 100644
File size: 41.6 KB
Line
1"""
2Define the resolution functions for the data.
3
4This defines classes for 1D and 2D resolution calculations.
5"""
6from __future__ import division
7
8import unittest
9
10from scipy.special import erf  # type: ignore
11from numpy import sqrt, log, log10, exp, pi  # type: ignore
12import numpy as np  # type: ignore
13
14__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D",
15           "apply_resolution_matrix", "pinhole_resolution", "slit_resolution",
16           "pinhole_extend_q", "slit_extend_q", "bin_edges",
17           "interpolate", "linear_extrapolation", "geometric_extrapolation",
18          ]
19
20MINIMUM_RESOLUTION = 1e-8
21MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data
22PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC
23
24class Resolution(object):
25    """
26    Abstract base class defining a 1D resolution function.
27
28    *q* is the set of q values at which the data is measured.
29
30    *q_calc* is the set of q values at which the theory needs to be evaluated.
31    This may extend and interpolate the q values.
32
33    *apply* is the method to call with I(q_calc) to compute the resolution
34    smeared theory I(q).
35    """
36    q = None  # type: np.ndarray
37    q_calc = None  # type: np.ndarray
38    def apply(self, theory):
39        """
40        Smear *theory* by the resolution function, returning *Iq*.
41        """
42        raise NotImplementedError("Subclass does not define the apply function")
43
44
45class Perfect1D(Resolution):
46    """
47    Resolution function to use when there is no actual resolution smearing
48    to be applied.  It has the same interface as the other resolution
49    functions, but returns the identity function.
50    """
51    def __init__(self, q):
52        self.q_calc = self.q = q
53
54    def apply(self, theory):
55        return theory
56
57
58class Pinhole1D(Resolution):
59    r"""
60    Pinhole aperture with q-dependent gaussian resolution.
61
62    *q* points at which the data is measured.
63
64    *q_width* gaussian 1-sigma resolution at each data point.
65
66    *q_calc* is the list of points to calculate, or None if this should
67    be estimated from the *q* and *q_width*.
68
69    *nsigma* is the width of the resolution function.  Should be 2.5.
70    See :func:pinhole_resolution for details.
71    """
72    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA):
73        #*min_step* is the minimum point spacing to use when computing the
74        #underlying model.  It should be on the order of
75        #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes
76        #are computed with sufficient density to avoid aliasing effects.
77
78        # Protect against calls with q_width=0.  The extend_q function will
79        # not extend the q if q_width is 0, but q_width must be non-zero when
80        # constructing the weight matrix to avoid division by zero errors.
81        # In practice this should never be needed, since resolution should
82        # default to Perfect1D if the pinhole geometry is not defined.
83        self.q, self.q_width = q, q_width
84        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
85                       if q_calc is None else np.sort(q_calc))
86
87        # Protect against models which are not defined for very low q.  Limit
88        # the smallest q value evaluated to 0.02*min.  Note that negative q
89        # values are trimmed even for broad resolution.  Although not possible
90        # from the geometry, they may appear since we are using a truncated
91        # gaussian to represent resolution rather than a skew distribution.
92        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
93        self.q_calc = self.q_calc[self.q_calc >= cutoff]
94
95        # Build weight matrix from calculated q values
96        self.weight_matrix = pinhole_resolution(
97            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION),
98            nsigma=nsigma)
99
100    def apply(self, theory):
101        return apply_resolution_matrix(self.weight_matrix, theory)
102
103
104class Slit1D(Resolution):
105    """
106    Slit aperture with resolution function.
107
108    *q* points at which the data is measured.
109
110    *qx_width* slit width in qx
111
112    *qy_width* slit height in qy
113
114    *q_calc* is the list of points to calculate, or None if this should
115    be estimated from the *q* and *q_width*.
116
117    The *weight_matrix* is computed by :func:slit1d_resolution
118    """
119    def __init__(self, q, qx_width, qy_width=0., q_calc=None):
120        # Remember what width/dqy was used even though we won't need them
121        # after the weight matrix is constructed
122        self.qx_width, self.qy_width = qx_width, qy_width
123
124        # Allow independent resolution on each point even though it is not
125        # needed in practice.
126        if np.isscalar(qx_width):
127            qx_width = np.ones(len(q))*qx_width
128        else:
129            qx_width = np.asarray(qx_width)
130        if np.isscalar(qy_width):
131            qy_width = np.ones(len(q))*qy_width
132        else:
133            qy_width = np.asarray(qy_width)
134
135        self.q = q.flatten()
136        self.q_calc = slit_extend_q(q, qx_width, qy_width) \
137            if q_calc is None else np.sort(q_calc)
138
139        # Protect against models which are not defined for very low q.  Limit
140        # the smallest q value evaluated (in absolute) to 0.02*min
141        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
142        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
143
144        # Build weight matrix from calculated q values
145        self.weight_matrix = \
146            slit_resolution(self.q_calc, self.q, qx_width, qy_width)
147        self.q_calc = abs(self.q_calc)
148
149    def apply(self, theory):
150        return apply_resolution_matrix(self.weight_matrix, theory)
151
152
153def apply_resolution_matrix(weight_matrix, theory):
154    """
155    Apply the resolution weight matrix to the computed theory function.
156    """
157    #print("apply shapes", theory.shape, weight_matrix.shape)
158    Iq = np.dot(theory[None, :], weight_matrix)
159    #print("result shape",Iq.shape)
160    return Iq.flatten()
161
162
163def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA):
164    r"""
165    Compute the convolution matrix *W* for pinhole resolution 1-D data.
166
167    Each row *W[i]* determines the normalized weight that the corresponding
168    points *q_calc* contribute to the resolution smeared point *q[i]*.  Given
169    *W*, the resolution smearing can be computed using *dot(W,q)*.
170
171    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution
172    function is a broadened triangle, and does not extend over the entire
173    range $(-\infty, +\infty)$.  It is important to impose this limitation
174    since some models fall so steeply that the weighted value in gaussian
175    tails would otherwise dominate the integral.
176
177    *q_calc* must be increasing.  *q_width* must be greater than zero.
178
179    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects
180    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and
181    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14.
182    https://doi.org/10.1107/S0021889894010095.
183    """
184    # The current algorithm is a midpoint rectangle rule.  In the test case,
185    # neither trapezoid nor Simpson's rule improved the accuracy.
186    edges = bin_edges(q_calc)
187    #edges[edges < 0.0] = 0.0 # clip edges below zero
188    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
189    weights = cdf[1:] - cdf[:-1]
190    # Limit q range to +/- 2.5 sigma
191    qhigh = q + nsigma*q_width
192    #qlow = q - nsigma*q_width  # linear limits
193    qlow = q*q/qhigh  # log limits
194    weights[q_calc[:, None] < qlow[None, :]] = 0.
195    weights[q_calc[:, None] > qhigh[None, :]] = 0.
196    weights /= np.sum(weights, axis=0)[None, :]
197    return weights
198
199
200def slit_resolution(q_calc, q, width, height, n_height=30):
201    r"""
202    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
203    $q_\perp$ = *width* and $q_\parallel$ = *height*.  *n_height* is
204    is the number of steps to use in the integration over $q_\parallel$
205    when both $q_\perp$ and $q_\parallel$ are non-zero.
206
207    Each $q$ can have an independent width and height value even though
208    current instruments use the same slit setting for all measured points.
209
210    If slit height is large relative to width, use:
211
212    .. math::
213
214        I_s(q_i) = \frac{1}{\Delta q_\perp}
215            \int_0^{\Delta q_\perp}
216                I\left(\sqrt{q_i^2 + q_\perp^2}\right) \,dq_\perp
217
218    If slit width is large relative to height, use:
219
220    .. math::
221
222        I_s(q_i) = \frac{1}{2 \Delta q_\parallel}
223            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
224                I\left(|q_i + q_\parallel|\right) \,dq_\parallel
225
226    For a mixture of slit width and height use:
227
228    .. math::
229
230        I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp}
231            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
232            \int_0^{\Delta q_\perp}
233                I\left(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}\right)
234                \,dq_\perp dq_\parallel
235
236    **Definition**
237
238    We are using the mid-point integration rule to assign weights to each
239    element of a weight matrix $W$ so that
240
241    .. math::
242
243        I_s(q) = W\,I(q_\text{calc})
244
245    If *q_calc* is at the mid-point, we can infer the bin edges from the
246    pairwise averages of *q_calc*, adding the missing edges before
247    *q_calc[0]* and after *q_calc[-1]*.
248
249    For $q_\parallel = 0$, the smeared value can be computed numerically
250    using the $u$ substitution
251
252    .. math::
253
254        u_j = \sqrt{q_j^2 - q^2}
255
256    This gives
257
258    .. math::
259
260        I_s(q) \approx \sum_j I(u_j) \Delta u_j
261
262    where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the
263    difference between consecutive edges which have been first converted
264    to $u$.  Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds
265    to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so
266
267    .. math::
268
269        W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j
270               = \frac{1}{\Delta q_\perp} \left(
271                    \sqrt{q_{j+1}^2 - q_i^2} - \sqrt{q_j^2 - q_i^2} \right)
272            \ \text{if}\  q_j \in \left[q_i, \sqrt{q_i^2 + q_\perp^2}\right]
273
274    where $I_s(q_i)$ is the theory function being computed and $q_j$ are the
275    mid-points between the calculated values in *q_calc*.  We tweak the
276    edges of the initial and final intervals so that they lie on integration
277    limits.
278
279    (To be precise, the transformed midpoint $u(q_j)$ is not necessarily the
280    midpoint of the edges $u((q_{j-1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$,
281    but it is at least in the interval, so the approximation is going to be
282    a little better than the left or right Riemann sum, and should be
283    good enough for our purposes.)
284
285    For $q_\perp = 0$, the $u$ substitution is simpler:
286
287    .. math::
288
289        u_j = \left|q_j - q\right|
290
291    so
292
293    .. math::
294
295        W_{ij} = \frac{1}{2 \Delta q_\parallel} \Delta u_j
296            = \frac{1}{2 \Delta q_\parallel} (q_{j+1} - q_j)
297            \ \text{if}\ q_j \in
298                \left[q-\Delta q_\parallel, q+\Delta q_\parallel\right]
299
300    However, we need to support cases were $u_j < 0$, which means using
301    $2 (q_{j+1} - q_j)$ when $q_j \in \left[0, q_\parallel-q_i\right]$.
302    This is not an issue for $q_i > q_\parallel$.
303
304    For both $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional
305    integration with
306
307    .. math::
308
309        u_{jk} = \sqrt{q_j^2 - (q + (k\Delta q_\parallel/L))^2}
310            \ \text{for}\ k = -L \ldots L
311
312    for $L$ = *n_height*.  This gives
313
314    .. math::
315
316        W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel}
317            \sum_{k=-L}^L \Delta u_{jk}
318                \left(\frac{\Delta q_\parallel}{2 L + 1}\right)
319
320
321    """
322    #np.set_printoptions(precision=6, linewidth=10000)
323
324    # The current algorithm is a midpoint rectangle rule.
325    q_edges = bin_edges(q_calc) # Note: requires q > 0
326    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
327    weights = np.zeros((len(q), len(q_calc)), 'd')
328
329    #print(q_calc)
330    for i, (qi, w, h) in enumerate(zip(q, width, height)):
331        if w == 0. and h == 0.:
332            # Perfect resolution, so return the theory value directly.
333            # Note: assumes that q is a subset of q_calc.  If qi need not be
334            # in q_calc, then we can do a weighted interpolation by looking
335            # up qi in q_calc, then weighting the result by the relative
336            # distance to the neighbouring points.
337            weights[i, :] = (q_calc == qi)
338        elif h == 0:
339            weights[i, :] = _q_perp_weights(q_edges, qi, w)
340        elif w == 0:
341            in_x = 1.0 * ((q_calc >= qi-h) & (q_calc <= qi+h))
342            abs_x = 1.0*(q_calc < abs(qi - h)) if qi < h else 0.
343            #print(qi - h, qi + h)
344            #print(in_x + abs_x)
345            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
346        else:
347            for k in range(-n_height, n_height+1):
348                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/n_height, w)
349            weights[i, :] /= 2*n_height + 1
350
351    return weights.T
352
353
354def _q_perp_weights(q_edges, qi, w):
355    # Convert bin edges from q to u
356    u_limit = np.sqrt(qi**2 + w**2)
357    u_edges = q_edges**2 - qi**2
358    u_edges[q_edges < abs(qi)] = 0.
359    u_edges[q_edges > u_limit] = u_limit**2 - qi**2
360    weights = np.diff(np.sqrt(u_edges))/w
361    #print("i, qi",i,qi,qi+width)
362    #print(q_calc)
363    #print(weights)
364    return weights
365
366
367def pinhole_extend_q(q, q_width, nsigma=3):
368    """
369    Given *q* and *q_width*, find a set of sampling points *q_calc* so
370    that each point $I(q)$ has sufficient support from the underlying
371    function.
372    """
373    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width)
374    return linear_extrapolation(q, q_min, q_max)
375
376
377def slit_extend_q(q, width, height):
378    """
379    Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
380    that each point I(q) has sufficient support from the underlying
381    function.
382    """
383    q_min, q_max = np.min(q-height), np.max(np.sqrt((q+height)**2 + width**2))
384
385    return geometric_extrapolation(q, q_min, q_max)
386
387
388def bin_edges(x):
389    """
390    Determine bin edges from bin centers, assuming that edges are centered
391    between the bins.
392
393    Note: this uses the arithmetic mean, which may not be appropriate for
394    log-scaled data.
395    """
396    if len(x) < 2 or (np.diff(x) < 0).any():
397        raise ValueError("Expected bins to be an increasing set")
398    edges = np.hstack([
399        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
400        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
401        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
402        ])
403    return edges
404
405
406def interpolate(q, max_step):
407    """
408    Returns *q_calc* with points spaced at most max_step apart.
409    """
410    step = np.diff(q)
411    index = step > max_step
412    if np.any(index):
413        inserts = []
414        for q_i, step_i in zip(q[:-1][index], step[index]):
415            n = np.ceil(step_i/max_step)
416            inserts.extend(q_i + np.arange(1, n)*(step_i/n))
417        # Extend a couple of fringes beyond the end of the data
418        inserts.extend(q[-1] + np.arange(1, 8)*max_step)
419        q_calc = np.sort(np.hstack((q, inserts)))
420    else:
421        q_calc = q
422    return q_calc
423
424
425def linear_extrapolation(q, q_min, q_max):
426    """
427    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as
428    a guide.  Extrapolation below uses about the same size as the first
429    interval.  Extrapolation above uses about the same size as the final
430    interval.
431
432    Note that extrapolated values may be negative.
433    """
434    q = np.sort(q)
435    if q_min + 2*MINIMUM_RESOLUTION < q[0]:
436        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
437        q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
438    else:
439        q_low = []
440    if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
441        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
442        q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
443    else:
444        q_high = []
445    return np.concatenate([q_low, q, q_high])
446
447
449    r"""
450    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the
451    average geometric step size in *q* as the step size.
452
453    if *q_min* is zero or less then *q[0]/10* is used instead.
454
455    *points_per_decade* sets the ratio between consecutive steps such
456    that there will be $n$ points used for every factor of 10 increase
457    in *q*.
458
459    If *points_per_decade* is not given, it will be estimated as follows.
460    Starting at $q_1$ and stepping geometrically by $\Delta q$ to $q_n$
461    in $n$ points gives a geometric average of:
462
463    .. math::
464
465         \log \Delta q = (\log q_n - \log q_1) / (n - 1)
466
467    From this we can compute the number of steps required to extend $q$
468    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
469
470    .. math::
471
472         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
473
474    Substituting:
475
476    .. math::
477
478         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
479            / (\log q_n - \log q_1)
480    """
481    q = np.sort(q)
483        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0]))
484    else:
485        log_delta_q = log(10.) / points_per_decade
486    if q_min < q[0]:
487        if q_min < 0:
488            q_min = q[0]*MINIMUM_ABSOLUTE_Q
489        n_low = log_delta_q * (log(q[0])-log(q_min))
490        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
491    else:
492        q_low = []
493    if q_max > q[-1]:
494        n_high = log_delta_q * (log(q_max)-log(q[-1]))
495        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
496    else:
497        q_high = []
498    return np.concatenate([q_low, q, q_high])
499
500
501############################################################################
502# unit tests
503############################################################################
504
505def eval_form(q, form, pars):
506    """
507    Return the SAS model evaluated at *q*.
508
509    *form* is the SAS model returned from :fun:core.load_model.
510
511    *pars* are the parameter values to use when evaluating.
512    """
513    from sasmodels import direct_model
514    kernel = form.make_kernel([q])
515    theory = direct_model.call_kernel(kernel, pars)
516    kernel.release()
517    return theory
518
519
520def gaussian(q, q0, dq, nsigma=2.5):
521    """
522    Return the truncated Gaussian resolution function.
523
524    *q0* is the center, *dq* is the width and *q* are the points to evaluate.
525    """
526    # Calculate the density of the tails; the resulting gaussian needs to be
527    # scaled by this amount in ordere to integrate to 1.0
528    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2
529    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density)
530
531
532def romberg_slit_1d(q, width, height, form, pars):
533    """
534    Romberg integration for slit resolution.
535
536    This is an adaptive integration technique.  It is called with settings
537    that make it slow to evaluate but give it good accuracy.
538    """
539    from scipy.integrate import romberg  # type: ignore
540
541    par_set = set([p.name for p in form.info.parameters.call_parameters])
542    if any(k not in par_set for k in pars.keys()):
543        extra = set(pars.keys()) - par_set
544        raise ValueError("bad parameters: [%s] not in [%s]"
545                         % (", ".join(sorted(extra)),
546                            ", ".join(sorted(pars.keys()))))
547
548    if np.isscalar(width):
549        width = [width]*len(q)
550    if np.isscalar(height):
551        height = [height]*len(q)
552    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
553    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
554    # If both width and height are defined, then it is too slow to use dblquad.
555    # Instead use trapz on a fixed grid, interpolated into the I(Q) for
556    # the extended Q range.
557    #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
558    q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height))
559    Iq = eval_form(q_calc, form, pars)
560    result = np.empty(len(q))
561    for i, (qi, w, h) in enumerate(zip(q, width, height)):
562        if h == 0.:
563            total = romberg(_int_w, 0, w, args=(qi,),
564                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
565            result[i] = total/w
566        elif w == 0.:
567            total = romberg(_int_h, -h, h, args=(qi,),
568                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
569            result[i] = total/(2*h)
570        else:
571            w_grid = np.linspace(0, w, 21)[None, :]
572            h_grid = np.linspace(-h, h, 23)[:, None]
573            u_sub = sqrt((qi+h_grid)**2 + w_grid**2)
574            f_at_u = np.interp(u_sub, q_calc, Iq)
575            #print(np.trapz(Iu, w_grid, axis=1))
576            total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), h_grid[:, 0])
577            result[i] = total / (2*h*w)
578            # from scipy.integrate import dblquad
579            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
580            #                  args=(qi,))
581            # result[i] = r/(w*2*h)
582
583    # r should be [float, ...], but it is [array([float]), array([float]),...]
584    return result
585
586
587def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5):
588    """
589    Romberg integration for pinhole resolution.
590
591    This is an adaptive integration technique.  It is called with settings
592    that make it slow to evaluate but give it good accuracy.
593    """
594    from scipy.integrate import romberg  # type: ignore
595
596    par_set = set([p.name for p in form.info.parameters.call_parameters])
597    if any(k not in par_set for k in pars.keys()):
598        extra = set(pars.keys()) - par_set
599        raise ValueError("bad parameters: [%s] not in [%s]"
600                         % (", ".join(sorted(extra)),
601                            ", ".join(sorted(pars.keys()))))
602
603    func = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
604    total = [romberg(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
605                     args=(qi, dqi), divmax=100, vec_func=True,
606                     tol=0, rtol=1e-8)
607             for qi, dqi in zip(q, q_width)]
608    return np.asarray(total).flatten()
609
610
611class ResolutionTest(unittest.TestCase):
612    """
613    Test the resolution calculations.
614    """
615
616    def setUp(self):
617        self.x = 0.001*np.arange(1, 11)
618        self.y = self.Iq(self.x)
619
620    def Iq(self, q):
621        "Linear function for resolution unit test"
622        return 12.0 - 1000.0*q
623
624    def test_perfect(self):
625        """
626        Perfect resolution and no smearing.
627        """
628        resolution = Perfect1D(self.x)
629        theory = self.Iq(resolution.q_calc)
630        output = resolution.apply(theory)
631        np.testing.assert_equal(output, self.y)
632
633    def test_slit_zero(self):
634        """
635        Slit smearing with perfect resolution.
636        """
637        resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x)
638        theory = self.Iq(resolution.q_calc)
639        output = resolution.apply(theory)
640        np.testing.assert_equal(output, self.y)
641
642    @unittest.skip("not yet supported")
643    def test_slit_high(self):
644        """
645        Slit smearing with height 0.005
646        """
647        resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x)
648        theory = self.Iq(resolution.q_calc)
649        output = resolution.apply(theory)
651            9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
652            5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
653            ]
655
656    @unittest.skip("not yet supported")
657    def test_slit_both_high(self):
658        """
659        Slit smearing with width < 100*height.
660        """
661        q = np.logspace(-4, -1, 10)
662        resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf)
663        theory = 1000*self.Iq(resolution.q_calc**4)
664        output = resolution.apply(theory)
666            8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
667            5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
668            ]
670
671    @unittest.skip("not yet supported")
672    def test_slit_wide(self):
673        """
674        Slit smearing with width 0.0002
675        """
676        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x)
677        theory = self.Iq(resolution.q_calc)
678        output = resolution.apply(theory)
680            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
681            ]
683
684    @unittest.skip("not yet supported")
685    def test_slit_both_wide(self):
686        """
687        Slit smearing with width > 100*height.
688        """
689        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001,
690                            q_calc=self.x)
691        theory = self.Iq(resolution.q_calc)
692        output = resolution.apply(theory)
694            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
695            ]
697
698    def test_pinhole_zero(self):
699        """
700        Pinhole smearing with perfect resolution
701        """
702        resolution = Pinhole1D(self.x, 0.0*self.x)
703        theory = self.Iq(resolution.q_calc)
704        output = resolution.apply(theory)
705        np.testing.assert_equal(output, self.y)
706
707    # TODO: turn pinhole/slit demos into tests
708
709    @unittest.skip("suppress comparison with old version; pinhole calc changed")
710    def test_pinhole(self):
711        """
712        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001]
713        """
714        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x),
715                               q_calc=self.x)
716        theory = 12.0-1000.0*resolution.q_calc
717        output = resolution.apply(theory)
718        # Note: answer came from output of previous run.  Non-integer
719        # values at ends come from the fact that q_calc does not
720        # extend beyond q, and so the weights don't balance.
722            10.47037734, 9.86925860,
723            9., 8., 7., 6., 5., 4.,
724            3.13074140, 2.52962266,
725            ]
727
728
729class IgorComparisonTest(unittest.TestCase):
730    """
731    Test resolution calculations against those returned by Igor.
732    """
733
734    def setUp(self):
735        self.pars = TEST_PARS_PINHOLE_SPHERE
736        from sasmodels import core
738
739    def _eval_sphere(self, pars, resolution):
740        from sasmodels import direct_model
741        kernel = self.model.make_kernel([resolution.q_calc])
742        theory = direct_model.call_kernel(kernel, pars)
743        result = resolution.apply(theory)
744        kernel.release()
745        return result
746
747    def _compare(self, q, output, answer, tolerance):
749        #idx = abs(err) >= tolerance
750        #problem = zip(q[idx], output[idx], answer[idx], err[idx])
751        #print("\n".join(str(v) for v in problem))
753
754    def test_perfect(self):
755        """
756        Compare sphere model with NIST Igor SANS, no resolution smearing.
757        """
758        pars = TEST_PARS_SLIT_SPHERE
759        data_string = TEST_DATA_SLIT_SPHERE
760
762        q, _, answer, _ = data
763        resolution = Perfect1D(q)
764        output = self._eval_sphere(pars, resolution)
766
767    @unittest.skip("suppress comparison with old version; pinhole calc changed")
768    def test_pinhole(self):
769        """
770        Compare pinhole resolution smearing with NIST Igor SANS
771        """
772        pars = TEST_PARS_PINHOLE_SPHERE
773        data_string = TEST_DATA_PINHOLE_SPHERE
774
776        q, q_width, answer = data
777        resolution = Pinhole1D(q, q_width)
778        output = self._eval_sphere(pars, resolution)
779        # TODO: relative error should be lower
781
782    @unittest.skip("suppress comparison with old version; pinhole calc changed")
783    def test_pinhole_romberg(self):
784        """
785        Compare pinhole resolution smearing with romberg integration result.
786        """
787        pars = TEST_PARS_PINHOLE_SPHERE
788        data_string = TEST_DATA_PINHOLE_SPHERE
790
792        q, q_width, answer = data
793        answer = romberg_pinhole_1d(q, q_width, self.model, pars)
794        ## Getting 0.1% requires 5 sigma and 200 points per fringe
795        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
797        #tol = 0.001
798        ## The default 2.5 sigma and no extra points gets 1%
799        q_calc = None  # type: np.ndarray
800        tol = 0.01
801        resolution = Pinhole1D(q, q_width, q_calc=q_calc)
802        output = self._eval_sphere(pars, resolution)
803        if 0: # debug plot
804            import matplotlib.pyplot as plt  # type: ignore
805            resolution = Perfect1D(q)
806            source = self._eval_sphere(pars, resolution)
807            plt.loglog(q, source, '.')
809            plt.loglog(q, output, '-', hold=True)
810            plt.show()
812
813    def test_slit(self):
814        """
815        Compare slit resolution smearing with NIST Igor SANS
816        """
817        pars = TEST_PARS_SLIT_SPHERE
818        data_string = TEST_DATA_SLIT_SPHERE
819
821        q, delta_qv, _, answer = data
822        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0)
823        output = self._eval_sphere(pars, resolution)
824        # TODO: eliminate Igor test since it is too inaccurate to be useful.
825        # This means we can eliminate the test data as well, and instead
826        # use a generated q vector.
828
829    def test_slit_romberg(self):
830        """
831        Compare slit resolution smearing with romberg integration result.
832        """
833        pars = TEST_PARS_SLIT_SPHERE
834        data_string = TEST_DATA_SLIT_SPHERE
835
837        q, delta_qv, _, answer = data
838        answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars)
840                               delta_qv[0], 0.)
841        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc)
842        output = self._eval_sphere(pars, resolution)
843        # TODO: relative error should be lower
845
846    def test_ellipsoid(self):
847        """
848        Compare romberg integration for ellipsoid model.
849        """
851        pars = {
852            'scale':0.05,
854            'sld':6, 'sld_solvent': 1,
855            }
857        q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
858        width, height = 0.117, 0.
859        resolution = Slit1D(q, qx_width=width, qy_width=height)
860        answer = romberg_slit_1d(q, width, height, form, pars)
861        output = resolution.apply(eval_form(resolution.q_calc, form, pars))
862        # TODO: 10% is too much error; use better algorithm
865
866    #TODO: can sas q spacing be too sparse for the resolution calculation?
867    @unittest.skip("suppress sparse data test; not supported by current code")
868    def test_pinhole_sparse(self):
869        """
870        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
871        """
872        pars = TEST_PARS_PINHOLE_SPHERE
873        data_string = TEST_DATA_PINHOLE_SPHERE
874
876        q, q_width, answer = data[:, ::20] # Take every nth point
877        resolution = Pinhole1D(q, q_width)
878        output = self._eval_sphere(pars, resolution)
880
881
882# pinhole sphere parameters
883TEST_PARS_PINHOLE_SPHERE = {
884    'scale': 1.0, 'background': 0.01,
885    'radius': 60.0, 'sld': 1, 'sld_solvent': 6.3,
886    }
887# Q, dQ, I(Q) calculated by NIST Igor SANS package
888TEST_DATA_PINHOLE_SPHERE = """\
8890.001278 0.0002847 2538.41176383
8900.001562 0.0002905 2536.91820405
8910.001846 0.0002956 2535.13182479
8920.002130 0.0003017 2533.06217813
8930.002414 0.0003087 2530.70378586
8940.002698 0.0003165 2528.05024192
8950.002982 0.0003249 2525.10408349
8960.003266 0.0003340 2521.86667499
8970.003550 0.0003437 2518.33907750
8980.003834 0.0003539 2514.52246995
8990.004118 0.0003646 2510.41798319
9000.004402 0.0003757 2506.02690988
9010.004686 0.0003872 2501.35067884
9020.004970 0.0003990 2496.38678318
9030.005253 0.0004112 2491.16237596
9040.005537 0.0004237 2485.63911673
9050.005821 0.0004365 2479.83657083
9060.006105 0.0004495 2473.75676948
9070.006389 0.0004628 2467.40145990
9080.006673 0.0004762 2460.77293372
9090.006957 0.0004899 2453.86724627
9100.007241 0.0005037 2446.69623838
9110.007525 0.0005177 2439.25775219
9120.007809 0.0005318 2431.55421398
9130.008093 0.0005461 2423.58785521
9140.008377 0.0005605 2415.36158137
9150.008661 0.0005750 2406.87009473
9160.008945 0.0005896 2398.12841186
9170.009229 0.0006044 2389.13360806
9180.009513 0.0006192 2379.88958042
9190.009797 0.0006341 2370.39776774
9200.010080 0.0006491 2360.69528793
9210.010360 0.0006641 2350.85169027
9220.010650 0.0006793 2340.42023633
9230.010930 0.0006945 2330.11206013
9240.011220 0.0007097 2319.20109972
9250.011500 0.0007251 2308.43503981
9260.011780 0.0007404 2297.44820179
9270.012070 0.0007558 2285.83853677
9280.012350 0.0007713 2274.41290746
9290.012640 0.0007868 2262.36219581
9300.012920 0.0008024 2250.51169731
9310.013200 0.0008180 2238.45596231
9320.013490 0.0008336 2225.76495666
9330.013770 0.0008493 2213.29618391
9340.014060 0.0008650 2200.19110751
9350.014340 0.0008807 2187.34050325
9360.014620 0.0008965 2174.30529864
9370.014910 0.0009123 2160.61632548
9380.015190 0.0009281 2147.21038112
9390.015470 0.0009440 2133.62023580
9400.015760 0.0009598 2119.37907426
9410.016040 0.0009757 2105.45234903
9420.016330 0.0009916 2090.86319102
9430.016610 0.0010080 2076.60576032
9440.016890 0.0010240 2062.19214565
9450.017180 0.0010390 2047.10550219
9460.017460 0.0010550 2032.38715621
9470.017740 0.0010710 2017.52560123
9480.018030 0.0010880 2001.99124318
9490.018310 0.0011040 1986.84662060
9500.018600 0.0011200 1971.03389745
9510.018880 0.0011360 1955.61395119
9520.019160 0.0011520 1940.08291563
9530.019450 0.0011680 1923.87672225
9540.019730 0.0011840 1908.10656374
9550.020020 0.0012000 1891.66297192
9560.020300 0.0012160 1875.66789021
9570.020580 0.0012320 1859.56357196
9580.020870 0.0012490 1842.79468290
9590.021150 0.0012650 1826.50064489
9600.021430 0.0012810 1810.11533702
9610.021720 0.0012970 1793.06840882
9620.022000 0.0013130 1776.51153580
9630.022280 0.0013290 1759.87201249
9640.022570 0.0013460 1742.57354412
9650.022850 0.0013620 1725.79397319
9660.023140 0.0013780 1708.35831550
9670.023420 0.0013940 1691.45256069
9680.023700 0.0014110 1674.48561783
9690.023990 0.0014270 1656.86525366
9700.024270 0.0014430 1639.79847285
9710.024550 0.0014590 1622.68887088
9720.024840 0.0014760 1604.96421100
9730.025120 0.0014920 1587.85768129
9740.025410 0.0015080 1569.99297335
9750.025690 0.0015240 1552.84580279
9760.025970 0.0015410 1535.54074115
9770.026260 0.0015570 1517.75249337
9780.026540 0.0015730 1500.40115023
9790.026820 0.0015900 1483.03632237
9800.027110 0.0016060 1465.05942429
9810.027390 0.0016220 1447.67682181
9820.027670 0.0016390 1430.46495191
9830.027960 0.0016550 1412.49232282
9840.028240 0.0016710 1395.13182318
9850.028520 0.0016880 1377.93439837
9860.028810 0.0017040 1359.99528971
9870.029090 0.0017200 1342.67274512
9880.029370 0.0017370 1325.55375609
989"""
990
991# Slit sphere parameters
992TEST_PARS_SLIT_SPHERE = {
993    'scale': 0.01, 'background': 0.01,
994    'radius': 60000, 'sld': 1, 'sld_solvent': 4,
995    }
996# Q dQ I(Q) I_smeared(Q)
997TEST_DATA_SLIT_SPHERE = """\
9982.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
9992.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
10002.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
10013.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
10023.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
10033.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
10044.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
10055.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
10065.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
10076.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
10086.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
10097.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
10107.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
10118.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
10128.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
10139.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
10141.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
10151.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
10161.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
10171.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
10181.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
10191.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
10201.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
10211.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
10221.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
10231.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
10242.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
10252.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
10262.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
10272.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
10282.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
10292.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
10302.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
10312.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
10323.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
10333.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
10343.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
10353.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
10364.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
10374.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
10384.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
10394.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
10405.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
10415.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
10426.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
10436.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
10447.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
10457.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
10468.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
10478.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
10489.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
10499.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
10501.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
10511.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
10521.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
10531.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
10541.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
10551.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
10561.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
10571.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
10581.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
10591.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
10602.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
10612.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
10622.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
10632.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
10642.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
10652.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
1066"""
1067
1068def main():
1069    """
1070    Run tests given is sys.argv.
1071
1072    Returns 0 if success or 1 if any tests fail.
1073    """
1074    import sys
1075    import xmlrunner  # type: ignore
1076
1077    suite = unittest.TestSuite()
1079
1080    runner = xmlrunner.XMLTestRunner(output='logs')
1081    result = runner.run(suite)
1082    return 1 if result.failures or result.errors else 0
1083
1084
1085############################################################################
1086# usage demo
1087############################################################################
1088
1089def _eval_demo_1d(resolution, title):
1090    import sys
1091    from sasmodels import core
1092    from sasmodels import direct_model
1093    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder'
1094
1095    if name == 'cylinder':
1096        pars = {'length':210, 'radius':500, 'background': 0}
1097    elif name == 'teubner_strey':
1098        pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643}
1099    elif name == 'sphere' or name == 'spherepy':
1100        pars = TEST_PARS_SLIT_SPHERE
1101    elif name == 'ellipsoid':
1102        pars = {
1103            'scale':0.05, 'background': 0,
1104            'r_polar':500, 'r_equatorial':15000,
1105            'sld':6, 'sld_solvent': 1,
1106            }
1107    else:
1108        pars = {}
1110    model = core.build_model(model_info)
1111
1112    kernel = model.make_kernel([resolution.q_calc])
1113    theory = direct_model.call_kernel(kernel, pars)
1114    Iq = resolution.apply(theory)
1115
1116    if isinstance(resolution, Slit1D):
1117        width, height = resolution.qx_width, resolution.qy_width
1118        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
1119    else:
1120        dq = resolution.q_width
1121        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars)
1122
1123    import matplotlib.pyplot as plt  # type: ignore
1124    plt.loglog(resolution.q_calc, theory, label='unsmeared')
1125    plt.loglog(resolution.q, Iq, label='smeared', hold=True)
1126    plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True)
1127    plt.legend()
1128    plt.title(title)
1129    plt.xlabel("Q (1/Ang)")
1130    plt.ylabel("I(Q) (1/cm)")
1131
1132def demo_pinhole_1d():
1133    """
1134    Show example of pinhole smearing.
1135    """
1136    q = np.logspace(-4, np.log10(0.2), 400)
1137    q_width = 0.1*q
1138    resolution = Pinhole1D(q, q_width)
1139    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
1140
1141def demo_slit_1d():
1142    """
1143    Show example of slit smearing.
1144    """
1145    q = np.logspace(-4, np.log10(0.2), 100)
1146    w = h = 0.
1147    #w = 0.000000277790
1148    w = 0.0277790
1149    #h = 0.00277790
1150    #h = 0.0277790
1151    resolution = Slit1D(q, w, h)
1152    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
1153
1154def demo():
1155    """
1156    Run the resolution demos.
1157    """
1158    import matplotlib.pyplot as plt  # type: ignore
1159    plt.subplot(121)
1160    demo_pinhole_1d()
1161    #plt.yscale('linear')
1162    plt.subplot(122)
1163    demo_slit_1d()
1164    #plt.yscale('linear')
1165    plt.show()
1166
1167
1168if __name__ == "__main__":
1169    #demo()
1170    main()
Note: See TracBrowser for help on using the repository browser.