source: sasmodels/sasmodels/resolution.py @ 3199b17

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3199b17 was 9e7837a, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

use asymmetric integration window for resolution, (-2.5,+3.0) sigma

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