source: sasmodels/explore/precision.py @ 5181ccc

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

improve accuracy of bessel function J1
  • Property mode set to 100755
File size: 17.8 KB
Line 
1#!/usr/bin/env python
2r"""
3Show numerical precision of $2 J_1(x)/x$.
4"""
5from __future__ import division, print_function
6
7import sys
8import os
9sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
10
11import numpy as np
12from numpy import pi, inf
13import scipy.special
14try:
15    from mpmath import mp
16except ImportError:
17    # CRUFT: mpmath split out into its own package
18    from sympy.mpmath import mp
19#import matplotlib; matplotlib.use('TkAgg')
20import pylab
21
22from sasmodels import core, data, direct_model, modelinfo
23
24class Comparator(object):
25    def __init__(self, name, mp_function, np_function, ocl_function, xaxis, limits):
26        self.name = name
27        self.mp_function = mp_function
28        self.np_function = np_function
29        self.ocl_function = ocl_function
30        self.xaxis = xaxis
31        self.limits = limits
32
33    def __repr__(self):
34        return "Comparator(%s)"%self.name
35
36    def call_mpmath(self, vec, bits=500):
37        """
38        Direct calculation using mpmath extended precision library.
39        """
40        with mp.workprec(bits):
41            return [self.mp_function(mp.mpf(x)) for x in vec]
42
43    def call_numpy(self, x, dtype):
44        """
45        Direct calculation using numpy/scipy.
46        """
47        x = np.asarray(x, dtype)
48        return self.np_function(x)
49
50    def call_ocl(self, x, dtype, platform='ocl'):
51        """
52        Calculation using sasmodels ocl libraries.
53        """
54        x = np.asarray(x, dtype)
55        model = core.build_model(self.ocl_function, dtype=dtype)
56        calculator = direct_model.DirectModel(data.empty_data1D(x), model)
57        return calculator(background=0)
58
59    def run(self, xrange="log", diff="relative"):
60        r"""
61        Compare accuracy of different methods for computing f.
62
63        *xrange* is::
64
65            log:    [10^-3,10^5]
66            logq:   [10^-4, 10^1]
67            linear: [1,1000]
68            zoom:   [1000,1010]
69            neg:    [-100,100]
70
71        *diff* is "relative", "absolute" or "none"
72
73        *x_bits* is the precision with which the x values are specified.  The
74        default 23 should reproduce the equivalent of a single precisio
75        """
76        linear = not xrange.startswith("log")
77        if xrange == "zoom":
78            lin_min, lin_max, lin_steps = 1000, 1010, 2000
79        elif xrange == "neg":
80            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
81        elif xrange == "linear":
82            lin_min, lin_max, lin_steps = 1, 1000, 2000
83        elif xrange == "log":
84            log_min, log_max, log_steps = -3, 5, 400
85        elif xrange == "logq":
86            log_min, log_max, log_steps = -4, 1, 400
87        else:
88            raise ValueError("unknown range "+xrange)
89        with mp.workprec(500):
90            # Note: we make sure that we are comparing apples to apples...
91            # The x points are set using single precision so that we are
92            # examining the accuracy of the transformation from x to f(x)
93            # rather than x to f(nearest(x)) where nearest(x) is the nearest
94            # value to x in the given precision.
95            if linear:
96                lin_min = max(lin_min, self.limits[0])
97                lin_max = min(lin_max, self.limits[1])
98                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
99                #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double')
100                qr = [mp.mpf(float(v)) for v in qrf]
101                #qr = mp.linspace(lin_min, lin_max, lin_steps)
102            else:
103                log_min = np.log10(max(10**log_min, self.limits[0]))
104                log_max = np.log10(min(10**log_max, self.limits[1]))
105                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
106                #qrf = np.logspace(log_min, log_max, log_steps, dtype='double')
107                qr = [mp.mpf(float(v)) for v in qrf]
108                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
109
110        target = self.call_mpmath(qr, bits=500)
111        pylab.subplot(121)
112        self.compare(qr, 'single', target, linear, diff)
113        pylab.legend(loc='best')
114        pylab.subplot(122)
115        self.compare(qr, 'double', target, linear, diff)
116        pylab.legend(loc='best')
117        pylab.suptitle(self.name + " compared to 500-bit mpmath")
118
119    def compare(self, x, precision, target, linear=False, diff="relative"):
120        r"""
121        Compare the different computation methods using the given precision.
122        """
123        if precision == 'single':
124            #n=11; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
125            #n=23; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
126            pass
127        elif precision == 'double':
128            #n=53; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
129            #n=83; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
130            pass
131        plotdiff(x, target, self.call_numpy(x, precision), 'numpy '+precision, diff=diff)
132        plotdiff(x, target, self.call_ocl(x, precision, 0), 'OpenCL '+precision, diff=diff)
133        pylab.xlabel(self.xaxis)
134        if diff == "relative":
135            pylab.ylabel("relative error")
136        elif diff == "absolute":
137            pylab.ylabel("absolute error")
138        else:
139            pylab.ylabel(self.name)
140            pylab.semilogx(x, target, '-', label="true value")
141        if linear:
142            pylab.xscale('linear')
143
144def plotdiff(x, target, actual, label, diff):
145    """
146    Plot the computed value.
147
148    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
149    """
150    if diff == "relative":
151        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
152        #err = np.clip(err, 0, 1)
153        pylab.loglog(x, err, '-', label=label)
154    elif diff == "absolute":
155        err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
156        pylab.loglog(x, err, '-', label=label)
157    else:
158        limits = np.min(target), np.max(target)
159        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
160
161def make_ocl(function, name, source=[]):
162    class Kernel(object):
163        pass
164    Kernel.__file__ = name+".py"
165    Kernel.name = name
166    Kernel.parameters = []
167    Kernel.source = source
168    Kernel.Iq = function
169    model_info = modelinfo.make_model_info(Kernel)
170    return model_info
171
172
173# =============== FUNCTION DEFINITIONS ================
174
175FUNCTIONS = {}
176def add_function(name, mp_function, np_function, ocl_function,
177                 shortname=None, xaxis="x", limits=(-inf, inf)):
178    if shortname is None:
179        shortname = name.replace('(x)', '').replace(' ', '')
180    FUNCTIONS[shortname] = Comparator(name, mp_function, np_function, ocl_function, xaxis, limits)
181
182add_function(
183    name="J0(x)",
184    mp_function=mp.j0,
185    np_function=scipy.special.j0,
186    ocl_function=make_ocl("return sas_J0(q);", "sas_J0", ["lib/polevl.c", "lib/sas_J0.c"]),
187)
188add_function(
189    name="J1(x)",
190    mp_function=mp.j1,
191    np_function=scipy.special.j1,
192    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
193)
194add_function(
195    name="JN(-3, x)",
196    mp_function=lambda x: mp.besselj(-3, x),
197    np_function=lambda x: scipy.special.jn(-3, x),
198    ocl_function=make_ocl("return sas_JN(-3, q);", "sas_JN",
199                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
200    shortname="J-3",
201)
202add_function(
203    name="JN(3, x)",
204    mp_function=lambda x: mp.besselj(3, x),
205    np_function=lambda x: scipy.special.jn(3, x),
206    ocl_function=make_ocl("return sas_JN(3, q);", "sas_JN",
207                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
208    shortname="J3",
209)
210add_function(
211    name="JN(2, x)",
212    mp_function=lambda x: mp.besselj(2, x),
213    np_function=lambda x: scipy.special.jn(2, x),
214    ocl_function=make_ocl("return sas_JN(2, q);", "sas_JN",
215                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
216    shortname="J2",
217)
218add_function(
219    name="2 J1(x)/x",
220    mp_function=lambda x: 2*mp.j1(x)/x,
221    np_function=lambda x: 2*scipy.special.j1(x)/x,
222    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
223)
224add_function(
225    name="J1(x)",
226    mp_function=mp.j1,
227    np_function=scipy.special.j1,
228    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
229)
230add_function(
231    name="Si(x)",
232    mp_function=mp.si,
233    np_function=lambda x: scipy.special.sici(x)[0],
234    ocl_function=make_ocl("return sas_Si(q);", "sas_Si", ["lib/sas_Si.c"]),
235)
236#import fnlib
237#add_function(
238#    name="fnlibJ1",
239#    mp_function=mp.j1,
240#    np_function=fnlib.J1,
241#    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
242#)
243add_function(
244    name="sin(x)",
245    mp_function=mp.sin,
246    np_function=np.sin,
247    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return sn;", "sas_sin"),
248    ocl_function=make_ocl("return sin(q);", "sas_sin"),
249)
250add_function(
251    name="sin(x)/x",
252    mp_function=lambda x: mp.sin(x)/x if x != 0 else 1,
253    ## scipy sinc function is inaccurate and has an implied pi*x term
254    #np_function=lambda x: scipy.special.sinc(x/pi),
255    ## numpy sin(x)/x needs to check for x=0
256    np_function=lambda x: np.sin(x)/x,
257    ocl_function=make_ocl("return sas_sinx_x(q);", "sas_sinc"),
258)
259add_function(
260    name="cos(x)",
261    mp_function=mp.cos,
262    np_function=np.cos,
263    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return cn;", "sas_cos"),
264    ocl_function=make_ocl("return cos(q);", "sas_cos"),
265)
266add_function(
267    name="gamma(x)",
268    mp_function=mp.gamma,
269    np_function=scipy.special.gamma,
270    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
271    limits=(-3.1,10),
272)
273add_function(
274    name="erf(x)",
275    mp_function=mp.erf,
276    np_function=scipy.special.erf,
277    ocl_function=make_ocl("return sas_erf(q);", "sas_erf", ["lib/polevl.c", "lib/sas_erf.c"]),
278    limits=(-5.,5.),
279)
280add_function(
281    name="erfc(x)",
282    mp_function=mp.erfc,
283    np_function=scipy.special.erfc,
284    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
285    limits=(-5.,5.),
286)
287add_function(
288    name="arctan(x)",
289    mp_function=mp.atan,
290    np_function=np.arctan,
291    ocl_function=make_ocl("return atan(q);", "sas_arctan"),
292)
293add_function(
294    name="3 j1(x)/x",
295    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
296    # Note: no taylor expansion near 0
297    np_function=lambda x: 3*(np.sin(x)/x - np.cos(x))/(x*x),
298    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
299)
300add_function(
301    name="fmod_2pi",
302    mp_function=lambda x: mp.fmod(x, 2*mp.pi),
303    np_function=lambda x: np.fmod(x, 2*np.pi),
304    ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
305)
306
307RADIUS=3000
308LENGTH=30
309THETA=45
310def mp_cyl(x):
311    f = mp.mpf
312    theta = f(THETA)*mp.pi/f(180)
313    qr = x * f(RADIUS)*mp.sin(theta)
314    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
315    be = f(2)*mp.j1(qr)/qr
316    si = mp.sin(qh)/qh
317    background = f(0)
318    #background = f(1)/f(1000)
319    volume = mp.pi*f(RADIUS)**f(2)*f(LENGTH)
320    contrast = f(5)
321    units = f(1)/f(10000)
322    #return be
323    #return si
324    return units*(volume*contrast*be*si)**f(2)/volume + background
325def np_cyl(x):
326    f = np.float64 if x.dtype == np.float64 else np.float32
327    theta = f(THETA)*f(np.pi)/f(180)
328    qr = x * f(RADIUS)*np.sin(theta)
329    qh = x * f(LENGTH)/f(2)*np.cos(theta)
330    be = f(2)*scipy.special.j1(qr)/qr
331    si = np.sin(qh)/qh
332    background = f(0)
333    #background = f(1)/f(1000)
334    volume = f(np.pi)*f(RADIUS)**2*f(LENGTH)
335    contrast = f(5)
336    units = f(1)/f(10000)
337    #return be
338    #return si
339    return units*(volume*contrast*be*si)**f(2)/volume + background
340ocl_cyl = """\
341    double THETA = %(THETA).15e*M_PI_180;
342    double qr = q*%(RADIUS).15e*sin(THETA);
343    double qh = q*0.5*%(LENGTH).15e*cos(THETA);
344    double be = sas_2J1x_x(qr);
345    double si = sas_sinx_x(qh);
346    double background = 0;
347    //double background = 0.001;
348    double volume = M_PI*square(%(RADIUS).15e)*%(LENGTH).15e;
349    double contrast = 5.0;
350    double units = 1e-4;
351    //return be;
352    //return si;
353    return units*square(volume*contrast*be*si)/volume + background;
354"""%{"LENGTH":LENGTH, "RADIUS": RADIUS, "THETA": THETA}
355add_function(
356    name="cylinder(r=%g, l=%g, theta=%g)"%(RADIUS, LENGTH, THETA),
357    mp_function=mp_cyl,
358    np_function=np_cyl,
359    ocl_function=make_ocl(ocl_cyl, "ocl_cyl", ["lib/polevl.c", "lib/sas_J1.c"]),
360    shortname="cylinder",
361    xaxis="$q/A^{-1}$",
362)
363
364lanczos_gamma = """\
365    const double coeff[] = {
366            76.18009172947146,     -86.50532032941677,
367            24.01409824083091,     -1.231739572450155,
368            0.1208650973866179e-2,-0.5395239384953e-5
369            };
370    const double x = q;
371    double tmp  = x + 5.5;
372    tmp -= (x + 0.5)*log(tmp);
373    double ser = 1.000000000190015;
374    for (int k=0; k < 6; k++) ser += coeff[k]/(x + k+1);
375    return -tmp + log(2.5066282746310005*ser/x);
376"""
377add_function(
378    name="log gamma(x)",
379    mp_function=mp.loggamma,
380    np_function=scipy.special.gammaln,
381    ocl_function=make_ocl(lanczos_gamma, "lgamma"),
382)
383
384# Alternate versions of 3 j1(x)/x, for posterity
385def taylor_3j1x_x(x):
386    """
387    Calculation using taylor series.
388    """
389    # Generate coefficients using the precision of the target value.
390    n = 5
391    cinv = [3991680, -45360, 840, -30, 3]
392    three = x.dtype.type(3)
393    p = three/np.array(cinv, x.dtype)
394    return np.polyval(p[-n:], x*x)
395add_function(
396    name="3 j1(x)/x: taylor",
397    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
398    np_function=taylor_3j1x_x,
399    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
400)
401def trig_3j1x_x(x):
402    r"""
403    Direct calculation using linear combination of sin/cos.
404
405    Use the following trig identity:
406
407    .. math::
408
409        a \sin(x) + b \cos(x) = c \sin(x + \phi)
410
411    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
412    numerator $\sin(x) - x\cos(x)$.
413    """
414    one = x.dtype.type(1)
415    three = x.dtype.type(3)
416    c = np.sqrt(one + x*x)
417    phi = np.arctan2(-x, one)
418    return three*(c*np.sin(x+phi))/(x*x*x)
419add_function(
420    name="3 j1(x)/x: trig",
421    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
422    np_function=trig_3j1x_x,
423    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
424)
425def np_2J1x_x(x):
426    """
427    numpy implementation of 2J1(x)/x using single precision algorithm
428    """
429    # pylint: disable=bad-continuation
430    f = x.dtype.type
431    ax = abs(x)
432    if ax < f(8.0):
433        y = x*x
434        ans1 = f(2)*(f(72362614232.0)
435                  + y*(f(-7895059235.0)
436                  + y*(f(242396853.1)
437                  + y*(f(-2972611.439)
438                  + y*(f(15704.48260)
439                  + y*(f(-30.16036606)))))))
440        ans2 = (f(144725228442.0)
441                  + y*(f(2300535178.0)
442                  + y*(f(18583304.74)
443                  + y*(f(99447.43394)
444                  + y*(f(376.9991397)
445                  + y)))))
446        return ans1/ans2
447    else:
448        y = f(64.0)/(ax*ax)
449        xx = ax - f(2.356194491)
450        ans1 = (f(1.0)
451                  + y*(f(0.183105e-2)
452                  + y*(f(-0.3516396496e-4)
453                  + y*(f(0.2457520174e-5)
454                  + y*f(-0.240337019e-6)))))
455        ans2 = (f(0.04687499995)
456                  + y*(f(-0.2002690873e-3)
457                  + y*(f(0.8449199096e-5)
458                  + y*(f(-0.88228987e-6)
459                  + y*f(0.105787412e-6)))))
460        sn, cn = np.sin(xx), np.cos(xx)
461        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
462        return -ans if (x < f(0.0)) else ans
463add_function(
464    name="2 J1(x)/x:alt",
465    mp_function=lambda x: 2*mp.j1(x)/x,
466    np_function=lambda x: np.asarray([np_2J1x_x(v) for v in x], x.dtype),
467    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
468)
469
470ALL_FUNCTIONS = set(FUNCTIONS.keys())
471ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
472ALL_FUNCTIONS.discard("3j1/x:taylor")
473ALL_FUNCTIONS.discard("3j1/x:trig")
474ALL_FUNCTIONS.discard("2J1/x:alt")
475
476# =============== MAIN PROGRAM ================
477
478def usage():
479    names = ", ".join(sorted(ALL_FUNCTIONS))
480    print("""\
481usage: precision.py [-f/a/r] [-x<range>] name...
482where
483    -f indicates that the function value should be plotted,
484    -a indicates that the absolute error should be plotted,
485    -r indicates that the relative error should be plotted (default),
486    -x<range> indicates the steps in x, where <range> is one of the following
487      log indicates log stepping in [10^-3, 10^5] (default)
488      logq indicates log stepping in [10^-4, 10^1]
489      linear indicates linear stepping in [1, 1000]
490      zoom indicates linear stepping in [1000, 1010]
491      neg indicates linear stepping in [-100.1, 100.1]
492and name is "all [first]" or one of:
493    """+names)
494    sys.exit(1)
495
496def main():
497    import sys
498    diff = "relative"
499    xrange = "log"
500    options = [v for v in sys.argv[1:] if v.startswith('-')]
501    for opt in options:
502        if opt == '-f':
503            diff = "none"
504        elif opt == '-r':
505            diff = "relative"
506        elif opt == '-a':
507            diff = "absolute"
508        elif opt.startswith('-x'):
509            xrange = opt[2:]
510        else:
511            usage()
512
513    names = [v for v in sys.argv[1:] if not v.startswith('-')]
514    if not names:
515        usage()
516
517    if names[0] == "all":
518        cutoff = names[1] if len(names) > 1 else ""
519        names = list(sorted(ALL_FUNCTIONS))
520        names = [k for k in names if k >= cutoff]
521    if any(k not in FUNCTIONS for k in names):
522        usage()
523    multiple = len(names) > 1
524    pylab.interactive(multiple)
525    for k in names:
526        pylab.clf()
527        comparator = FUNCTIONS[k]
528        comparator.run(xrange=xrange, diff=diff)
529        if multiple:
530            raw_input()
531    if not multiple:
532        pylab.show()
533
534if __name__ == "__main__":
535    main()
Note: See TracBrowser for help on using the repository browser.