source: sasmodels/sasmodels/resolution.py @ 355ee8b

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 355ee8b was 355ee8b, checked in by Paul Kienzle <pkienzle@…>, 3 years ago

resolution: trim extrapolated q values in [-0.02*min(q), 0.02*min(q)]

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