Ticket #312: slit_test.py

File slit_test.py, 18.7 KB (added by pkienzle, 10 years ago)
Line 
1r"""
2Explore different integration schemes for slit smearing.
3
4Introduction
5============
6
7Slit smearing evaluates the following integral:
8
9.. math::
10
11    I_s(q_o) = \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2}) du
12
13where $\Delta q_v$ is the vertical slit contribution. See
14`http://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_61.pdf`_
15
16Procedure
17=========
18
19This program is not intended as an end user application, but rather as a
20platform for exploring implementations of the resolution integral.
21Different experiments can be performed by modifying the :func:`demo` function.
22The integration functions have the name smear_METHOD, for a variety of
23different methods.  Each one takes a form factor with parameters, a set of
24q values, a resolution, and maybe a set of q_calc values at which to evaluate
25the theory function.   The *ms* and *mt* variables in :func:`demo`, determine
26which methods are being compared.  :func:`interpolate` allows you to fill
27in the integral between q points since the demo function is rapidly
28oscillating.  *TEST_DATA_SLIT_SMEAR* contains the values calculated in
29Igor for a sphere model with a radius of 6 microns and a contrast of 3e-6
30inverse Angstroms.
31
32:func:`smear_romberg` is a direct implementation of the integral using
33adaptive romberg integration from scipy.  This method would add considerable
34complexity to sasmodels since we could no longer rely on evaluating the
35models at a fixed set of q points for all parameter sets in a fit.  A
36"good enough" method with fixed q, although less accurate and requiring
37more q points on average to evaluate, will be a win in complexity and perhaps
38performance.  The romberg integration parameters are set to produce a
39pretty good value for the integral.
40
41The remaining methods used a fixed set of *q* points, *q_calc*, to evaluate the
42integral.   They are :func:`smear_simpsons`,  :func:`smear_trapezoid`,
43and :func:`smear_rectangular`.  Except for :func:`smear_simpsons_pow`,
44*q_calc* is extrapolated to the full width of the resolution function.
45The extrapolation uses exponential spacing between the points, much like *q*
46spacing within the measured data range, though :func:`smear_simpsons_lin`
47uses linear spacing.
48
49:func:`smear_rectangular` is implemented as a non-square weight matrix which
50can be multiplied by a set of computed values *I(q_calc)* to produce
51*I_smeared(q)*.  The same style could be used for trapezoid and Simpson's
52rule with the appropriate bookkeeping.  This was not done because at least
53for the demo problem, rectangular performed better than the other
54non-adaptive methods, and so wins on both simplicity and accuracy.
55
56:func:`smear_simpsons_pow` uses simpson's rule to estimate the integral within
57the range of the data, and a power law approximation for the tail above the
58maximum q value.  This is akin to the Igor implementation of slit smearing,
59though it uses the theory function itself to estimate the power law rather
60than data.
61
62The power law approximation has two problems:
63
641. the power is estimated from the data, which means that this will only
65   work if the data has a flat region at the tail, and
662. any features such as lamellar peaks that occur beyond the end of the
67   data will not be correctly modeled.
68
69The demo shows that this method performs poorly for high q when oscillations
70in the theory have not completely decayed by the end of the q range, which
71is not really a surprise.  It is not a problem in practice for usans data on
72igor since polydispersity will likely dampen any ringing that a monodisperse
73model might exhibit.  The demo does not simulate a model containing a peak at
74q beyond the end of the data.
75
76
77Conclusions
78===========
79
80At least for this theory function, the simple rectangular integral is as
81good as the more complicated trapezoid and Simpson's rule integration.
82Due to the oscillations in the theory function, the relative error is
83huge (50%) at high q.  Oversampling by a factor of 10 (requiring 1020
84evaluations rather than 122), reduces the error to 2.5%.  Another factor
85of 10 (9130 evaluations) reduces error to 0.15%. Getting error down to
865e-5 requires about 90,000 evaluations.
87
88In comparison, romberg integration requires only 1360 evaluations.
89
90"""
91import numpy as np
92from numpy import sqrt, cos, sin, pi, log10, log, exp
93
94def smear_romberg(q, form, pars, delta_qv, q_calc=None):
95    """
96    Romberg integration.
97
98    This is an adaptive integration technique.  It is called with settings
99    that make it slow to evaluate but give it good accuracy.
100    """
101    from scipy.integrate import romberg
102    evaluations = [0]
103    def _fn(u, q0):
104        evaluations[0] += 1 if np.isscalar(q0) else len(q0)
105        return form(sqrt(q0**2 + u**2), *pars)
106    r = [romberg(_fn, 0, delta_qv, args=(qi,),
107                 divmax=100, vec_func=True, tol=0, rtol=1e-14)
108         for qi in q]
109    print "romberg evaluations:", evaluations[0]
110    return r
111
112
113def smear_simpsons_pow(q, form, pars, delta_qv, q_calc=None):
114    """
115    Use Simpson's rule with power low approximation to compute the integral.
116
117    The *q* range between the end of the data and *delta_qv* is not sampled.
118    Instead it is approximated by a power law computed from the (I(delta_qv)) and
119    the final 25\% of the theory function.  Igor fits the data to the a power
120    law and assumes that it will be an adequate estimate for the fitted
121    theory function.
122    """
123    from scipy.integrate import simps
124    if q_calc is None: q_calc = q
125    x = np.hstack((q_calc, q_calc[-1]+delta_qv))
126    print "#q:",q.size, "#q_calc:",x.size
127    y = form(x, *pars)
128    A, m = fit_power_law(x, y, portion=0.25)
129    return [simps(y[i:-1], q_to_u(x[i:-1], q[i]))
130            + power_law_residual(delta_qv, q[i], q[-1], A, m)
131            for i in range(len(q))]
132
133
134def smear_simpsons_lin(q, form, pars, delta_qv, q_calc=None):
135    """
136    Use Simpson's rule with linear extrapolation to compute the integral.
137
138    The evaluation points *q* are extrapolated to *delta_qv+q[-1]* using
139    the final *q* step as the step size.
140
141    For a typical NCNR USANS dataset this requires over 10000 extra evaluations
142    of the function.
143    """
144    from scipy.integrate import simps
145    if q_calc is None: q_calc = q
146    x = linear_extrapolation(q_calc, q_calc[-1]+delta_qv)
147    print "#q:",q.size, "#q_calc:",x.size
148    y = form(x, *pars)
149    return [simps(y[j:], q_to_u(x[j:], qi))
150            for qi in q for j in [np.searchsorted(x,qi)]]
151
152
153def smear_simpsons(q, form, pars, delta_qv, q_calc=None):
154    """
155    Use Simpson's rule with geometric extrapolation to compute the integral.
156
157    The evaluation points *q* are extrapolated to *delta_qv+q[-1]* using
158    the geometric average of *q* as the step size.
159
160    For a typical NCNR USANS dataset the function is evaluated at twice as
161    many points as there are points in q.
162    """
163    from scipy.integrate import simps
164    if q_calc is None: q_calc = q
165    x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv)
166    print "#q:",q.size, "#q_calc:",x.size
167    y = form(x, *pars)
168    return [simps(y[j:], q_to_u(x[j:], qi))
169            for qi in q for j in [np.searchsorted(x,qi)]]
170
171
172def smear_trapezoid(q, form, pars, delta_qv, q_calc=None):
173    """
174    Use the trapezoid rule with geometric extrapolation to compute the integral.
175    """
176    from scipy.integrate import trapz
177
178    if q_calc is None: q_calc = q
179    x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv)
180    print "#q:",q.size, "#q_calc:",x.size
181    y = form(x, *pars)
182    return [trapz(y[j:], q_to_u(x[j:], qi))
183            for qi in q for j in [np.searchsorted(x,qi)]]
184
185
186def smear_rectangular(q, form, pars, delta_qv, q_calc=None):
187    """
188    Use the trapezoid rule with geometric extrapolation to compute the integral.
189    """
190    from scipy.integrate import trapz
191    if q_calc is None: q_calc = q
192    x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv)
193    print "#q:",q.size, "#q_calc:",x.size
194    edges = bin_edges(x)
195    W = np.zeros((len(q), len(x)), 'd')
196    for i, qi in enumerate(q):
197        W[i,:] = np.diff(q_to_u(edges, qi))
198
199    y = form(x, *pars)
200    return np.dot(W,y)
201
202
203###########################################################################
204# Helper functions
205
206def sphere(q, radius, contrast):
207    """
208    Sphere form factor with no polydispersity
209    """
210    qr = q * radius
211    sn, cn = sin(qr), cos(qr)
212    bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
213    #bes[qr == 0] = 1
214    fq = bes * contrast
215    return 4*pi*radius**3/3 * 1.0e-4 * fq ** 2
216
217
218
219def q_to_u(q, q0):
220    """
221    Convert *q* values to *u* values for the integral computed at *q0*.
222    """
223    # array([value])**2 - value**2 is not always zero
224    qpsq = q**2 - q0**2
225    qpsq[qpsq<0] = 0
226    return sqrt(qpsq)
227
228
229def interpolate(q, max_step):
230    """
231    Returns *q_calc* with points spaced at most max_step apart.
232    """
233    step = np.diff(q)
234    index = step>max_step
235    if np.any(index):
236        inserts = []
237        for q_i,step_i in zip(q[:-1][index],step[index]):
238            n = np.ceil(step_i/max_step)
239            inserts.extend(q_i + np.arange(1,n)*(step_i/n))
240        # Extend a couple of fringes beyond the end of the data
241        inserts.extend(q[-1] + np.arange(1,16)*max_step)
242        q_calc = np.sort(np.hstack((q,inserts)))
243    else:
244        q_calc = q
245    return q_calc
246
247def linear_extrapolation(q, qmax):
248    """
249    Extrapolate *q* out to *qmax* using the final step in *q* as the step
250    size.
251    """
252    extrapolate = np.arange(q[-1], qmax, q[-1]-q[-2])
253    x = np.concatenate((q, extrapolate[1:]))
254
255
256def geometric_extrapolation(q, qmax):
257    r"""
258    Extrapolate *q* out to *qmax* using geometric steps, with the average
259    geometric step size in *q* as the step size.
260
261    Starting at $q_1$ and stepping geometrically by $\Delta q$
262    to $q_n$ in $n$ points gives a geometric average of:
263
264    .. math::
265
266         \log \Delta q = (\log q_n - log q_1) / (n - 1)
267
268    From this we can compute the number of steps required to extend $q$
269    from $q_n$ to $q_\text{max}$ by $\Delta q$ as:
270
271    .. math::
272
273         n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q
274
275    Substituting:
276
277        n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
278            / (\log q_n - log q_1)
279    """
280    n_extend = (len(q)-1)*(log(qmax) - log(q[-1]))/(log(q[-1]) - log(q[0]))
281    extrapolate = np.logspace(log10(q[-1]), log10(qmax), np.ceil(n_extend)+1)
282    x = np.concatenate((q, extrapolate[1:]))
283    return x
284
285
286def fit_power_law(x, y, portion=1.0):
287    """
288    Fit a power law $y=Ax^{-m}$ against *x*, *y*.  If *portion* is given, use
289    only the given portion from the tail of the data in the fit.
290    """
291    fitting_points = max(2, np.ceil(len(x)*portion))
292    p = np.polyfit(log(x[-fitting_points:]), log(y[-fitting_points:]), 1)
293    A,m = exp(p[1]), -p[0]
294    #print "A,m",A,m
295
296
297def power_law_residual(delta_qv, q0, qmax, A, m):
298    r"""
299    Residual term for the integration, assuming I(q) follows a power law for
300    high q.
301
302    This integrates $I_p(\sqrt(q_o^2 + u^2})$ by $u$ for $u$ corresponding to
303    $q_\text{max}$ to $u=\Delta q_v$, where $I_p$ is:
304
305    .. math::
306
307        I_p(q) = A q^{-m}
308    """
309    from scipy.integrate import romberg
310    a, b = sqrt(qmax**2 - q0**2), delta_qv
311    _residual_fn = lambda u, q0, A, m: A*(q0**2 + u**2)**(-m/2)
312    r = romberg(_residual_fn, a, b, args=(q0, A, m),
313                divmax=100, vec_func=True, tol=0, rtol=1e-8)
314    return r
315
316def bin_edges(x):
317    """
318    Determine bin edges from bin centers, assuming that edges are centered
319    between the bins.
320
321    Note: this uses the arithmetic mean, which may not be appropriate for
322    log-scaled data.
323    """
324    if len(x) < 2 or (np.diff(x)<0).any():
325        raise ValueError("Expected bins to be an increasing set")
326    edges = np.hstack([
327        x[0]  - 0.5*(x[1]  - x[0]),  # first point minus half first interval
328        0.5*(x[1:] + x[:-1]),        # mid points of all central intervals
329        x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval
330        ])
331    return edges
332
333def Iq_smeared(q, form, pars, scale=1, background=0., delta_qv=0.117,
334               method='trapezoid', max_step=None):
335    """
336    Compute slit-smeared version of *form(q;pars)*
337    """
338    q_calc = interpolate(q, max_step) if max_step is not None else q
339    integrator = globals()['smear_'+method]
340    res = integrator(q, form, pars, delta_qv, q_calc)
341    return scale/delta_qv*np.asarray(res) + background
342
343
344def Iq(q, form, pars, scale=1, background=0):
345    """
346    Compute unsmeared version of *form(q;pars)*
347    """
348    return scale*form(q, *pars) + background
349
350
351TEST_DATA_SLIT_SPHERE = """\
3522.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06
3532.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06
3542.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06
3553.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06
3563.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05
3573.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05
3584.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05
3595.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05
3605.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05
3616.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04
3626.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04
3637.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04
3647.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04
3658.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04
3668.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04
3679.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04
3681.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04
3691.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04
3701.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04
3711.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03
3721.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03
3731.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03
3741.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03
3751.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03
3761.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03
3771.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03
3782.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03
3792.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03
3802.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03
3812.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03
3822.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03
3832.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03
3842.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03
3852.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02
3863.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02
3873.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02
3883.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02
3893.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02
3904.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02
3914.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02
3924.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02
3934.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02
3945.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02
3955.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01
3966.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01
3976.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01
3987.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01
3997.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01
4008.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01
4018.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01
4029.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01
4039.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01
4041.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01
4051.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01
4061.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00
4071.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00
4081.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00
4091.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00
4101.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00
4111.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00
4121.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00
4131.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00
4142.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00
4152.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00
4162.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00
4172.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00
4182.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00
4192.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00
420"""
421
422def demo():
423    from matplotlib import pyplot as plt
424
425    ## ==== Set data and model parameters
426    delta_qv = 0.117
427    q = np.logspace(-5,-3,400)
428    #q = np.logspace(log10(2.3e-5),log10(2.7e-3),68*10)
429    #q = np.logspace(log10(2.3e-5),log10(2.7e-3),68)
430    if 1:  # use values from Igor
431        data = np.loadtxt(TEST_DATA_SLIT_SPHERE.split('\n')).T
432        q, dq, _, answer = data
433        delta_qv = dq[0]
434    else:
435        answer = None
436    scale, background = 0.01, 0. #0.01
437    pars = (60000, 3) # radius, contrast
438
439    ## Set the maximum step size between q values for numerical integration
440    ## This can be based on the biggest dimension in the model, the visible
441    ## oscillation in the data, or None if the data is assumed to be j
442    d_max = pars[0]
443    points_per_fringe = 2000
444    max_step = 2*pi/(d_max*points_per_fringe)
445    #max_step = None
446    #print "max_step",max_step
447
448    ## ==== Select which models are being compared
449    ## Models are smear_MODEL functions above.
450    ## romberg is the slow accurate model computed with adaptive integration
451    ## rectangular is the simplest model to implement, and hopefully good
452    ## enough.  It is the only one that uses max_step at the time of this
453    ## writing, and is what we will use in sasmodels since it seems to be
454    ## good enough.
455    #ms,mt = "simpsons","rectangular"
456    ms,mt = "romberg","rectangular"
457    #ms,mt = "romberg","trapezoid"
458    #ms,mt = "romberg","simpsons"
459
460    ## ==== Evaluate the model and its comparison
461    q_calc = interpolate(q, max_step) if max_step is not None else q
462    y = Iq(q_calc, sphere, pars, scale=scale, background=background)
463    ys = Iq_smeared(q, sphere, pars, delta_qv=delta_qv, max_step=max_step,
464                    scale=scale, background=background, method=ms)
465    yt = Iq_smeared(q, sphere, pars, delta_qv=delta_qv, max_step=max_step,
466                    scale=scale, background=background, method=mt)
467
468    ## ==== plots
469    plt.subplot(211)
470    plt.loglog(q_calc, y, '-', label="unsmeared")
471    plt.loglog(q, ys, '.', label=ms, hold=True)
472    plt.loglog(q, yt, '.', label=mt, hold=True)
473    if answer is not None:
474        plt.loglog(q, answer, '-', label='igor', hold=True)
475    plt.ylabel("y")
476    ## Uncomment to see linear fringing in the model; this works particularly
477    ## well when q_calc is set to extrapolated and interpolated points.
478    #plt.xscale('linear')
479    plt.grid(True)
480    plt.legend()
481    plt.subplot(212)
482    err = (ys-yt)/ys
483    print mt, "err max:", max(abs(err)), " median:", np.median(abs(err))
484    plt.semilogx(q, err, '-', label=mt)
485    ## Uncomment to compare test model to igor
486    #yt=answer # use igor output for comparison
487    #plt.semilogx(q, (ys-yt)/ys, '-', label='Igor', hold=True)
488    plt.xlabel("q")
489    plt.ylabel("(%s - %s)/%s"%(ms,mt,ms))
490    plt.legend()
491
492
493if __name__ == "__main__":
494    demo()
495    import pylab; pylab.show()