source: sasmodels/sasmodels/resolution.py @ 7b7fcf0

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

allow negative q in resolution calculation

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