source: sasmodels/sasmodels/resolution.py @ b397165

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b397165 was b397165, checked in by GitHub <noreply@…>, 4 years ago

Revert "Jurtest2"

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