source: sasmodels/explore/precision.py @ 57eb6a4

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

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