source: sasmodels/sasmodels/resolution.py @ 2d81cfe

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

lint

  • Property mode set to 100644
File size: 39.7 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
22
23class Resolution(object):
24    """
25    Abstract base class defining a 1D resolution function.
26
27    *q* is the set of q values at which the data is measured.
28
29    *q_calc* is the set of q values at which the theory needs to be evaluated.
30    This may extend and interpolate the q values.
31
32    *apply* is the method to call with I(q_calc) to compute the resolution
33    smeared theory I(q).
34    """
35    q = None  # type: np.ndarray
36    q_calc = None  # type: np.ndarray
37    def apply(self, theory):
38        """
39        Smear *theory* by the resolution function, returning *Iq*.
40        """
41        raise NotImplementedError("Subclass does not define the apply function")
42
43
44class Perfect1D(Resolution):
45    """
46    Resolution function to use when there is no actual resolution smearing
47    to be applied.  It has the same interface as the other resolution
48    functions, but returns the identity function.
49    """
50    def __init__(self, q):
51        self.q_calc = self.q = q
52
53    def apply(self, theory):
54        return theory
55
56
57class Pinhole1D(Resolution):
58    r"""
59    Pinhole aperture with q-dependent gaussian resolution.
60
61    *q* points at which the data is measured.
62
63    *q_width* gaussian 1-sigma resolution at each data point.
64
65    *q_calc* is the list of points to calculate, or None if this should
66    be estimated from the *q* and *q_width*.
67    """
68    def __init__(self, q, q_width, q_calc=None, nsigma=3):
69        #*min_step* is the minimum point spacing to use when computing the
70        #underlying model.  It should be on the order of
71        #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes
72        #are computed with sufficient density to avoid aliasing effects.
73
74        # Protect against calls with q_width=0.  The extend_q function will
75        # not extend the q if q_width is 0, but q_width must be non-zero when
76        # constructing the weight matrix to avoid division by zero errors.
77        # In practice this should never be needed, since resolution should
78        # default to Perfect1D if the pinhole geometry is not defined.
79        self.q, self.q_width = q, q_width
80        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
81                       if q_calc is None else np.sort(q_calc))
82
83        # Protect against models which are not defined for very low q.  Limit
84        # the smallest q value evaluated (in absolute) to 0.02*min
85        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
86        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
87
88        # Build weight matrix from calculated q values
89        self.weight_matrix = pinhole_resolution(
90            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION))
91        self.q_calc = abs(self.q_calc)
92
93    def apply(self, theory):
94        return apply_resolution_matrix(self.weight_matrix, theory)
95
96
97class Slit1D(Resolution):
98    """
99    Slit aperture with resolution function.
100
101    *q* points at which the data is measured.
102
103    *dqx* slit width in qx
104
105    *dqy* slit height in qy
106
107    *q_calc* is the list of points to calculate, or None if this should
108    be estimated from the *q* and *q_width*.
109
110    The *weight_matrix* is computed by :func:`slit1d_resolution`
111    """
112    def __init__(self, q, qx_width, qy_width=0., q_calc=None):
113        # Remember what width/dqy was used even though we won't need them
114        # after the weight matrix is constructed
115        self.qx_width, self.qy_width = qx_width, qy_width
116
117        # Allow independent resolution on each point even though it is not
118        # needed in practice.
119        if np.isscalar(qx_width):
120            qx_width = np.ones(len(q))*qx_width
121        else:
122            qx_width = np.asarray(qx_width)
123        if np.isscalar(qy_width):
124            qy_width = np.ones(len(q))*qy_width
125        else:
126            qy_width = np.asarray(qy_width)
127
128        self.q = q.flatten()
129        self.q_calc = slit_extend_q(q, qx_width, qy_width) \
130            if q_calc is None else np.sort(q_calc)
131
132        # Protect against models which are not defined for very low q.  Limit
133        # the smallest q value evaluated (in absolute) to 0.02*min
134        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
135        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
136
137        # Build weight matrix from calculated q values
138        self.weight_matrix = \
139            slit_resolution(self.q_calc, self.q, qx_width, qy_width)
140        self.q_calc = abs(self.q_calc)
141
142    def apply(self, theory):
143        return apply_resolution_matrix(self.weight_matrix, theory)
144
145
146def apply_resolution_matrix(weight_matrix, theory):
147    """
148    Apply the resolution weight matrix to the computed theory function.
149    """
150    #print("apply shapes", theory.shape, weight_matrix.shape)
151    Iq = np.dot(theory[None, :], weight_matrix)
152    #print("result shape",Iq.shape)
153    return Iq.flatten()
154
155
156def pinhole_resolution(q_calc, q, q_width):
157    """
158    Compute the convolution matrix *W* for pinhole resolution 1-D data.
159
160    Each row *W[i]* determines the normalized weight that the corresponding
161    points *q_calc* contribute to the resolution smeared point *q[i]*.  Given
162    *W*, the resolution smearing can be computed using *dot(W,q)*.
163
164    *q_calc* must be increasing.  *q_width* must be greater than zero.
165    """
166    # The current algorithm is a midpoint rectangle rule.  In the test case,
167    # neither trapezoid nor Simpson's rule improved the accuracy.
168    edges = bin_edges(q_calc)
169    #edges[edges < 0.0] = 0.0 # clip edges below zero
170    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
171    weights = cdf[1:] - cdf[:-1]
172    weights /= np.sum(weights, axis=0)[None, :]
173    return weights
174
175
176def slit_resolution(q_calc, q, width, height, n_height=30):
177    r"""
178    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
179    $q_\perp$ = *width* and $q_\parallel$ = *height*.  *n_height* is
180    is the number of steps to use in the integration over $q_\parallel$
181    when both $q_\perp$ and $q_\parallel$ are non-zero.
182
183    Each $q$ can have an independent width and height value even though
184    current instruments use the same slit setting for all measured points.
185
186    If slit height is large relative to width, use:
187
188    .. math::
189
190        I_s(q_i) = \frac{1}{\Delta q_\perp}
191            \int_0^{\Delta q_\perp}
192                I\left(\sqrt{q_i^2 + q_\perp^2}\right) \,dq_\perp
193
194    If slit width is large relative to height, use:
195
196    .. math::
197
198        I_s(q_i) = \frac{1}{2 \Delta q_\parallel}
199            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
200                I\left(|q_i + q_\parallel|\right) \,dq_\parallel
201
202    For a mixture of slit width and height use:
203
204    .. math::
205
206        I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp}
207            \int_{-\Delta q_\parallel}^{\Delta q_\parallel}
208            \int_0^{\Delta q_\perp}
209                I\left(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}\right)
210                \,dq_\perp dq_\parallel
211
212    **Definition**
213
214    We are using the mid-point integration rule to assign weights to each
215    element of a weight matrix $W$ so that
216
217    .. math::
218
219        I_s(q) = W\,I(q_\text{calc})
220
221    If *q_calc* is at the mid-point, we can infer the bin edges from the
222    pairwise averages of *q_calc*, adding the missing edges before
223    *q_calc[0]* and after *q_calc[-1]*.
224
225    For $q_\parallel = 0$, the smeared value can be computed numerically
226    using the $u$ substitution
227
228    .. math::
229
230        u_j = \sqrt{q_j^2 - q^2}
231
232    This gives
233
234    .. math::
235
236        I_s(q) \approx \sum_j I(u_j) \Delta u_j
237
238    where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the
239    difference between consecutive edges which have been first converted
240    to $u$.  Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds
241    to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so
242
243    .. math::
244
245        W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j
246               = \frac{1}{\Delta q_\perp} \left(
247                    \sqrt{q_{j+1}^2 - q_i^2} - \sqrt{q_j^2 - q_i^2} \right)
248            \ \text{if}\  q_j \in \left[q_i, \sqrt{q_i^2 + q_\perp^2}\right]
249
250    where $I_s(q_i)$ is the theory function being computed and $q_j$ are the
251    mid-points between the calculated values in *q_calc*.  We tweak the
252    edges of the initial and final intervals so that they lie on integration
253    limits.
254
255    (To be precise, the transformed midpoint $u(q_j)$ is not necessarily the
256    midpoint of the edges $u((q_{j-1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$,
257    but it is at least in the interval, so the approximation is going to be
258    a little better than the left or right Riemann sum, and should be
259    good enough for our purposes.)
260
261    For $q_\perp = 0$, the $u$ substitution is simpler:
262
263    .. math::
264
265        u_j = \left|q_j - q\right|
266
267    so
268
269    .. math::
270
271        W_{ij} = \frac{1}{2 \Delta q_\parallel} \Delta u_j
272            = \frac{1}{2 \Delta q_\parallel} (q_{j+1} - q_j)
273            \ \text{if}\ q_j \in
274                \left[q-\Delta q_\parallel, q+\Delta q_\parallel\right]
275
276    However, we need to support cases were $u_j < 0$, which means using
277    $2 (q_{j+1} - q_j)$ when $q_j \in \left[0, q_\parallel-q_i\right]$.
278    This is not an issue for $q_i > q_\parallel$.
279
280    For both $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional
281    integration with
282
283    .. math::
284
285        u_{jk} = \sqrt{q_j^2 - (q + (k\Delta q_\parallel/L))^2}
286            \ \text{for}\ k = -L \ldots L
287
288    for $L$ = *n_height*.  This gives
289
290    .. math::
291
292        W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel}
293            \sum_{k=-L}^L \Delta u_{jk}
294                \left(\frac{\Delta q_\parallel}{2 L + 1}\right)
295
296
297    """
298    #np.set_printoptions(precision=6, linewidth=10000)
299
300    # The current algorithm is a midpoint rectangle rule.
301    q_edges = bin_edges(q_calc) # Note: requires q > 0
302    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
303    weights = np.zeros((len(q), len(q_calc)), 'd')
304
305    #print(q_calc)
306    for i, (qi, w, h) in enumerate(zip(q, width, height)):
307        if w == 0. and h == 0.:
308            # Perfect resolution, so return the theory value directly.
309            # Note: assumes that q is a subset of q_calc.  If qi need not be
310            # in q_calc, then we can do a weighted interpolation by looking
311            # up qi in q_calc, then weighting the result by the relative
312            # distance to the neighbouring points.
313            weights[i, :] = (q_calc == qi)
314        elif h == 0:
315            weights[i, :] = _q_perp_weights(q_edges, qi, w)
316        elif w == 0:
317            in_x = 1.0 * ((q_calc >= qi-h) & (q_calc <= qi+h))
318            abs_x = 1.0*(q_calc < abs(qi - h)) if qi < h else 0.
319            #print(qi - h, qi + h)
320            #print(in_x + abs_x)
321            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
322        else:
323            for k in range(-n_height, n_height+1):
324                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/n_height, w)
325            weights[i, :] /= 2*n_height + 1
326
327    return weights.T
328
329
330def _q_perp_weights(q_edges, qi, w):
331    # Convert bin edges from q to u
332    u_limit = np.sqrt(qi**2 + w**2)
333    u_edges = q_edges**2 - qi**2
334    u_edges[q_edges < abs(qi)] = 0.
335    u_edges[q_edges > u_limit] = u_limit**2 - qi**2
336    weights = np.diff(np.sqrt(u_edges))/w
337    #print("i, qi",i,qi,qi+width)
338    #print(q_calc)
339    #print(weights)
340    return weights
341
342
343def pinhole_extend_q(q, q_width, nsigma=3):
344    """
345    Given *q* and *q_width*, find a set of sampling points *q_calc* so
346    that each point $I(q)$ has sufficient support from the underlying
347    function.
348    """
349    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width)
350    return linear_extrapolation(q, q_min, q_max)
351
352
353def slit_extend_q(q, width, height):
354    """
355    Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
356    that each point I(q) has sufficient support from the underlying
357    function.
358    """
359    q_min, q_max = np.min(q-height), np.max(np.sqrt((q+height)**2 + width**2))
360
361    return geometric_extrapolation(q, q_min, q_max)
362
363
364def bin_edges(x):
365    """
366    Determine bin edges from bin centers, assuming that edges are centered
367    between the bins.
368
369    Note: this uses the arithmetic mean, which may not be appropriate for
370    log-scaled data.
371    """
372    if len(x) < 2 or (np.diff(x) < 0).any():
373        raise ValueError("Expected bins to be an increasing set")
374    edges = np.hstack([
375        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
376        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
377        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
378        ])
379    return edges
380
381
382def interpolate(q, max_step):
383    """
384    Returns *q_calc* with points spaced at most max_step apart.
385    """
386    step = np.diff(q)
387    index = step > max_step
388    if np.any(index):
389        inserts = []
390        for q_i, step_i in zip(q[:-1][index], step[index]):
391            n = np.ceil(step_i/max_step)
392            inserts.extend(q_i + np.arange(1, n)*(step_i/n))
393        # Extend a couple of fringes beyond the end of the data
394        inserts.extend(q[-1] + np.arange(1, 8)*max_step)
395        q_calc = np.sort(np.hstack((q, inserts)))
396    else:
397        q_calc = q
398    return q_calc
399
400
401def linear_extrapolation(q, q_min, q_max):
402    """
403    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as
404    a guide.  Extrapolation below uses about the same size as the first
405    interval.  Extrapolation above uses about the same size as the final
406    interval.
407
408    Note that extrapolated values may be negative.
409    """
410    q = np.sort(q)
411    if q_min + 2*MINIMUM_RESOLUTION < q[0]:
412        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
413        q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
414    else:
415        q_low = []
416    if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
417        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
418        q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
419    else:
420        q_high = []
421    return np.concatenate([q_low, q, q_high])
422
423
424def geometric_extrapolation(q, q_min, q_max, points_per_decade=None):
425    r"""
426    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the
427    average geometric step size in *q* as the step size.
428
429    if *q_min* is zero or less then *q[0]/10* is used instead.
430
431    *points_per_decade* sets the ratio between consecutive steps such
432    that there will be $n$ points used for every factor of 10 increase
433    in *q*.
434
435    If *points_per_decade* is not given, it will be estimated as follows.
436    Starting at $q_1$ and stepping geometrically by $\Delta q$ to $q_n$
437    in $n$ points gives a geometric average of:
438
439    .. math::
440
441         \log \Delta q = (\log q_n - \log q_1) / (n - 1)
442
443    From this we can compute the number of steps required to extend $q$
444    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
445
446    .. math::
447
448         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
449
450    Substituting:
451
452    .. math::
453
454         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
455            / (\log q_n - \log q_1)
456    """
457    q = np.sort(q)
458    if points_per_decade is None:
459        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0]))
460    else:
461        log_delta_q = log(10.) / points_per_decade
462    if q_min < q[0]:
463        if q_min < 0:
464            q_min = q[0]*MINIMUM_ABSOLUTE_Q
465        n_low = log_delta_q * (log(q[0])-log(q_min))
466        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
467    else:
468        q_low = []
469    if q_max > q[-1]:
470        n_high = log_delta_q * (log(q_max)-log(q[-1]))
471        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
472    else:
473        q_high = []
474    return np.concatenate([q_low, q, q_high])
475
476
477############################################################################
478# unit tests
479############################################################################
480
481def eval_form(q, form, pars):
482    """
483    Return the SAS model evaluated at *q*.
484
485    *form* is the SAS model returned from :fun:`core.load_model`.
486
487    *pars* are the parameter values to use when evaluating.
488    """
489    from sasmodels import direct_model
490    kernel = form.make_kernel([q])
491    theory = direct_model.call_kernel(kernel, pars)
492    kernel.release()
493    return theory
494
495
496def gaussian(q, q0, dq):
497    """
498    Return the Gaussian resolution function.
499
500    *q0* is the center, *dq* is the width and *q* are the points to evaluate.
501    """
502    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)
503
504
505def romberg_slit_1d(q, width, height, form, pars):
506    """
507    Romberg integration for slit resolution.
508
509    This is an adaptive integration technique.  It is called with settings
510    that make it slow to evaluate but give it good accuracy.
511    """
512    from scipy.integrate import romberg  # type: ignore
513
514    par_set = set([p.name for p in form.info.parameters.call_parameters])
515    if any(k not in par_set for k in pars.keys()):
516        extra = set(pars.keys()) - par_set
517        raise ValueError("bad parameters: [%s] not in [%s]"
518                         % (", ".join(sorted(extra)),
519                            ", ".join(sorted(pars.keys()))))
520
521    if np.isscalar(width):
522        width = [width]*len(q)
523    if np.isscalar(height):
524        height = [height]*len(q)
525    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
526    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
527    # If both width and height are defined, then it is too slow to use dblquad.
528    # Instead use trapz on a fixed grid, interpolated into the I(Q) for
529    # the extended Q range.
530    #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
531    q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height))
532    Iq = eval_form(q_calc, form, pars)
533    result = np.empty(len(q))
534    for i, (qi, w, h) in enumerate(zip(q, width, height)):
535        if h == 0.:
536            total = romberg(_int_w, 0, w, args=(qi,),
537                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
538            result[i] = total/w
539        elif w == 0.:
540            total = romberg(_int_h, -h, h, args=(qi,),
541                            divmax=100, vec_func=True, tol=0, rtol=1e-8)
542            result[i] = total/(2*h)
543        else:
544            w_grid = np.linspace(0, w, 21)[None, :]
545            h_grid = np.linspace(-h, h, 23)[:, None]
546            u_sub = sqrt((qi+h_grid)**2 + w_grid**2)
547            f_at_u = np.interp(u_sub, q_calc, Iq)
548            #print(np.trapz(Iu, w_grid, axis=1))
549            total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), h_grid[:, 0])
550            result[i] = total / (2*h*w)
551            # from scipy.integrate import dblquad
552            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
553            #                  args=(qi,))
554            # result[i] = r/(w*2*h)
555
556    # r should be [float, ...], but it is [array([float]), array([float]),...]
557    return result
558
559
560def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5):
561    """
562    Romberg integration for pinhole resolution.
563
564    This is an adaptive integration technique.  It is called with settings
565    that make it slow to evaluate but give it good accuracy.
566    """
567    from scipy.integrate import romberg  # type: ignore
568
569    par_set = set([p.name for p in form.info.parameters.call_parameters])
570    if any(k not in par_set for k in pars.keys()):
571        extra = set(pars.keys()) - par_set
572        raise ValueError("bad parameters: [%s] not in [%s]"
573                         % (", ".join(sorted(extra)),
574                            ", ".join(sorted(pars.keys()))))
575
576    func = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
577    total = [romberg(func, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
578                     args=(qi, dqi), divmax=100, vec_func=True,
579                     tol=0, rtol=1e-8)
580             for qi, dqi in zip(q, q_width)]
581    return np.asarray(total).flatten()
582
583
584class ResolutionTest(unittest.TestCase):
585    """
586    Test the resolution calculations.
587    """
588
589    def setUp(self):
590        self.x = 0.001*np.arange(1, 11)
591        self.y = self.Iq(self.x)
592
593    def Iq(self, q):
594        "Linear function for resolution unit test"
595        return 12.0 - 1000.0*q
596
597    def test_perfect(self):
598        """
599        Perfect resolution and no smearing.
600        """
601        resolution = Perfect1D(self.x)
602        theory = self.Iq(resolution.q_calc)
603        output = resolution.apply(theory)
604        np.testing.assert_equal(output, self.y)
605
606    def test_slit_zero(self):
607        """
608        Slit smearing with perfect resolution.
609        """
610        resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x)
611        theory = self.Iq(resolution.q_calc)
612        output = resolution.apply(theory)
613        np.testing.assert_equal(output, self.y)
614
615    @unittest.skip("not yet supported")
616    def test_slit_high(self):
617        """
618        Slit smearing with height 0.005
619        """
620        resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x)
621        theory = self.Iq(resolution.q_calc)
622        output = resolution.apply(theory)
623        answer = [
624            9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
625            5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
626            ]
627        np.testing.assert_allclose(output, answer, atol=1e-4)
628
629    @unittest.skip("not yet supported")
630    def test_slit_both_high(self):
631        """
632        Slit smearing with width < 100*height.
633        """
634        q = np.logspace(-4, -1, 10)
635        resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf)
636        theory = 1000*self.Iq(resolution.q_calc**4)
637        output = resolution.apply(theory)
638        answer = [
639            8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
640            5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
641            ]
642        np.testing.assert_allclose(output, answer, atol=1e-4)
643
644    @unittest.skip("not yet supported")
645    def test_slit_wide(self):
646        """
647        Slit smearing with width 0.0002
648        """
649        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x)
650        theory = self.Iq(resolution.q_calc)
651        output = resolution.apply(theory)
652        answer = [
653            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
654            ]
655        np.testing.assert_allclose(output, answer, atol=1e-4)
656
657    @unittest.skip("not yet supported")
658    def test_slit_both_wide(self):
659        """
660        Slit smearing with width > 100*height.
661        """
662        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001,
663                            q_calc=self.x)
664        theory = self.Iq(resolution.q_calc)
665        output = resolution.apply(theory)
666        answer = [
667            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
668            ]
669        np.testing.assert_allclose(output, answer, atol=1e-4)
670
671    def test_pinhole_zero(self):
672        """
673        Pinhole smearing with perfect resolution
674        """
675        resolution = Pinhole1D(self.x, 0.0*self.x)
676        theory = self.Iq(resolution.q_calc)
677        output = resolution.apply(theory)
678        np.testing.assert_equal(output, self.y)
679
680    def test_pinhole(self):
681        """
682        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001]
683        """
684        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x),
685                               q_calc=self.x)
686        theory = 12.0-1000.0*resolution.q_calc
687        output = resolution.apply(theory)
688        answer = [
689            10.44785079, 9.84991299, 8.98101708,
690            7.99906585, 6.99998311, 6.00001689,
691            5.00093415, 4.01898292, 3.15008701, 2.55214921,
692            ]
693        np.testing.assert_allclose(output, answer, atol=1e-8)
694
695
696class IgorComparisonTest(unittest.TestCase):
697    """
698    Test resolution calculations against those returned by Igor.
699    """
700
701    def setUp(self):
702        self.pars = TEST_PARS_PINHOLE_SPHERE
703        from sasmodels import core
704        self.model = core.load_model("sphere", dtype='double')
705
706    def _eval_sphere(self, pars, resolution):
707        from sasmodels import direct_model
708        kernel = self.model.make_kernel([resolution.q_calc])
709        theory = direct_model.call_kernel(kernel, pars)
710        result = resolution.apply(theory)
711        kernel.release()
712        return result
713
714    def _compare(self, q, output, answer, tolerance):
715        #err = (output - answer)/answer
716        #idx = abs(err) >= tolerance
717        #problem = zip(q[idx], output[idx], answer[idx], err[idx])
718        #print("\n".join(str(v) for v in problem))
719        np.testing.assert_allclose(output, answer, rtol=tolerance)
720
721    def test_perfect(self):
722        """
723        Compare sphere model with NIST Igor SANS, no resolution smearing.
724        """
725        pars = TEST_PARS_SLIT_SPHERE
726        data_string = TEST_DATA_SLIT_SPHERE
727
728        data = np.loadtxt(data_string.split('\n')).T
729        q, _, answer, _ = data
730        resolution = Perfect1D(q)
731        output = self._eval_sphere(pars, resolution)
732        self._compare(q, output, answer, 1e-6)
733
734    def test_pinhole(self):
735        """
736        Compare pinhole resolution smearing with NIST Igor SANS
737        """
738        pars = TEST_PARS_PINHOLE_SPHERE
739        data_string = TEST_DATA_PINHOLE_SPHERE
740
741        data = np.loadtxt(data_string.split('\n')).T
742        q, q_width, answer = data
743        resolution = Pinhole1D(q, q_width)
744        output = self._eval_sphere(pars, resolution)
745        # TODO: relative error should be lower
746        self._compare(q, output, answer, 3e-4)
747
748    def test_pinhole_romberg(self):
749        """
750        Compare pinhole resolution smearing with romberg integration result.
751        """
752        pars = TEST_PARS_PINHOLE_SPHERE
753        data_string = TEST_DATA_PINHOLE_SPHERE
754        pars['radius'] *= 5
755
756        data = np.loadtxt(data_string.split('\n')).T
757        q, q_width, answer = data
758        answer = romberg_pinhole_1d(q, q_width, self.model, pars)
759        ## Getting 0.1% requires 5 sigma and 200 points per fringe
760        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
761        #                     2*np.pi/pars['radius']/200)
762        #tol = 0.001
763        ## The default 3 sigma and no extra points gets 1%
764        q_calc = None  # type: np.ndarray
765        tol = 0.01
766        resolution = Pinhole1D(q, q_width, q_calc=q_calc)
767        output = self._eval_sphere(pars, resolution)
768        if 0: # debug plot
769            import matplotlib.pyplot as plt  # type: ignore
770            resolution = Perfect1D(q)
771            source = self._eval_sphere(pars, resolution)
772            plt.loglog(q, source, '.')
773            plt.loglog(q, answer, '-', hold=True)
774            plt.loglog(q, output, '-', hold=True)
775            plt.show()
776        self._compare(q, output, answer, tol)
777
778    def test_slit(self):
779        """
780        Compare slit resolution smearing with NIST Igor SANS
781        """
782        pars = TEST_PARS_SLIT_SPHERE
783        data_string = TEST_DATA_SLIT_SPHERE
784
785        data = np.loadtxt(data_string.split('\n')).T
786        q, delta_qv, _, answer = data
787        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0)
788        output = self._eval_sphere(pars, resolution)
789        # TODO: eliminate Igor test since it is too inaccurate to be useful.
790        # This means we can eliminate the test data as well, and instead
791        # use a generated q vector.
792        self._compare(q, output, answer, 0.5)
793
794    def test_slit_romberg(self):
795        """
796        Compare slit resolution smearing with romberg integration result.
797        """
798        pars = TEST_PARS_SLIT_SPHERE
799        data_string = TEST_DATA_SLIT_SPHERE
800
801        data = np.loadtxt(data_string.split('\n')).T
802        q, delta_qv, _, answer = data
803        answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars)
804        q_calc = slit_extend_q(interpolate(q, 2*np.pi/pars['radius']/20),
805                               delta_qv[0], 0.)
806        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc)
807        output = self._eval_sphere(pars, resolution)
808        # TODO: relative error should be lower
809        self._compare(q, output, answer, 0.025)
810
811    def test_ellipsoid(self):
812        """
813        Compare romberg integration for ellipsoid model.
814        """
815        from .core import load_model
816        pars = {
817            'scale':0.05,
818            'radius_polar':500, 'radius_equatorial':15000,
819            'sld':6, 'sld_solvent': 1,
820            }
821        form = load_model('ellipsoid', dtype='double')
822        q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
823        width, height = 0.117, 0.
824        resolution = Slit1D(q, qx_width=width, qy_width=height)
825        answer = romberg_slit_1d(q, width, height, form, pars)
826        output = resolution.apply(eval_form(resolution.q_calc, form, pars))
827        # TODO: 10% is too much error; use better algorithm
828        #print(np.max(abs(answer-output)/answer))
829        self._compare(q, output, answer, 0.1)
830
831    #TODO: can sas q spacing be too sparse for the resolution calculation?
832    @unittest.skip("suppress sparse data test; not supported by current code")
833    def test_pinhole_sparse(self):
834        """
835        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
836        """
837        pars = TEST_PARS_PINHOLE_SPHERE
838        data_string = TEST_DATA_PINHOLE_SPHERE
839
840        data = np.loadtxt(data_string.split('\n')).T
841        q, q_width, answer = data[:, ::20] # Take every nth point
842        resolution = Pinhole1D(q, q_width)
843        output = self._eval_sphere(pars, resolution)
844        self._compare(q, output, answer, 1e-6)
845
846
847# pinhole sphere parameters
848TEST_PARS_PINHOLE_SPHERE = {
849    'scale': 1.0, 'background': 0.01,
850    'radius': 60.0, 'sld': 1, 'sld_solvent': 6.3,
851    }
852# Q, dQ, I(Q) calculated by NIST Igor SANS package
853TEST_DATA_PINHOLE_SPHERE = """\
8540.001278 0.0002847 2538.41176383
8550.001562 0.0002905 2536.91820405
8560.001846 0.0002956 2535.13182479
8570.002130 0.0003017 2533.06217813
8580.002414 0.0003087 2530.70378586
8590.002698 0.0003165 2528.05024192
8600.002982 0.0003249 2525.10408349
8610.003266 0.0003340 2521.86667499
8620.003550 0.0003437 2518.33907750
8630.003834 0.0003539 2514.52246995
8640.004118 0.0003646 2510.41798319
8650.004402 0.0003757 2506.02690988
8660.004686 0.0003872 2501.35067884
8670.004970 0.0003990 2496.38678318
8680.005253 0.0004112 2491.16237596
8690.005537 0.0004237 2485.63911673
8700.005821 0.0004365 2479.83657083
8710.006105 0.0004495 2473.75676948
8720.006389 0.0004628 2467.40145990
8730.006673 0.0004762 2460.77293372
8740.006957 0.0004899 2453.86724627
8750.007241 0.0005037 2446.69623838
8760.007525 0.0005177 2439.25775219
8770.007809 0.0005318 2431.55421398
8780.008093 0.0005461 2423.58785521
8790.008377 0.0005605 2415.36158137
8800.008661 0.0005750 2406.87009473
8810.008945 0.0005896 2398.12841186
8820.009229 0.0006044 2389.13360806
8830.009513 0.0006192 2379.88958042
8840.009797 0.0006341 2370.39776774
8850.010080 0.0006491 2360.69528793
8860.010360 0.0006641 2350.85169027
8870.010650 0.0006793 2340.42023633
8880.010930 0.0006945 2330.11206013
8890.011220 0.0007097 2319.20109972
8900.011500 0.0007251 2308.43503981
8910.011780 0.0007404 2297.44820179
8920.012070 0.0007558 2285.83853677
8930.012350 0.0007713 2274.41290746
8940.012640 0.0007868 2262.36219581
8950.012920 0.0008024 2250.51169731
8960.013200 0.0008180 2238.45596231
8970.013490 0.0008336 2225.76495666
8980.013770 0.0008493 2213.29618391
8990.014060 0.0008650 2200.19110751
9000.014340 0.0008807 2187.34050325
9010.014620 0.0008965 2174.30529864
9020.014910 0.0009123 2160.61632548
9030.015190 0.0009281 2147.21038112
9040.015470 0.0009440 2133.62023580
9050.015760 0.0009598 2119.37907426
9060.016040 0.0009757 2105.45234903
9070.016330 0.0009916 2090.86319102
9080.016610 0.0010080 2076.60576032
9090.016890 0.0010240 2062.19214565
9100.017180 0.0010390 2047.10550219
9110.017460 0.0010550 2032.38715621
9120.017740 0.0010710 2017.52560123
9130.018030 0.0010880 2001.99124318
9140.018310 0.0011040 1986.84662060
9150.018600 0.0011200 1971.03389745
9160.018880 0.0011360 1955.61395119
9170.019160 0.0011520 1940.08291563
9180.019450 0.0011680 1923.87672225
9190.019730 0.0011840 1908.10656374
9200.020020 0.0012000 1891.66297192
9210.020300 0.0012160 1875.66789021
9220.020580 0.0012320 1859.56357196
9230.020870 0.0012490 1842.79468290
9240.021150 0.0012650 1826.50064489
9250.021430 0.0012810 1810.11533702
9260.021720 0.0012970 1793.06840882
9270.022000 0.0013130 1776.51153580
9280.022280 0.0013290 1759.87201249
9290.022570 0.0013460 1742.57354412
9300.022850 0.0013620 1725.79397319
9310.023140 0.0013780 1708.35831550
9320.023420 0.0013940 1691.45256069
9330.023700 0.0014110 1674.48561783
9340.023990 0.0014270 1656.86525366
9350.024270 0.0014430 1639.79847285
9360.024550 0.0014590 1622.68887088
9370.024840 0.0014760 1604.96421100
9380.025120 0.0014920 1587.85768129
9390.025410 0.0015080 1569.99297335
9400.025690 0.0015240 1552.84580279
9410.025970 0.0015410 1535.54074115
9420.026260 0.0015570 1517.75249337
9430.026540 0.0015730 1500.40115023
9440.026820 0.0015900 1483.03632237
9450.027110 0.0016060 1465.05942429
9460.027390 0.0016220 1447.67682181
9470.027670 0.0016390 1430.46495191
9480.027960 0.0016550 1412.49232282
9490.028240 0.0016710 1395.13182318
9500.028520 0.0016880 1377.93439837
9510.028810 0.0017040 1359.99528971
9520.029090 0.0017200 1342.67274512
9530.029370 0.0017370 1325.55375609
954"""
955
956# Slit sphere parameters
957TEST_PARS_SLIT_SPHERE = {
958    'scale': 0.01, 'background': 0.01,
959    'radius': 60000, 'sld': 1, 'sld_solvent': 4,
960    }
961# Q dQ I(Q) I_smeared(Q)
962TEST_DATA_SLIT_SPHERE = """\
9632.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
9642.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
9652.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
9663.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
9673.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
9683.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
9694.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
9705.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
9715.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
9726.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
9736.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
9747.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
9757.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
9768.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
9778.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
9789.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
9791.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
9801.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
9811.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
9821.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
9831.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
9841.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
9851.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
9861.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
9871.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
9881.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
9892.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
9902.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
9912.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
9922.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
9932.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
9942.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
9952.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
9962.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
9973.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
9983.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
9993.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
10003.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
10014.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
10024.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
10034.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
10044.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
10055.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
10065.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
10076.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
10086.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
10097.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
10107.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
10118.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
10128.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
10139.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
10149.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
10151.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
10161.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
10171.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
10181.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
10191.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
10201.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
10211.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
10221.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
10231.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
10241.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
10252.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
10262.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
10272.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
10282.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
10292.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
10302.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
1031"""
1032
1033def main():
1034    """
1035    Run tests given is sys.argv.
1036
1037    Returns 0 if success or 1 if any tests fail.
1038    """
1039    import sys
1040    import xmlrunner  # type: ignore
1041
1042    suite = unittest.TestSuite()
1043    suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__]))
1044
1045    runner = xmlrunner.XMLTestRunner(output='logs')
1046    result = runner.run(suite)
1047    return 1 if result.failures or result.errors else 0
1048
1049
1050############################################################################
1051# usage demo
1052############################################################################
1053
1054def _eval_demo_1d(resolution, title):
1055    import sys
1056    from sasmodels import core
1057    from sasmodels import direct_model
1058    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder'
1059
1060    if name == 'cylinder':
1061        pars = {'length':210, 'radius':500, 'background': 0}
1062    elif name == 'teubner_strey':
1063        pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643}
1064    elif name == 'sphere' or name == 'spherepy':
1065        pars = TEST_PARS_SLIT_SPHERE
1066    elif name == 'ellipsoid':
1067        pars = {
1068            'scale':0.05, 'background': 0,
1069            'r_polar':500, 'r_equatorial':15000,
1070            'sld':6, 'sld_solvent': 1,
1071            }
1072    else:
1073        pars = {}
1074    model_info = core.load_model_info(name)
1075    model = core.build_model(model_info)
1076
1077    kernel = model.make_kernel([resolution.q_calc])
1078    theory = direct_model.call_kernel(kernel, pars)
1079    Iq = resolution.apply(theory)
1080
1081    if isinstance(resolution, Slit1D):
1082        width, height = resolution.dqx, resolution.dqy
1083        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
1084    else:
1085        dq = resolution.q_width
1086        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars)
1087
1088    import matplotlib.pyplot as plt  # type: ignore
1089    plt.loglog(resolution.q_calc, theory, label='unsmeared')
1090    plt.loglog(resolution.q, Iq, label='smeared', hold=True)
1091    plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True)
1092    plt.legend()
1093    plt.title(title)
1094    plt.xlabel("Q (1/Ang)")
1095    plt.ylabel("I(Q) (1/cm)")
1096
1097def demo_pinhole_1d():
1098    """
1099    Show example of pinhole smearing.
1100    """
1101    q = np.logspace(-4, np.log10(0.2), 400)
1102    q_width = 0.1*q
1103    resolution = Pinhole1D(q, q_width)
1104    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
1105
1106def demo_slit_1d():
1107    """
1108    Show example of slit smearing.
1109    """
1110    q = np.logspace(-4, np.log10(0.2), 100)
1111    w = h = 0.
1112    #w = 0.000000277790
1113    w = 0.0277790
1114    #h = 0.00277790
1115    #h = 0.0277790
1116    resolution = Slit1D(q, w, h)
1117    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
1118
1119def demo():
1120    """
1121    Run the resolution demos.
1122    """
1123    import matplotlib.pyplot as plt  # type: ignore
1124    plt.subplot(121)
1125    demo_pinhole_1d()
1126    #plt.yscale('linear')
1127    plt.subplot(122)
1128    demo_slit_1d()
1129    #plt.yscale('linear')
1130    plt.show()
1131
1132
1133if __name__ == "__main__":
1134    #demo()
1135    main()
Note: See TracBrowser for help on using the repository browser.