source: sasmodels/explore/precision.py @ 1ceb951

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1ceb951 was 237c9cf, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

flexible_cylinder: improve low q Sdebye calculation; adjust code so different branches use the same equation numbers
  • Property mode set to 100755
File size: 20.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            lin_min, lin_max, lin_steps = 0.001, 2, 2000
84        elif xrange == "log":
85            log_min, log_max, log_steps = -3, 5, 400
86        elif xrange == "logq":
87            log_min, log_max, log_steps = -4, 1, 400
88        else:
89            raise ValueError("unknown range "+xrange)
90        with mp.workprec(500):
91            # Note: we make sure that we are comparing apples to apples...
92            # The x points are set using single precision so that we are
93            # examining the accuracy of the transformation from x to f(x)
94            # rather than x to f(nearest(x)) where nearest(x) is the nearest
95            # value to x in the given precision.
96            if linear:
97                lin_min = max(lin_min, self.limits[0])
98                lin_max = min(lin_max, self.limits[1])
99                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
100                #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double')
101                qr = [mp.mpf(float(v)) for v in qrf]
102                #qr = mp.linspace(lin_min, lin_max, lin_steps)
103            else:
104                log_min = np.log10(max(10**log_min, self.limits[0]))
105                log_max = np.log10(min(10**log_max, self.limits[1]))
106                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
107                #qrf = np.logspace(log_min, log_max, log_steps, dtype='double')
108                qr = [mp.mpf(float(v)) for v in qrf]
109                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
110
111        target = self.call_mpmath(qr, bits=500)
112        pylab.subplot(121)
113        self.compare(qr, 'single', target, linear, diff)
114        pylab.legend(loc='best')
115        pylab.subplot(122)
116        self.compare(qr, 'double', target, linear, diff)
117        pylab.legend(loc='best')
118        pylab.suptitle(self.name + " compared to 500-bit mpmath")
119
120    def compare(self, x, precision, target, linear=False, diff="relative"):
121        r"""
122        Compare the different computation methods using the given precision.
123        """
124        if precision == 'single':
125            #n=11; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
126            #n=23; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
127            pass
128        elif precision == 'double':
129            #n=53; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
130            #n=83; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff)
131            pass
132        plotdiff(x, target, self.call_numpy(x, precision), 'numpy '+precision, diff=diff)
133        plotdiff(x, target, self.call_ocl(x, precision, 0), 'OpenCL '+precision, diff=diff)
134        pylab.xlabel(self.xaxis)
135        if diff == "relative":
136            pylab.ylabel("relative error")
137        elif diff == "absolute":
138            pylab.ylabel("absolute error")
139        else:
140            pylab.ylabel(self.name)
141            pylab.semilogx(x, target, '-', label="true value")
142        if linear:
143            pylab.xscale('linear')
144
145def plotdiff(x, target, actual, label, diff):
146    """
147    Plot the computed value.
148
149    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
150    """
151    if diff == "relative":
152        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
153        #err = np.clip(err, 0, 1)
154        pylab.loglog(x, err, '-', label=label)
155    elif diff == "absolute":
156        err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
157        pylab.loglog(x, err, '-', label=label)
158    else:
159        limits = np.min(target), np.max(target)
160        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
161
162def make_ocl(function, name, source=[]):
163    class Kernel(object):
164        pass
165    Kernel.__file__ = name+".py"
166    Kernel.name = name
167    Kernel.parameters = []
168    Kernel.source = source
169    Kernel.Iq = function
170    model_info = modelinfo.make_model_info(Kernel)
171    return model_info
172
173
174# =============== FUNCTION DEFINITIONS ================
175
176FUNCTIONS = {}
177def add_function(name, mp_function, np_function, ocl_function,
178                 shortname=None, xaxis="x", limits=(-inf, inf)):
179    if shortname is None:
180        shortname = name.replace('(x)', '').replace(' ', '')
181    FUNCTIONS[shortname] = Comparator(name, mp_function, np_function, ocl_function, xaxis, limits)
182
183add_function(
184    name="J0(x)",
185    mp_function=mp.j0,
186    np_function=scipy.special.j0,
187    ocl_function=make_ocl("return sas_J0(q);", "sas_J0", ["lib/polevl.c", "lib/sas_J0.c"]),
188)
189add_function(
190    name="J1(x)",
191    mp_function=mp.j1,
192    np_function=scipy.special.j1,
193    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
194)
195add_function(
196    name="JN(-3, x)",
197    mp_function=lambda x: mp.besselj(-3, x),
198    np_function=lambda x: scipy.special.jn(-3, x),
199    ocl_function=make_ocl("return sas_JN(-3, q);", "sas_JN",
200                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
201    shortname="J-3",
202)
203add_function(
204    name="JN(3, x)",
205    mp_function=lambda x: mp.besselj(3, x),
206    np_function=lambda x: scipy.special.jn(3, x),
207    ocl_function=make_ocl("return sas_JN(3, q);", "sas_JN",
208                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
209    shortname="J3",
210)
211add_function(
212    name="JN(2, x)",
213    mp_function=lambda x: mp.besselj(2, x),
214    np_function=lambda x: scipy.special.jn(2, x),
215    ocl_function=make_ocl("return sas_JN(2, q);", "sas_JN",
216                          ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]),
217    shortname="J2",
218)
219add_function(
220    name="2 J1(x)/x",
221    mp_function=lambda x: 2*mp.j1(x)/x,
222    np_function=lambda x: 2*scipy.special.j1(x)/x,
223    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
224)
225add_function(
226    name="J1(x)",
227    mp_function=mp.j1,
228    np_function=scipy.special.j1,
229    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
230)
231add_function(
232    name="Si(x)",
233    mp_function=mp.si,
234    np_function=lambda x: scipy.special.sici(x)[0],
235    ocl_function=make_ocl("return sas_Si(q);", "sas_Si", ["lib/sas_Si.c"]),
236)
237#import fnlib
238#add_function(
239#    name="fnlibJ1",
240#    mp_function=mp.j1,
241#    np_function=fnlib.J1,
242#    ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]),
243#)
244add_function(
245    name="sin(x)",
246    mp_function=mp.sin,
247    np_function=np.sin,
248    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return sn;", "sas_sin"),
249    ocl_function=make_ocl("return sin(q);", "sas_sin"),
250)
251add_function(
252    name="sin(x)/x",
253    mp_function=lambda x: mp.sin(x)/x if x != 0 else 1,
254    ## scipy sinc function is inaccurate and has an implied pi*x term
255    #np_function=lambda x: scipy.special.sinc(x/pi),
256    ## numpy sin(x)/x needs to check for x=0
257    np_function=lambda x: np.sin(x)/x,
258    ocl_function=make_ocl("return sas_sinx_x(q);", "sas_sinc"),
259)
260add_function(
261    name="cos(x)",
262    mp_function=mp.cos,
263    np_function=np.cos,
264    #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return cn;", "sas_cos"),
265    ocl_function=make_ocl("return cos(q);", "sas_cos"),
266)
267add_function(
268    name="gamma(x)",
269    mp_function=mp.gamma,
270    np_function=scipy.special.gamma,
271    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
272    limits=(-3.1, 10),
273)
274add_function(
275    name="erf(x)",
276    mp_function=mp.erf,
277    np_function=scipy.special.erf,
278    ocl_function=make_ocl("return sas_erf(q);", "sas_erf", ["lib/polevl.c", "lib/sas_erf.c"]),
279    limits=(-5., 5.),
280)
281add_function(
282    name="erfc(x)",
283    mp_function=mp.erfc,
284    np_function=scipy.special.erfc,
285    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
286    limits=(-5., 5.),
287)
288add_function(
289    name="arctan(x)",
290    mp_function=mp.atan,
291    np_function=np.arctan,
292    ocl_function=make_ocl("return atan(q);", "sas_arctan"),
293)
294add_function(
295    name="3 j1(x)/x",
296    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
297    # Note: no taylor expansion near 0
298    np_function=lambda x: 3*(np.sin(x)/x - np.cos(x))/(x*x),
299    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
300)
301add_function(
302    name="(1-cos(x))/x^2",
303    mp_function=lambda x: (1 - mp.cos(x))/(x*x),
304    np_function=lambda x: (1 - np.cos(x))/(x*x),
305    ocl_function=make_ocl("return (1-cos(q))/q/q;", "sas_1mcosx_x2"),
306)
307add_function(
308    name="(1-sin(x)/x)/x",
309    mp_function=lambda x: 1/x - mp.sin(x)/(x*x),
310    np_function=lambda x: 1/x - np.sin(x)/(x*x),
311    ocl_function=make_ocl("return (1-sas_sinx_x(q))/q;", "sas_1msinx_x_x"),
312)
313add_function(
314    name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x",
315    mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x,
316    np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x,
317    ocl_function=make_ocl("return (0.5-sin(q)/q + (1-cos(q))/q/q)/q;", "sas_T2"),
318)
319add_function(
320    name="fmod_2pi",
321    mp_function=lambda x: mp.fmod(x, 2*mp.pi),
322    np_function=lambda x: np.fmod(x, 2*np.pi),
323    ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
324)
325add_function(
326    name="debye",
327    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
328    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
329    ocl_function=make_ocl("""
330    const double qsq = q*q;
331    if (qsq < 1.0) { // Pade approximation
332        const double x = qsq;
333        if (0) { // 0.36 single
334            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}]
335            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1);
336        } else if (0) { // 1.0 for single
337            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
338            const double A1=1./24., A2=1./84, A3=-1./3360;
339            const double B1=3./8., B2=3./56., B3=1./336.;
340            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
341        } else if (1) { // 1.0 for single, 0.25 for double
342            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
343            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
344            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
345            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
346                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
347        } else { // 1.0 for single, 0.5 for double
348            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
349            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
350            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
351            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
352                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
353        }
354    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
355        const double x = qsq;
356        const double C0 = +1.;
357        const double C1 = -1./3.;
358        const double C2 = +1./12.;
359        const double C3 = -1./60.;
360        const double C4 = +1./360.;
361        const double C5 = -1./2520.;
362        const double C6 = +1./20160.;
363        const double C7 = -1./181440.;
364        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
365        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
366        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
367    } else {
368        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
369    }
370    """, "sas_debye"),
371)
372
373RADIUS=3000
374LENGTH=30
375THETA=45
376def mp_cyl(x):
377    f = mp.mpf
378    theta = f(THETA)*mp.pi/f(180)
379    qr = x * f(RADIUS)*mp.sin(theta)
380    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
381    be = f(2)*mp.j1(qr)/qr
382    si = mp.sin(qh)/qh
383    background = f(0)
384    #background = f(1)/f(1000)
385    volume = mp.pi*f(RADIUS)**f(2)*f(LENGTH)
386    contrast = f(5)
387    units = f(1)/f(10000)
388    #return be
389    #return si
390    return units*(volume*contrast*be*si)**f(2)/volume + background
391def np_cyl(x):
392    f = np.float64 if x.dtype == np.float64 else np.float32
393    theta = f(THETA)*f(np.pi)/f(180)
394    qr = x * f(RADIUS)*np.sin(theta)
395    qh = x * f(LENGTH)/f(2)*np.cos(theta)
396    be = f(2)*scipy.special.j1(qr)/qr
397    si = np.sin(qh)/qh
398    background = f(0)
399    #background = f(1)/f(1000)
400    volume = f(np.pi)*f(RADIUS)**2*f(LENGTH)
401    contrast = f(5)
402    units = f(1)/f(10000)
403    #return be
404    #return si
405    return units*(volume*contrast*be*si)**f(2)/volume + background
406ocl_cyl = """\
407    double THETA = %(THETA).15e*M_PI_180;
408    double qr = q*%(RADIUS).15e*sin(THETA);
409    double qh = q*0.5*%(LENGTH).15e*cos(THETA);
410    double be = sas_2J1x_x(qr);
411    double si = sas_sinx_x(qh);
412    double background = 0;
413    //double background = 0.001;
414    double volume = M_PI*square(%(RADIUS).15e)*%(LENGTH).15e;
415    double contrast = 5.0;
416    double units = 1e-4;
417    //return be;
418    //return si;
419    return units*square(volume*contrast*be*si)/volume + background;
420"""%{"LENGTH":LENGTH, "RADIUS": RADIUS, "THETA": THETA}
421add_function(
422    name="cylinder(r=%g, l=%g, theta=%g)"%(RADIUS, LENGTH, THETA),
423    mp_function=mp_cyl,
424    np_function=np_cyl,
425    ocl_function=make_ocl(ocl_cyl, "ocl_cyl", ["lib/polevl.c", "lib/sas_J1.c"]),
426    shortname="cylinder",
427    xaxis="$q/A^{-1}$",
428)
429
430lanczos_gamma = """\
431    const double coeff[] = {
432            76.18009172947146,     -86.50532032941677,
433            24.01409824083091,     -1.231739572450155,
434            0.1208650973866179e-2,-0.5395239384953e-5
435            };
436    const double x = q;
437    double tmp  = x + 5.5;
438    tmp -= (x + 0.5)*log(tmp);
439    double ser = 1.000000000190015;
440    for (int k=0; k < 6; k++) ser += coeff[k]/(x + k+1);
441    return -tmp + log(2.5066282746310005*ser/x);
442"""
443add_function(
444    name="log gamma(x)",
445    mp_function=mp.loggamma,
446    np_function=scipy.special.gammaln,
447    ocl_function=make_ocl(lanczos_gamma, "lgamma"),
448)
449
450# Alternate versions of 3 j1(x)/x, for posterity
451def taylor_3j1x_x(x):
452    """
453    Calculation using taylor series.
454    """
455    # Generate coefficients using the precision of the target value.
456    n = 5
457    cinv = [3991680, -45360, 840, -30, 3]
458    three = x.dtype.type(3)
459    p = three/np.array(cinv, x.dtype)
460    return np.polyval(p[-n:], x*x)
461add_function(
462    name="3 j1(x)/x: taylor",
463    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
464    np_function=taylor_3j1x_x,
465    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
466)
467def trig_3j1x_x(x):
468    r"""
469    Direct calculation using linear combination of sin/cos.
470
471    Use the following trig identity:
472
473    .. math::
474
475        a \sin(x) + b \cos(x) = c \sin(x + \phi)
476
477    where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the
478    numerator $\sin(x) - x\cos(x)$.
479    """
480    one = x.dtype.type(1)
481    three = x.dtype.type(3)
482    c = np.sqrt(one + x*x)
483    phi = np.arctan2(-x, one)
484    return three*(c*np.sin(x+phi))/(x*x*x)
485add_function(
486    name="3 j1(x)/x: trig",
487    mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x),
488    np_function=trig_3j1x_x,
489    ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]),
490)
491def np_2J1x_x(x):
492    """
493    numpy implementation of 2J1(x)/x using single precision algorithm
494    """
495    # pylint: disable=bad-continuation
496    f = x.dtype.type
497    ax = abs(x)
498    if ax < f(8.0):
499        y = x*x
500        ans1 = f(2)*(f(72362614232.0)
501                  + y*(f(-7895059235.0)
502                  + y*(f(242396853.1)
503                  + y*(f(-2972611.439)
504                  + y*(f(15704.48260)
505                  + y*(f(-30.16036606)))))))
506        ans2 = (f(144725228442.0)
507                  + y*(f(2300535178.0)
508                  + y*(f(18583304.74)
509                  + y*(f(99447.43394)
510                  + y*(f(376.9991397)
511                  + y)))))
512        return ans1/ans2
513    else:
514        y = f(64.0)/(ax*ax)
515        xx = ax - f(2.356194491)
516        ans1 = (f(1.0)
517                  + y*(f(0.183105e-2)
518                  + y*(f(-0.3516396496e-4)
519                  + y*(f(0.2457520174e-5)
520                  + y*f(-0.240337019e-6)))))
521        ans2 = (f(0.04687499995)
522                  + y*(f(-0.2002690873e-3)
523                  + y*(f(0.8449199096e-5)
524                  + y*(f(-0.88228987e-6)
525                  + y*f(0.105787412e-6)))))
526        sn, cn = np.sin(xx), np.cos(xx)
527        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
528        return -ans if (x < f(0.0)) else ans
529add_function(
530    name="2 J1(x)/x:alt",
531    mp_function=lambda x: 2*mp.j1(x)/x,
532    np_function=lambda x: np.asarray([np_2J1x_x(v) for v in x], x.dtype),
533    ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
534)
535
536ALL_FUNCTIONS = set(FUNCTIONS.keys())
537ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
538ALL_FUNCTIONS.discard("3j1/x:taylor")
539ALL_FUNCTIONS.discard("3j1/x:trig")
540ALL_FUNCTIONS.discard("2J1/x:alt")
541
542# =============== MAIN PROGRAM ================
543
544def usage():
545    names = ", ".join(sorted(ALL_FUNCTIONS))
546    print("""\
547usage: precision.py [-f/a/r] [-x<range>] name...
548where
549    -f indicates that the function value should be plotted,
550    -a indicates that the absolute error should be plotted,
551    -r indicates that the relative error should be plotted (default),
552    -x<range> indicates the steps in x, where <range> is one of the following
553      log indicates log stepping in [10^-3, 10^5] (default)
554      logq indicates log stepping in [10^-4, 10^1]
555      linear indicates linear stepping in [1, 1000]
556      zoom indicates linear stepping in [1000, 1010]
557      neg indicates linear stepping in [-100.1, 100.1]
558and name is "all [first]" or one of:
559    """+names)
560    sys.exit(1)
561
562def main():
563    import sys
564    diff = "relative"
565    xrange = "log"
566    options = [v for v in sys.argv[1:] if v.startswith('-')]
567    for opt in options:
568        if opt == '-f':
569            diff = "none"
570        elif opt == '-r':
571            diff = "relative"
572        elif opt == '-a':
573            diff = "absolute"
574        elif opt.startswith('-x'):
575            xrange = opt[2:]
576        else:
577            usage()
578
579    names = [v for v in sys.argv[1:] if not v.startswith('-')]
580    if not names:
581        usage()
582
583    if names[0] == "all":
584        cutoff = names[1] if len(names) > 1 else ""
585        names = list(sorted(ALL_FUNCTIONS))
586        names = [k for k in names if k >= cutoff]
587    if any(k not in FUNCTIONS for k in names):
588        usage()
589    multiple = len(names) > 1
590    pylab.interactive(multiple)
591    for k in names:
592        pylab.clf()
593        comparator = FUNCTIONS[k]
594        comparator.run(xrange=xrange, diff=diff)
595        if multiple:
596            raw_input()
597    if not multiple:
598        pylab.show()
599
600if __name__ == "__main__":
601    main()
Note: See TracBrowser for help on using the repository browser.