source: sasmodels/explore/precision.py @ 3448301

Last change on this file since 3448301 was cd28947, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

code cleanup and clarification comments
  • Property mode set to 100755
File size: 25.2 KB
Line 
1#!/usr/bin/env python
2r"""
3Show numerical precision of various expressions.
4
5Evaluates the same function(s) in single and double precision and compares
6the results to 500 digit mpmath evaluation of the same function.
7
8Note: a quick way to generation C and python code for taylor series
9expansions from sympy:
10
11    import sympy as sp
12    x = sp.var("x")
13    f = sp.sin(x)/x
14    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term
15    p = sp.horner(t)   # Horner representation
16    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero
17    p = p.n(15)  # evaluate coefficients to 15 digits (optional)
18    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code
19    py_code = c[:-1]  # strip semicolon to convert c to python
20
21    # mpmath has pade() rational function approximation, which might work
22    # better than the taylor series for some functions:
23    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M)
24    P = sum(a*x**n for n,a in enumerate(reversed(P)))
25    Q = sum(a*x**n for n,a in enumerate(reversed(Q)))
26    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p"))
27
28    # There are richardson and shanks series accelerators in both sympy
29    # and mpmath that may be helpful.
30"""
31from __future__ import division, print_function
32
33import sys
34import os
35sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
36
37import numpy as np
38from numpy import pi, inf
39import scipy.special
40try:
41    from mpmath import mp
42except ImportError:
43    # CRUFT: mpmath split out into its own package
44    from sympy.mpmath import mp
45#import matplotlib; matplotlib.use('TkAgg')
46import pylab
47
48from sasmodels import core, data, direct_model, modelinfo
49
50class Comparator(object):
51    def __init__(self, name, mp_function, np_function, ocl_function, xaxis, limits):
52        self.name = name
53        self.mp_function = mp_function
54        self.np_function = np_function
55        self.ocl_function = ocl_function
56        self.xaxis = xaxis
57        self.limits = limits
58
59    def __repr__(self):
60        return "Comparator(%s)"%self.name
61
62    def call_mpmath(self, vec, bits=500):
63        """
64        Direct calculation using mpmath extended precision library.
65        """
66        with mp.workprec(bits):
67            return [self.mp_function(mp.mpf(x)) for x in vec]
68
69    def call_numpy(self, x, dtype):
70        """
71        Direct calculation using numpy/scipy.
72        """
73        x = np.asarray(x, dtype)
74        return self.np_function(x)
75
76    def call_ocl(self, x, dtype, platform='ocl'):
77        """
78        Calculation using sasmodels ocl libraries.
79        """
80        x = np.asarray(x, dtype)
81        model = core.build_model(self.ocl_function, dtype=dtype)
82        calculator = direct_model.DirectModel(data.empty_data1D(x), model)
83        return calculator(background=0)
84
85    def run(self, xrange="log", diff="relative"):
86        r"""
87        Compare accuracy of different methods for computing f.
88
89        *xrange* is::
90
91            log:    [10^-3,10^5]
92            logq:   [10^-4, 10^1]
93            linear: [1,1000]
94            zoom:   [1000,1010]
95            neg:    [-100,100]
96
97        For arbitrary range use "start:stop:steps:scale" where scale is
98        one of log, lin, or linear.
99
100        *diff* is "relative", "absolute" or "none"
101
102        *x_bits* is the precision with which the x values are specified.  The
103        default 23 should reproduce the equivalent of a single precisio
104        """
105        linear = not xrange.startswith("log")
106        if xrange == "zoom":
107            start, stop, steps = 1000, 1010, 2000
108        elif xrange == "neg":
109            start, stop, steps = -100.1, 100.1, 2000
110        elif xrange == "linear":
111            start, stop, steps = 1, 1000, 2000
112            start, stop, steps = 0.001, 2, 2000
113        elif xrange == "log":
114            start, stop, steps = -3, 5, 400
115        elif xrange == "logq":
116            start, stop, steps = -4, 1, 400
117        elif ':' in xrange:
118            parts = xrange.split(':')
119            linear = parts[3] != "log" if len(parts) == 4 else True
120            steps = int(parts[2]) if len(parts) > 2 else 400
121            start = float(parts[0])
122            stop = float(parts[1])
123
124        else:
125            raise ValueError("unknown range "+xrange)
126        with mp.workprec(500):
127            # Note: we make sure that we are comparing apples to apples...
128            # The x points are set using single precision so that we are
129            # examining the accuracy of the transformation from x to f(x)
130            # rather than x to f(nearest(x)) where nearest(x) is the nearest
131            # value to x in the given precision.
132            if linear:
133                start = max(start, self.limits[0])
134                stop = min(stop, self.limits[1])
135                qrf = np.linspace(start, stop, steps, dtype='single')
136                #qrf = np.linspace(start, stop, steps, dtype='double')
137                qr = [mp.mpf(float(v)) for v in qrf]
138                #qr = mp.linspace(start, stop, steps)
139            else:
140                start = np.log10(max(10**start, self.limits[0]))
141                stop = np.log10(min(10**stop, self.limits[1]))
142                qrf = np.logspace(start, stop, steps, dtype='single')
143                #qrf = np.logspace(start, stop, steps, dtype='double')
144                qr = [mp.mpf(float(v)) for v in qrf]
145                #qr = [10**v for v in mp.linspace(start, stop, steps)]
146
147        target = self.call_mpmath(qr, bits=500)
148        pylab.subplot(121)
149        self.compare(qr, 'single', target, linear, diff)
150        pylab.legend(loc='best')
151        pylab.subplot(122)
152        self.compare(qr, 'double', target, linear, diff)
153        pylab.legend(loc='best')
154        pylab.suptitle(self.name + " compared to 500-bit mpmath")
155
156    def compare(self, x, precision, target, linear=False, diff="relative"):
157        r"""
158        Compare the different computation methods using the given precision.
159        """
160        if precision == 'single':
161            #n=11; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
162            #n=23; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
163            pass
164        elif precision == 'double':
165            #n=53; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
166            #n=83; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
167            pass
168        plotdiff(x, target, self.call_numpy(x, precision), 'numpy '+precision, diff=diff)
169        plotdiff(x, target, self.call_ocl(x, precision, 0), 'OpenCL '+precision, diff=diff)
170        pylab.xlabel(self.xaxis)
171        if diff == "relative":
172            pylab.ylabel("relative error")
173        elif diff == "absolute":
174            pylab.ylabel("absolute error")
175        else:
176            pylab.ylabel(self.name)
177            pylab.semilogx(x, target, '-', label="true value")
178        if linear:
179            pylab.xscale('linear')
180
181def plotdiff(x, target, actual, label, diff):
182    """
183    Plot the computed value.
184
185    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
186    """
187    if diff == "relative":
188        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
189        #err = np.clip(err, 0, 1)
190        pylab.loglog(x, err, '-', label=label)
191    elif diff == "absolute":
192        err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
193        pylab.loglog(x, err, '-', label=label)
194    else:
195        limits = np.min(target), np.max(target)
196        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
197
198def make_ocl(function, name, source=[]):
199    class Kernel(object):
200        pass
201    Kernel.__file__ = name+".py"
202    Kernel.name = name
203    Kernel.parameters = []
204    Kernel.source = source
205    Kernel.Iq = function
206    model_info = modelinfo.make_model_info(Kernel)
207    return model_info
208
209# Hack to allow second parameter A in the gammainc and gammaincc functions.
210# Create a 2-D variant of the precision test if we need to handle other two
211# parameter functions.
212A = 1
213def parse_extra_pars():
214    """
215    Parse the command line looking for the second parameter "A=..." for the
216    gammainc/gammaincc functions.
217    """
218    global A
219
220    A_str = str(A)
221    pop = []
222    for k, v in enumerate(sys.argv[1:]):
223        if v.startswith("A="):
224            A_str = v[2:]
225            pop.append(k+1)
226    if pop:
227        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop]
228        A = float(A_str)
229
230parse_extra_pars()
231
232
233# =============== FUNCTION DEFINITIONS ================
234
235FUNCTIONS = {}
236def add_function(name, mp_function, np_function, ocl_function,
237                 shortname=None, xaxis="x", limits=(-inf, inf)):
238    if shortname is None:
239        shortname = name.replace('(x)', '').replace(' ', '')
240    FUNCTIONS[shortname] = Comparator(name, mp_function, np_function, ocl_function, xaxis, limits)
241
242add_function(
243    name="J0(x)",
244    mp_function=mp.j0,
245    np_function=scipy.special.j0,
246    ocl_function=make_ocl("return sas_J0(q);", "sas_J0", ["lib/polevl.c", "lib/sas_J0.c"]),
247)
248add_function(
249    name="J1(x)",
250    mp_function=mp.j1,
251    np_function=scipy.special.j1,
252    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
253)
254add_function(
255    name="JN(-3, x)",
256    mp_function=lambda x: mp.besselj(-3, x),
257    np_function=lambda x: scipy.special.jn(-3, x),
258    ocl_function=make_ocl("return sas_JN(-3, q);", "sas_JN",
259                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
260    shortname="J-3",
261)
262add_function(
263    name="JN(3, x)",
264    mp_function=lambda x: mp.besselj(3, x),
265    np_function=lambda x: scipy.special.jn(3, x),
266    ocl_function=make_ocl("return sas_JN(3, q);", "sas_JN",
267                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
268    shortname="J3",
269)
270add_function(
271    name="JN(2, x)",
272    mp_function=lambda x: mp.besselj(2, x),
273    np_function=lambda x: scipy.special.jn(2, x),
274    ocl_function=make_ocl("return sas_JN(2, q);", "sas_JN",
275                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
276    shortname="J2",
277)
278add_function(
279    name="2 J1(x)/x",
280    mp_function=lambda x: 2*mp.j1(x)/x,
281    np_function=lambda x: 2*scipy.special.j1(x)/x,
282    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
283)
284add_function(
285    name="J1(x)",
286    mp_function=mp.j1,
287    np_function=scipy.special.j1,
288    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
289)
290add_function(
291    name="Si(x)",
292    mp_function=mp.si,
293    np_function=lambda x: scipy.special.sici(x)[0],
294    ocl_function=make_ocl("return sas_Si(q);", "sas_Si", ["lib/sas_Si.c"]),
295)
296#import fnlib
297#add_function(
298#    name="fnlibJ1",
299#    mp_function=mp.j1,
300#    np_function=fnlib.J1,
301#    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
302#)
303add_function(
304    name="sin(x)",
305    mp_function=mp.sin,
306    np_function=np.sin,
307    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return sn;", "sas_sin"),
308    ocl_function=make_ocl("return sin(q);", "sas_sin"),
309)
310add_function(
311    name="sin(x)/x",
312    mp_function=lambda x: mp.sin(x)/x if x != 0 else 1,
313    ## scipy sinc function is inaccurate and has an implied pi*x term
314    #np_function=lambda x: scipy.special.sinc(x/pi),
315    ## numpy sin(x)/x needs to check for x=0
316    np_function=lambda x: np.sin(x)/x,
317    ocl_function=make_ocl("return sas_sinx_x(q);", "sas_sinc"),
318)
319add_function(
320    name="cos(x)",
321    mp_function=mp.cos,
322    np_function=np.cos,
323    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return cn;", "sas_cos"),
324    ocl_function=make_ocl("return cos(q);", "sas_cos"),
325)
326add_function(
327    name="gamma(x)",
328    mp_function=mp.gamma,
329    np_function=scipy.special.gamma,
330    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
331    limits=(-3.1, 10),
332)
333add_function(
334    name="gammaln(x)",
335    mp_function=mp.loggamma,
336    np_function=scipy.special.gammaln,
337    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]),
338    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"),
339)
340add_function(
341    # Note: "a" is given as A=... on the command line via parse_extra_pars
342    name="gammainc(x)",
343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a),
344    np_function=lambda x, a=A: scipy.special.gammainc(a, x),
345    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]),
346)
347add_function(
348    # Note: "a" is given as A=... on the command line via parse_extra_pars
349    name="gammaincc(x)",
350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a),
351    np_function=lambda x, a=A: scipy.special.gammaincc(a, x),
352    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]),
353)
354add_function(
355    name="erf(x)",
356    mp_function=mp.erf,
357    np_function=scipy.special.erf,
358    ocl_function=make_ocl("return sas_erf(q);", "sas_erf", ["lib/polevl.c", "lib/sas_erf.c"]),
359    limits=(-5., 5.),
360)
361add_function(
362    name="erfc(x)",
363    mp_function=mp.erfc,
364    np_function=scipy.special.erfc,
365    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
366    limits=(-5., 5.),
367)
368add_function(
369    name="expm1(x)",
370    mp_function=mp.expm1,
371    np_function=np.expm1,
372    ocl_function=make_ocl("return expm1(q);", "sas_expm1"),
373    limits=(-5., 5.),
374)
375add_function(
376    name="arctan(x)",
377    mp_function=mp.atan,
378    np_function=np.arctan,
379    ocl_function=make_ocl("return atan(q);", "sas_arctan"),
380)
381add_function(
382    name="3 j1(x)/x",
383    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
384    # Note: no taylor expansion near 0
385    np_function=lambda x: 3*(np.sin(x)/x - np.cos(x))/(x*x),
386    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
387)
388add_function(
389    name="(1-cos(x))/x^2",
390    mp_function=lambda x: (1 - mp.cos(x))/(x*x),
391    np_function=lambda x: (1 - np.cos(x))/(x*x),
392    ocl_function=make_ocl("return (1-cos(q))/q/q;", "sas_1mcosx_x2"),
393)
394add_function(
395    name="(1-sin(x)/x)/x",
396    mp_function=lambda x: 1/x - mp.sin(x)/(x*x),
397    np_function=lambda x: 1/x - np.sin(x)/(x*x),
398    ocl_function=make_ocl("return (1-sas_sinx_x(q))/q;", "sas_1msinx_x_x"),
399)
400add_function(
401    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x",
402    mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x,
403    np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x,
404    ocl_function=make_ocl("return (0.5-sin(q)/q + (1-cos(q))/q/q)/q;", "sas_T2"),
405)
406add_function(
407    name="fmod_2pi",
408    mp_function=lambda x: mp.fmod(x, 2*mp.pi),
409    np_function=lambda x: np.fmod(x, 2*np.pi),
410    ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
411)
412add_function(
413    name="gauss_coil",
414    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
415    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
416    ocl_function=make_ocl("""
417    const double qsq = q*q;
418    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
419    // For single: use O(7) Taylor with 0.8 cutoff (7 mad)
420    if (qsq < 0.0) {
421        const double x = qsq;
422        if (0) { // 0.36 single
423            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}]
424            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1);
425        } else if (0) { // 1.0 for single
426            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
427            const double A1=1./24., A2=1./84, A3=-1./3360;
428            const double B1=3./8., B2=3./56., B3=1./336.;
429            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
430        } else if (0) { // 1.0 for single, 0.25 for double
431            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
432            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
433            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
434            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
435                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
436        } else { // 1.0 for single, 0.5 for double
437            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
438            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
439            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
440            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
441                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
442        }
443    } else if (qsq < 0.8) {
444        const double x = qsq;
445        const double C0 = +1.;
446        const double C1 = -1./3.;
447        const double C2 = +1./12.;
448        const double C3 = -1./60.;
449        const double C4 = +1./360.;
450        const double C5 = -1./2520.;
451        const double C6 = +1./20160.;
452        const double C7 = -1./181440.;
453        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
454        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
455        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
456    } else {
457        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
458    }
459    """, "sas_debye"),
460)
461
462RADIUS=3000
463LENGTH=30
464THETA=45
465def mp_cyl(x):
466    f = mp.mpf
467    theta = f(THETA)*mp.pi/f(180)
468    qr = x * f(RADIUS)*mp.sin(theta)
469    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
470    be = f(2)*mp.j1(qr)/qr
471    si = mp.sin(qh)/qh
472    background = f(0)
473    #background = f(1)/f(1000)
474    volume = mp.pi*f(RADIUS)**f(2)*f(LENGTH)
475    contrast = f(5)
476    units = f(1)/f(10000)
477    #return be
478    #return si
479    return units*(volume*contrast*be*si)**f(2)/volume + background
480def np_cyl(x):
481    f = np.float64 if x.dtype == np.float64 else np.float32
482    theta = f(THETA)*f(np.pi)/f(180)
483    qr = x * f(RADIUS)*np.sin(theta)
484    qh = x * f(LENGTH)/f(2)*np.cos(theta)
485    be = f(2)*scipy.special.j1(qr)/qr
486    si = np.sin(qh)/qh
487    background = f(0)
488    #background = f(1)/f(1000)
489    volume = f(np.pi)*f(RADIUS)**2*f(LENGTH)
490    contrast = f(5)
491    units = f(1)/f(10000)
492    #return be
493    #return si
494    return units*(volume*contrast*be*si)**f(2)/volume + background
495ocl_cyl = """\
496    double THETA = %(THETA).15e*M_PI_180;
497    double qr = q*%(RADIUS).15e*sin(THETA);
498    double qh = q*0.5*%(LENGTH).15e*cos(THETA);
499    double be = sas_2J1x_x(qr);
500    double si = sas_sinx_x(qh);
501    double background = 0;
502    //double background = 0.001;
503    double volume = M_PI*square(%(RADIUS).15e)*%(LENGTH).15e;
504    double contrast = 5.0;
505    double units = 1e-4;
506    //return be;
507    //return si;
508    return units*square(volume*contrast*be*si)/volume + background;
509"""%{"LENGTH":LENGTH, "RADIUS": RADIUS, "THETA": THETA}
510add_function(
511    name="cylinder(r=%g, l=%g, theta=%g)"%(RADIUS, LENGTH, THETA),
512    mp_function=mp_cyl,
513    np_function=np_cyl,
514    ocl_function=make_ocl(ocl_cyl, "ocl_cyl", ["lib/polevl.c", "lib/sas_J1.c"]),
515    shortname="cylinder",
516    xaxis="$q/A^{-1}$",
517)
518
519lanczos_gamma = """\
520    const double coeff[] = {
521            76.18009172947146, -86.50532032941677,
522            24.01409824083091, -1.231739572450155,
523            0.1208650973866179e-2,-0.5395239384953e-5
524            };
525    const double x = q;
526    double tmp  = x + 5.5;
527    tmp -= (x + 0.5)*log(tmp);
528    double ser = 1.000000000190015;
529    for (int k=0; k < 6; k++) ser += coeff[k]/(x + k+1);
530    return -tmp + log(2.5066282746310005*ser/x);
531"""
532add_function(
533    name="loggamma(x)",
534    mp_function=mp.loggamma,
535    np_function=scipy.special.gammaln,
536    ocl_function=make_ocl(lanczos_gamma, "lgamma"),
537)
538
539replacement_expm1 = """\
540      double x = (double)q;  // go back to float for single precision kernels
541      // Adapted from the cephes math library.
542      // Copyright 1984 - 1992 by Stephen L. Moshier
543      if (x != x || x == 0.0) {
544         return x; // NaN and +/- 0
545      } else if (x < -0.5 || x > 0.5) {
546         return exp(x) - 1.0;
547      } else {
548         const double xsq = x*x;
549         const double p = (((
550            +1.2617719307481059087798E-4)*xsq
551            +3.0299440770744196129956E-2)*xsq
552            +9.9999999999999999991025E-1);
553         const double q = ((((
554            +3.0019850513866445504159E-6)*xsq
555            +2.5244834034968410419224E-3)*xsq
556            +2.2726554820815502876593E-1)*xsq
557            +2.0000000000000000000897E0);
558         double r = x * p;
559         r =  r / (q - r);
560         return r+r;
561       }
562"""
563add_function(
564    name="sas_expm1(x)",
565    mp_function=mp.expm1,
566    np_function=np.expm1,
567    ocl_function=make_ocl(replacement_expm1, "sas_expm1"),
568)
569
570# Alternate versions of 3 j1(x)/x, for posterity
571def taylor_3j1x_x(x):
572    """
573    Calculation using taylor series.
574    """
575    # Generate coefficients using the precision of the target value.
576    n = 5
577    cinv = [3991680, -45360, 840, -30, 3]
578    three = x.dtype.type(3)
579    p = three/np.array(cinv, x.dtype)
580    return np.polyval(p[-n:], x*x)
581add_function(
582    name="3 j1(x)/x: taylor",
583    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
584    np_function=taylor_3j1x_x,
585    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
586)
587def trig_3j1x_x(x):
588    r"""
589    Direct calculation using linear combination of sin/cos.
590
591    Use the following trig identity:
592
593    .. math::
594
595        a \sin(x) + b \cos(x) = c \sin(x + \phi)
596
597    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
598    numerator $\sin(x) - x\cos(x)$.
599    """
600    one = x.dtype.type(1)
601    three = x.dtype.type(3)
602    c = np.sqrt(one + x*x)
603    phi = np.arctan2(-x, one)
604    return three*(c*np.sin(x+phi))/(x*x*x)
605add_function(
606    name="3 j1(x)/x: trig",
607    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
608    np_function=trig_3j1x_x,
609    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
610)
611def np_2J1x_x(x):
612    """
613    numpy implementation of 2J1(x)/x using single precision algorithm
614    """
615    # pylint: disable=bad-continuation
616    f = x.dtype.type
617    ax = abs(x)
618    if ax < f(8.0):
619        y = x*x
620        ans1 = f(2)*(f(72362614232.0)
621                  + y*(f(-7895059235.0)
622                  + y*(f(242396853.1)
623                  + y*(f(-2972611.439)
624                  + y*(f(15704.48260)
625                  + y*(f(-30.16036606)))))))
626        ans2 = (f(144725228442.0)
627                  + y*(f(2300535178.0)
628                  + y*(f(18583304.74)
629                  + y*(f(99447.43394)
630                  + y*(f(376.9991397)
631                  + y)))))
632        return ans1/ans2
633    else:
634        y = f(64.0)/(ax*ax)
635        xx = ax - f(2.356194491)
636        ans1 = (f(1.0)
637                  + y*(f(0.183105e-2)
638                  + y*(f(-0.3516396496e-4)
639                  + y*(f(0.2457520174e-5)
640                  + y*f(-0.240337019e-6)))))
641        ans2 = (f(0.04687499995)
642                  + y*(f(-0.2002690873e-3)
643                  + y*(f(0.8449199096e-5)
644                  + y*(f(-0.88228987e-6)
645                  + y*f(0.105787412e-6)))))
646        sn, cn = np.sin(xx), np.cos(xx)
647        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
648        return -ans if (x < f(0.0)) else ans
649add_function(
650    name="2 J1(x)/x:alt",
651    mp_function=lambda x: 2*mp.j1(x)/x,
652    np_function=lambda x: np.asarray([np_2J1x_x(v) for v in x], x.dtype),
653    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
654)
655
656ALL_FUNCTIONS = set(FUNCTIONS.keys())
657ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
658ALL_FUNCTIONS.discard("3j1/x:taylor")
659ALL_FUNCTIONS.discard("3j1/x:trig")
660ALL_FUNCTIONS.discard("2J1/x:alt")
661
662# =============== MAIN PROGRAM ================
663
664def usage():
665    names = ", ".join(sorted(ALL_FUNCTIONS))
666    print("""\
667usage: precision.py [-f/a/r] [-x<range>] "name" ...
668where
669    -f indicates that the function value should be plotted,
670    -a indicates that the absolute error should be plotted,
671    -r indicates that the relative error should be plotted (default),
672    -x<range> indicates the steps in x, where <range> is one of the following
673        log indicates log stepping in [10^-3, 10^5] (default)
674        logq indicates log stepping in [10^-4, 10^1]
675        linear indicates linear stepping in [1, 1000]
676        zoom indicates linear stepping in [1000, 1010]
677        neg indicates linear stepping in [-100.1, 100.1]
678        start:stop:n[:stepping] indicates an n-step plot in [start, stop]
679            or [10^start, 10^stop] if stepping is "log" (default n=400)
680Some functions (notably gammainc/gammaincc) have an additional parameter A
681which can be set from the command line as A=value.  Default is A=1.
682
683Name is one of:
684    """+names)
685    sys.exit(1)
686
687def main():
688    import sys
689    diff = "relative"
690    xrange = "log"
691    options = [v for v in sys.argv[1:] if v.startswith('-')]
692    for opt in options:
693        if opt == '-f':
694            diff = "none"
695        elif opt == '-r':
696            diff = "relative"
697        elif opt == '-a':
698            diff = "absolute"
699        elif opt.startswith('-x'):
700            xrange = opt[2:]
701        else:
702            usage()
703
704    names = [v for v in sys.argv[1:] if not v.startswith('-')]
705    if not names:
706        usage()
707
708    if names[0] == "all":
709        cutoff = names[1] if len(names) > 1 else ""
710        names = list(sorted(ALL_FUNCTIONS))
711        names = [k for k in names if k >= cutoff]
712    if any(k not in FUNCTIONS for k in names):
713        usage()
714    multiple = len(names) > 1
715    pylab.interactive(multiple)
716    for k in names:
717        pylab.clf()
718        comparator = FUNCTIONS[k]
719        comparator.run(xrange=xrange, diff=diff)
720        if multiple:
721            raw_input()
722    if not multiple:
723        pylab.show()
724
725if __name__ == "__main__":
726    main()
Note: See TracBrowser for help on using the repository browser.