source: sasmodels/sasmodels/resolution.py @ 346bc88

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 346bc88 was 346bc88, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

Add 2D resolution smearing. Split BumpsModel? into Experiment and Model and fix up examples.

  • Property mode set to 100644
File size: 31.4 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
7from scipy.special import erf
8from numpy import sqrt, log, log10
9import numpy as np
10
11MINIMUM_RESOLUTION = 1e-8
12
13class Resolution(object):
14    """
15    Abstract base class defining a 1D resolution function.
16
17    *q* is the set of q values at which the data is measured.
18
19    *q_calc* is the set of q values at which the theory needs to be evaluated.
20    This may extend and interpolate the q values.
21
22    *apply* is the method to call with I(q_calc) to compute the resolution
23    smeared theory I(q).
24    """
25    q = None
26    q_calc = None
27    def apply(self, theory):
28        """
29        Smear *theory* by the resolution function, returning *Iq*.
30        """
31        raise NotImplementedError("Subclass does not define the apply function")
32
33
34class Perfect1D(Resolution):
35    """
36    Resolution function to use when there is no actual resolution smearing
37    to be applied.  It has the same interface as the other resolution
38    functions, but returns the identity function.
39    """
40    def __init__(self, q):
41        self.q_calc = self.q = q
42
43    def apply(self, theory):
44        return theory
45
46
47class Pinhole1D(Resolution):
48    r"""
49    Pinhole aperture with q-dependent gaussian resolution.
50
51    *q* points at which the data is measured.
52
53    *q_width* gaussian 1-sigma resolution at each data point.
54
55    *q_calc* is the list of points to calculate, or None if this should
56    be estimated from the *q* and *q_width*.
57    """
58    def __init__(self, q, q_width, q_calc=None, nsigma=3):
59        #*min_step* is the minimum point spacing to use when computing the
60        #underlying model.  It should be on the order of
61        #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes
62        #are computed with sufficient density to avoid aliasing effects.
63
64        # Protect against calls with q_width=0.  The extend_q function will
65        # not extend the q if q_width is 0, but q_width must be non-zero when
66        # constructing the weight matrix to avoid division by zero errors.
67        # In practice this should never be needed, since resolution should
68        # default to Perfect1D if the pinhole geometry is not defined.
69        self.q, self.q_width = q, q_width
70        self.q_calc = pinhole_extend_q(q, q_width, nsigma=nsigma) \
71            if q_calc is None else np.sort(q_calc)
72        self.weight_matrix = pinhole_resolution(self.q_calc,
73                self.q, np.maximum(q_width, MINIMUM_RESOLUTION))
74
75    def apply(self, theory):
76        print "q calc", self.q_calc
77        print "weights", self.weight_matrix.shape
78        return apply_resolution_matrix(self.weight_matrix, theory)
79
80
81class Slit1D(Resolution):
82    """
83    Slit aperture with a complicated resolution function.
84
85    *q* points at which the data is measured.
86
87    *qx_width* slit width
88
89    *qy_height* slit height
90
91    *q_calc* is the list of points to calculate, or None if this should
92    be estimated from the *q* and *q_width*.
93
94    The *weight_matrix* is computed by :func:`slit1d_resolution`
95    """
96    def __init__(self, q, width, height, q_calc=None):
97        # TODO: maybe issue warnings rather than raising errors
98        if not np.isscalar(width):
99            if np.any(np.diff(width) > 0.0):
100                raise ValueError("Slit resolution requires fixed width slits")
101            width = width[0]
102        if not np.isscalar(height):
103            if np.any(np.diff(height) > 0.0):
104                raise ValueError("Slit resolution requires fixed height slits")
105            height = height[0]
106
107        # Remember what width/height was used even though we won't need them
108        # after the weight matrix is constructed
109        self.width, self.height = width, height
110
111        self.q = q.flatten()
112        self.q_calc = slit_extend_q(q, width, height) \
113            if q_calc is None else np.sort(q_calc)
114        self.weight_matrix = \
115            slit_resolution(self.q_calc, self.q, width, height)
116
117    def apply(self, theory):
118        return apply_resolution_matrix(self.weight_matrix, theory)
119
120
121def apply_resolution_matrix(weight_matrix, theory):
122    """
123    Apply the resolution weight matrix to the computed theory function.
124    """
125    #print "apply shapes", theory.shape, weight_matrix.shape
126    Iq = np.dot(theory[None,:], weight_matrix)
127    #print "result shape",Iq.shape
128    return Iq.flatten()
129
130
131def pinhole_resolution(q_calc, q, q_width):
132    """
133    Compute the convolution matrix *W* for pinhole resolution 1-D data.
134
135    Each row *W[i]* determines the normalized weight that the corresponding
136    points *q_calc* contribute to the resolution smeared point *q[i]*.  Given
137    *W*, the resolution smearing can be computed using *dot(W,q)*.
138
139    *q_calc* must be increasing.  *q_width* must be greater than zero.
140    """
141    # The current algorithm is a midpoint rectangle rule.  In the test case,
142    # neither trapezoid nor Simpson's rule improved the accuracy.
143    edges = bin_edges(q_calc)
144    edges[edges<0.0] = 0.0 # clip edges below zero
145    G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] )
146    weights = G[1:] - G[:-1]
147    weights /= np.sum(weights, axis=0)[None,:]
148    return weights
149
150
151def slit_resolution(q_calc, q, width, height):
152    r"""
153    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given
154    $q_v$ = *width* and $q_h$ = *height*.
155
156    *width* and *height* are scalars since current instruments use the
157    same slit settings for all measurement points.
158
159    If slit height is large relative to width, use:
160
161    .. math::
162
163        I_s(q_o) = \frac{1}{\Delta q_v}
164            \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2} du
165
166    If slit width is large relative to height, use:
167
168    .. math::
169
170        I_s(q_o) = \frac{1}{2 \Delta q_v}
171            \int_{-\Delta q_v}^{\Delta q_v} I(u) du
172    """
173    if width == 0.0 and height == 0.0:
174        #print "condition zero"
175        return 1
176
177    q_edges = bin_edges(q_calc) # Note: requires q > 0
178    q_edges[q_edges<0.0] = 0.0 # clip edges below zero
179
180    #np.set_printoptions(linewidth=10000)
181    if width <= 100.0 * height or height == 0:
182        # The current algorithm is a midpoint rectangle rule.  In the test case,
183        # neither trapezoid nor Simpson's rule improved the accuracy.
184        #print "condition h", q_edges.shape, q.shape, q_calc.shape
185        weights = np.zeros((len(q), len(q_calc)), 'd')
186        for i, qi in enumerate(q):
187            weights[i, :] = np.diff(q_to_u(q_edges, qi))
188        weights /= width
189        weights = weights.T
190    else:
191        #print "condition w"
192        # Make q_calc into a row vector, and q into a column vector
193        q, q_calc = q[None,:], q_calc[:,None]
194        in_x = (q_calc >= q-width) * (q_calc <= q+width)
195        weights = np.diff(q_edges)[:,None] * in_x
196
197    return weights
198
199
200def pinhole_extend_q(q, q_width, nsigma=3):
201    """
202    Given *q* and *q_width*, find a set of sampling points *q_calc* so
203    that each point I(q) has sufficient support from the underlying
204    function.
205    """
206    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width)
207    return linear_extrapolation(q, q_min, q_max)
208
209
210def slit_extend_q(q, width, height):
211    """
212    Given *q*, *width* and *height*, find a set of sampling points *q_calc* so
213    that each point I(q) has sufficient support from the underlying
214    function.
215    """
216    height # keep lint happy
217    q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2))
218    return geometric_extrapolation(q, q_min, q_max)
219
220
221def bin_edges(x):
222    """
223    Determine bin edges from bin centers, assuming that edges are centered
224    between the bins.
225
226    Note: this uses the arithmetic mean, which may not be appropriate for
227    log-scaled data.
228    """
229    if len(x) < 2 or (np.diff(x)<0).any():
230        raise ValueError("Expected bins to be an increasing set")
231    edges = np.hstack([
232        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
233        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
234        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
235        ])
236    return edges
237
238
239def q_to_u(q, q0):
240    """
241    Convert *q* values to *u* values for the integral computed at *q0*.
242    """
243    # array([value])**2 - value**2 is not always zero
244    qpsq = q**2 - q0**2
245    qpsq[qpsq<0] = 0
246    return sqrt(qpsq)
247
248
249def interpolate(q, max_step):
250    """
251    Returns *q_calc* with points spaced at most max_step apart.
252    """
253    step = np.diff(q)
254    index = step>max_step
255    if np.any(index):
256        inserts = []
257        for q_i,step_i in zip(q[:-1][index],step[index]):
258            n = np.ceil(step_i/max_step)
259            inserts.extend(q_i + np.arange(1,n)*(step_i/n))
260        # Extend a couple of fringes beyond the end of the data
261        inserts.extend(q[-1] + np.arange(1,8)*max_step)
262        q_calc = np.sort(np.hstack((q,inserts)))
263    else:
264        q_calc = q
265    return q_calc
266
267
268def linear_extrapolation(q, q_min, q_max):
269    """
270    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as
271    a guide.  Extrapolation below uses about the same size as the first
272    interval.  Extrapolation above uses about the same size as the final
273    interval.
274
275    if *q_min* is zero or less then *q[0]/10* is used instead.
276    """
277    q = np.sort(q)
278    if q_min < q[0]:
279        if q_min <= 0: q_min = q[0]/10
280        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15
281        q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
282    else:
283        q_low = []
284    if q_max > q[-1]:
285        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1]>q[-2] else 15
286        q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
287    else:
288        q_high = []
289    return np.concatenate([q_low, q, q_high])
290
291
292def geometric_extrapolation(q, q_min, q_max, points_per_decade=None):
293    r"""
294    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the
295    average geometric step size in *q* as the step size.
296
297    if *q_min* is zero or less then *q[0]/10* is used instead.
298
299    *points_per_decade* sets the ratio between consecutive steps such
300    that there will be $n$ points used for every factor of 10 increase
301    in *q*.
302
303    If *points_per_decade* is not given, it will be estimated as follows.
304    Starting at $q_1$ and stepping geometrically by $\Delta q$ to $q_n$
305    in $n$ points gives a geometric average of:
306
307    .. math::
308
309         \log \Delta q = (\log q_n - log q_1) / (n - 1)
310
311    From this we can compute the number of steps required to extend $q$
312    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
313
314    .. math::
315
316         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
317
318    Substituting:
319
320        n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
321            / (\log q_n - log q_1)
322    """
323    q = np.sort(q)
324    if points_per_decade is None:
325        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0]))
326    else:
327        log_delta_q = log(10.) / points_per_decade
328    if q_min < q[0]:
329        if q_min < 0: q_min = q[0]/10
330        n_low = log_delta_q * (log(q[0])-log(q_min))
331        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
332    else:
333        q_low = []
334    if q_max > q[-1]:
335        n_high = log_delta_q * (log(q_max)-log(q[-1]))
336        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
337    else:
338        q_high = []
339    return np.concatenate([q_low, q, q_high])
340
341
342############################################################################
343# unit tests
344############################################################################
345import unittest
346
347
348def eval_form(q, form, pars):
349    from sasmodels import core
350    kernel = core.make_kernel(form, [q])
351    theory = core.call_kernel(kernel, pars)
352    kernel.release()
353    return theory
354
355
356def gaussian(q, q0, dq):
357    from numpy import exp, pi
358    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)
359
360
361def romberg_slit_1d(q, delta_qv, form, pars):
362    """
363    Romberg integration for slit resolution.
364
365    This is an adaptive integration technique.  It is called with settings
366    that make it slow to evaluate but give it good accuracy.
367    """
368    from scipy.integrate import romberg
369
370    if any(k not in form.info['defaults'] for k in pars.keys()):
371        keys = set(form.info['defaults'].keys())
372        extra = set(pars.keys()) - keys
373        raise ValueError("bad parameters: [%s] not in [%s]"%
374                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
375
376    _fn = lambda u, q0: eval_form(sqrt(q0**2 + u**2), form, pars)
377    r = [romberg(_fn, 0, delta_qv[0], args=(qi,),
378                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
379         for qi in q]
380    # r should be [float, ...], but it is [array([float]), array([float]),...]
381    return np.asarray(r).flatten()/delta_qv[0]
382
383
384def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5):
385    """
386    Romberg integration for pinhole resolution.
387
388    This is an adaptive integration technique.  It is called with settings
389    that make it slow to evaluate but give it good accuracy.
390    """
391    from scipy.integrate import romberg
392
393    if any(k not in form.info['defaults'] for k in pars.keys()):
394        keys = set(form.info['defaults'].keys())
395        extra = set(pars.keys()) - keys
396        raise ValueError("bad parameters: [%s] not in [%s]"%
397                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
398
399    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
400    r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi),
401                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
402         for qi,dqi in zip(q,q_width)]
403    return np.asarray(r).flatten()
404
405
406class ResolutionTest(unittest.TestCase):
407
408    def setUp(self):
409        self.x = 0.001*np.arange(1, 11)
410        self.y = self.Iq(self.x)
411
412    def Iq(self, q):
413        "Linear function for resolution unit test"
414        return 12.0 - 1000.0*q
415
416    def test_perfect(self):
417        """
418        Perfect resolution and no smearing.
419        """
420        resolution = Perfect1D(self.x)
421        theory = self.Iq(resolution.q_calc)
422        output = resolution.apply(theory)
423        np.testing.assert_equal(output, self.y)
424
425    def test_slit_zero(self):
426        """
427        Slit smearing with perfect resolution.
428        """
429        resolution = Slit1D(self.x, width=0, height=0, q_calc=self.x)
430        theory = self.Iq(resolution.q_calc)
431        output = resolution.apply(theory)
432        np.testing.assert_equal(output, self.y)
433
434    @unittest.skip("not yet supported")
435    def test_slit_high(self):
436        """
437        Slit smearing with height 0.005
438        """
439        resolution = Slit1D(self.x, width=0, height=0.005, q_calc=self.x)
440        theory = self.Iq(resolution.q_calc)
441        output = resolution.apply(theory)
442        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
443                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ]
444        np.testing.assert_allclose(output, answer, atol=1e-4)
445
446    @unittest.skip("not yet supported")
447    def test_slit_both_high(self):
448        """
449        Slit smearing with width < 100*height.
450        """
451        q = np.logspace(-4, -1, 10)
452        resolution = Slit1D(q, width=0.2, height=np.inf)
453        theory = 1000*self.Iq(resolution.q_calc**4)
454        output = resolution.apply(theory)
455        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
456                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ]
457        np.testing.assert_allclose(output, answer, atol=1e-4)
458
459    @unittest.skip("not yet supported")
460    def test_slit_wide(self):
461        """
462        Slit smearing with width 0.0002
463        """
464        resolution = Slit1D(self.x, width=0.0002, height=0, q_calc=self.x)
465        theory = self.Iq(resolution.q_calc)
466        output = resolution.apply(theory)
467        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
468        np.testing.assert_allclose(output, answer, atol=1e-4)
469
470    @unittest.skip("not yet supported")
471    def test_slit_both_wide(self):
472        """
473        Slit smearing with width > 100*height.
474        """
475        resolution = Slit1D(self.x, width=0.0002, height=0.000001,
476                            q_calc=self.x)
477        theory = self.Iq(resolution.q_calc)
478        output = resolution.apply(theory)
479        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
480        np.testing.assert_allclose(output, answer, atol=1e-4)
481
482    def test_pinhole_zero(self):
483        """
484        Pinhole smearing with perfect resolution
485        """
486        resolution = Pinhole1D(self.x, 0.0*self.x)
487        theory = self.Iq(resolution.q_calc)
488        output = resolution.apply(theory)
489        np.testing.assert_equal(output, self.y)
490
491    def test_pinhole(self):
492        """
493        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001]
494        """
495        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x),
496                               q_calc=self.x)
497        theory = 12.0-1000.0*resolution.q_calc
498        output = resolution.apply(theory)
499        answer = [ 10.44785079, 9.84991299, 8.98101708,
500                  7.99906585, 6.99998311, 6.00001689,
501                  5.00093415, 4.01898292, 3.15008701, 2.55214921]
502        np.testing.assert_allclose(output, answer, atol=1e-8)
503
504
505class IgorComparisonTest(unittest.TestCase):
506
507    def setUp(self):
508        self.pars = TEST_PARS_PINHOLE_SPHERE
509        from sasmodels import core
510        from sasmodels.models import sphere
511        self.model = core.load_model(sphere, dtype='double')
512
513    def Iq_sphere(self, pars, resolution):
514        from sasmodels import core
515        kernel = core.make_kernel(self.model, [resolution.q_calc])
516        theory = core.call_kernel(kernel, pars)
517        result = resolution.apply(theory)
518        kernel.release()
519        return result
520
521    def compare(self, q, output, answer, tolerance):
522        err = (output - answer)/answer
523        idx = abs(err) >= tolerance
524        problem = zip(q[idx], output[idx], answer[idx], err[idx])
525        print "\n".join(str(v) for v in problem)
526        np.testing.assert_allclose(output, answer, rtol=tolerance)
527
528    def test_perfect(self):
529        """
530        Compare sphere model with NIST Igor SANS, no resolution smearing.
531        """
532        pars = TEST_PARS_SLIT_SPHERE
533        data_string = TEST_DATA_SLIT_SPHERE
534
535        data = np.loadtxt(data_string.split('\n')).T
536        q, width, answer, _ = data
537        resolution = Perfect1D(q)
538        output = self.Iq_sphere(pars, resolution)
539        self.compare(q, output, answer, 1e-6)
540
541    def test_pinhole(self):
542        """
543        Compare pinhole resolution smearing with NIST Igor SANS
544        """
545        pars = TEST_PARS_PINHOLE_SPHERE
546        data_string = TEST_DATA_PINHOLE_SPHERE
547
548        data = np.loadtxt(data_string.split('\n')).T
549        q, q_width, answer = data
550        resolution = Pinhole1D(q, q_width)
551        output = self.Iq_sphere(pars, resolution)
552        # TODO: relative error should be lower
553        self.compare(q, output, answer, 3e-4)
554
555    def test_pinhole_romberg(self):
556        """
557        Compare pinhole resolution smearing with romberg integration result.
558        """
559        pars = TEST_PARS_PINHOLE_SPHERE
560        data_string = TEST_DATA_PINHOLE_SPHERE
561        pars['radius'] *= 5
562        radius = pars['radius']
563
564        data = np.loadtxt(data_string.split('\n')).T
565        q, q_width, answer = data
566        answer = romberg_pinhole_1d(q, q_width, self.model, pars)
567        ## Getting 0.1% requires 5 sigma and 200 points per fringe
568        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5),
569        #                     2*np.pi/radius/200)
570        #tol = 0.001
571        ## The default 3 sigma and no extra points gets 1%
572        q_calc, tol = None, 0.01
573        resolution = Pinhole1D(q, q_width, q_calc=q_calc)
574        output = self.Iq_sphere(pars, resolution)
575        if 0: # debug plot
576            import matplotlib.pyplot as plt
577            resolution = Perfect1D(q)
578            source = self.Iq_sphere(pars, resolution)
579            plt.loglog(q, source, '.')
580            plt.loglog(q, answer, '-', hold=True)
581            plt.loglog(q, output, '-', hold=True)
582            plt.show()
583        self.compare(q, output, answer, tol)
584
585    def test_slit(self):
586        """
587        Compare slit resolution smearing with NIST Igor SANS
588        """
589        pars = TEST_PARS_SLIT_SPHERE
590        data_string = TEST_DATA_SLIT_SPHERE
591
592        data = np.loadtxt(data_string.split('\n')).T
593        q, delta_qv, _, answer = data
594        resolution = Slit1D(q, width=delta_qv, height=0)
595        output = self.Iq_sphere(pars, resolution)
596        # TODO: eliminate Igor test since it is too inaccurate to be useful.
597        # This means we can eliminate the test data as well, and instead
598        # use a generated q vector.
599        self.compare(q, output, answer, 0.5)
600
601    def test_slit_romberg(self):
602        """
603        Compare slit resolution smearing with romberg integration result.
604        """
605        pars = TEST_PARS_SLIT_SPHERE
606        data_string = TEST_DATA_SLIT_SPHERE
607        radius = pars['radius']
608
609        data = np.loadtxt(data_string.split('\n')).T
610        q, delta_qv, _, answer = data
611        answer = romberg_slit_1d(q, delta_qv, self.model, pars)
612        q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20),
613                               delta_qv[0], delta_qv[0])
614        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc)
615        output = self.Iq_sphere(pars, resolution)
616        # TODO: relative error should be lower
617        self.compare(q, output, answer, 0.025)
618
619    def test_ellipsoid(self):
620        """
621        Compare romberg integration for ellipsoid model.
622        """
623        from .core import load_model
624        pars = {
625            'scale':0.05,
626            'rpolar':500, 'requatorial':15000,
627            'sld':6, 'solvent_sld': 1,
628            }
629        form = load_model('ellipsoid', dtype='double')
630        q = np.logspace(log10(4e-5),log10(2.5e-2), 68)
631        delta_qv = [0.117]
632        resolution = Slit1D(q, width=delta_qv, height=0)
633        answer = romberg_slit_1d(q, delta_qv, form, pars)
634        output = resolution.apply(eval_form(resolution.q_calc, form, pars))
635        # TODO: 10% is too much error; use better algorithm
636        self.compare(q, output, answer, 0.1)
637
638    #TODO: can sas q spacing be too sparse for the resolution calculation?
639    @unittest.skip("suppress sparse data test; not supported by current code")
640    def test_pinhole_sparse(self):
641        """
642        Compare pinhole resolution smearing with NIST Igor SANS on sparse data
643        """
644        pars = TEST_PARS_PINHOLE_SPHERE
645        data_string = TEST_DATA_PINHOLE_SPHERE
646
647        data = np.loadtxt(data_string.split('\n')).T
648        q, q_width, answer = data[:, ::20] # Take every nth point
649        resolution = Pinhole1D(q, q_width)
650        output = self.Iq_sphere(pars, resolution)
651        self.compare(q, output, answer, 1e-6)
652
653
654# pinhole sphere parameters
655TEST_PARS_PINHOLE_SPHERE = {
656    'scale': 1.0, 'background': 0.01,
657    'radius': 60.0, 'sld': 1, 'solvent_sld': 6.3,
658    }
659# Q, dQ, I(Q) calculated by NIST Igor SANS package
660TEST_DATA_PINHOLE_SPHERE = """\
6610.001278 0.0002847 2538.41176383
6620.001562 0.0002905 2536.91820405
6630.001846 0.0002956 2535.13182479
6640.002130 0.0003017 2533.06217813
6650.002414 0.0003087 2530.70378586
6660.002698 0.0003165 2528.05024192
6670.002982 0.0003249 2525.10408349
6680.003266 0.0003340 2521.86667499
6690.003550 0.0003437 2518.33907750
6700.003834 0.0003539 2514.52246995
6710.004118 0.0003646 2510.41798319
6720.004402 0.0003757 2506.02690988
6730.004686 0.0003872 2501.35067884
6740.004970 0.0003990 2496.38678318
6750.005253 0.0004112 2491.16237596
6760.005537 0.0004237 2485.63911673
6770.005821 0.0004365 2479.83657083
6780.006105 0.0004495 2473.75676948
6790.006389 0.0004628 2467.40145990
6800.006673 0.0004762 2460.77293372
6810.006957 0.0004899 2453.86724627
6820.007241 0.0005037 2446.69623838
6830.007525 0.0005177 2439.25775219
6840.007809 0.0005318 2431.55421398
6850.008093 0.0005461 2423.58785521
6860.008377 0.0005605 2415.36158137
6870.008661 0.0005750 2406.87009473
6880.008945 0.0005896 2398.12841186
6890.009229 0.0006044 2389.13360806
6900.009513 0.0006192 2379.88958042
6910.009797 0.0006341 2370.39776774
6920.010080 0.0006491 2360.69528793
6930.010360 0.0006641 2350.85169027
6940.010650 0.0006793 2340.42023633
6950.010930 0.0006945 2330.11206013
6960.011220 0.0007097 2319.20109972
6970.011500 0.0007251 2308.43503981
6980.011780 0.0007404 2297.44820179
6990.012070 0.0007558 2285.83853677
7000.012350 0.0007713 2274.41290746
7010.012640 0.0007868 2262.36219581
7020.012920 0.0008024 2250.51169731
7030.013200 0.0008180 2238.45596231
7040.013490 0.0008336 2225.76495666
7050.013770 0.0008493 2213.29618391
7060.014060 0.0008650 2200.19110751
7070.014340 0.0008807 2187.34050325
7080.014620 0.0008965 2174.30529864
7090.014910 0.0009123 2160.61632548
7100.015190 0.0009281 2147.21038112
7110.015470 0.0009440 2133.62023580
7120.015760 0.0009598 2119.37907426
7130.016040 0.0009757 2105.45234903
7140.016330 0.0009916 2090.86319102
7150.016610 0.0010080 2076.60576032
7160.016890 0.0010240 2062.19214565
7170.017180 0.0010390 2047.10550219
7180.017460 0.0010550 2032.38715621
7190.017740 0.0010710 2017.52560123
7200.018030 0.0010880 2001.99124318
7210.018310 0.0011040 1986.84662060
7220.018600 0.0011200 1971.03389745
7230.018880 0.0011360 1955.61395119
7240.019160 0.0011520 1940.08291563
7250.019450 0.0011680 1923.87672225
7260.019730 0.0011840 1908.10656374
7270.020020 0.0012000 1891.66297192
7280.020300 0.0012160 1875.66789021
7290.020580 0.0012320 1859.56357196
7300.020870 0.0012490 1842.79468290
7310.021150 0.0012650 1826.50064489
7320.021430 0.0012810 1810.11533702
7330.021720 0.0012970 1793.06840882
7340.022000 0.0013130 1776.51153580
7350.022280 0.0013290 1759.87201249
7360.022570 0.0013460 1742.57354412
7370.022850 0.0013620 1725.79397319
7380.023140 0.0013780 1708.35831550
7390.023420 0.0013940 1691.45256069
7400.023700 0.0014110 1674.48561783
7410.023990 0.0014270 1656.86525366
7420.024270 0.0014430 1639.79847285
7430.024550 0.0014590 1622.68887088
7440.024840 0.0014760 1604.96421100
7450.025120 0.0014920 1587.85768129
7460.025410 0.0015080 1569.99297335
7470.025690 0.0015240 1552.84580279
7480.025970 0.0015410 1535.54074115
7490.026260 0.0015570 1517.75249337
7500.026540 0.0015730 1500.40115023
7510.026820 0.0015900 1483.03632237
7520.027110 0.0016060 1465.05942429
7530.027390 0.0016220 1447.67682181
7540.027670 0.0016390 1430.46495191
7550.027960 0.0016550 1412.49232282
7560.028240 0.0016710 1395.13182318
7570.028520 0.0016880 1377.93439837
7580.028810 0.0017040 1359.99528971
7590.029090 0.0017200 1342.67274512
7600.029370 0.0017370 1325.55375609
761"""
762
763# Slit sphere parameters
764TEST_PARS_SLIT_SPHERE = {
765    'scale': 0.01, 'background': 0.01,
766    'radius': 60000, 'sld': 1, 'solvent_sld': 4,
767    }
768# Q dQ I(Q) I_smeared(Q)
769TEST_DATA_SLIT_SPHERE = """\
7702.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
7712.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
7722.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
7733.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
7743.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
7753.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
7764.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
7775.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
7785.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
7796.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
7806.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
7817.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
7827.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
7838.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
7848.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
7859.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
7861.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
7871.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
7881.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
7891.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
7901.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
7911.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
7921.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
7931.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
7941.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
7951.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
7962.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
7972.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
7982.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
7992.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
8002.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
8012.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
8022.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
8032.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
8043.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
8053.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
8063.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
8073.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
8084.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
8094.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
8104.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
8114.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
8125.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
8135.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
8146.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
8156.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
8167.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
8177.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
8188.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
8198.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
8209.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
8219.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
8221.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
8231.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
8241.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
8251.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
8261.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
8271.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
8281.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
8291.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
8301.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
8311.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
8322.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
8332.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
8342.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
8352.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
8362.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
8372.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
838"""
839
840def main():
841    """
842    Run tests given is sys.argv.
843
844    Returns 0 if success or 1 if any tests fail.
845    """
846    import sys
847    import xmlrunner
848
849    suite = unittest.TestSuite()
850    suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__]))
851
852    runner = xmlrunner.XMLTestRunner(output='logs')
853    result = runner.run(suite)
854    return 1 if result.failures or result.errors else 0
855
856
857############################################################################
858# usage demo
859############################################################################
860
861def _eval_demo_1d(resolution, title):
862    from sasmodels import core
863    from sasmodels.models import cylinder
864    ## or alternatively:
865    # cylinder = core.load_model_definition('cylinder')
866    model = core.load_model(cylinder)
867
868    kernel = core.make_kernel(model, [resolution.q_calc])
869    theory = core.call_kernel(kernel, {'length':210, 'radius':500})
870    Iq = resolution.apply(theory)
871
872    import matplotlib.pyplot as plt
873    plt.loglog(resolution.q_calc, theory, label='unsmeared')
874    plt.loglog(resolution.q, Iq, label='smeared', hold=True)
875    plt.legend()
876    plt.title(title)
877    plt.xlabel("Q (1/Ang)")
878    plt.ylabel("I(Q) (1/cm)")
879
880def demo_pinhole_1d():
881    q = np.logspace(-3, -1, 400)
882    q_width = 0.1*q
883    resolution = Pinhole1D(q, q_width)
884    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution")
885
886def demo_slit_1d():
887    q = np.logspace(-3, -1, 400)
888    qx_width = 0.005
889    qy_width = 0.0
890    resolution = Slit1D(q, qx_width, qy_width)
891    _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution")
892
893def demo():
894    import matplotlib.pyplot as plt
895    plt.subplot(121)
896    demo_pinhole_1d()
897    plt.subplot(122)
898    demo_slit_1d()
899    plt.show()
900
901
902if __name__ == "__main__":
903    #demo()
904    main()
Note: See TracBrowser for help on using the repository browser.