source: sasmodels/explore/J1.py @ 0e0e867

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0e0e867 was 95ce773, checked in by wojciech, 9 years ago

Bessel functions clean-up

  • Property mode set to 100644
File size: 2.9 KB
Line 
1r"""
2Show numerical precision of $2 J_1(x)/x$.
3"""
4import sys; sys.path.insert(0, '..')
5
6import numpy as np
7try:
8    from mpmath import mp
9except:
10    # CRUFT: mpmath split out into its own package
11    from sympy.mpmath import mp
12#import matplotlib; matplotlib.use('TkAgg')
13import pylab
14
15
16SHOW_DIFF = True # True if show diff rather than function value
17#SHOW_DIFF = False # True if show diff rather than function value
18LINEAR_X = False  # True if q is linearly spaced instead of log spaced
19#LINEAR_X = True # True if q is linearly spaced instead of log spaced
20FUNCTION = "J1(x)"
21
22def mp_fn(vec, bits=500):
23    """
24    Direct calculation using sympy multiprecision library.
25    """
26    with mp.workprec(bits):
27        return [_mp_fn(mp.mpf(x)) for x in vec]
28
29def _mp_fn(x):
30    """
31    Actual function that gets evaluated.  The caller just vectorizes.
32    """
33    #return mp.mpf(2)*mp.j1(x)/x
34    return mp.j1(x)
35
36def np_fn(x, dtype):
37    """
38    Direct calculation using scipy.
39    """
40    from scipy.special import j1 as J1
41    x = np.asarray(x, dtype)
42    #return np.asarray(2, dtype)*J1(x)/x
43    return J1(x)
44
45def sasmodels_fn(x, dtype, platform='ocl'):
46    """
47    Calculation using pade approximant.
48    """
49    from sasmodels import core, data, direct_model
50    model = core.load_model('bessel', dtype=dtype)
51    calculator = direct_model.DirectModel(data.empty_data1D(x), model)
52
53    return calculator(background=0)
54
55def plotdiff(x, target, actual, label):
56    """
57    Plot the computed value.
58
59    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
60    """
61    if SHOW_DIFF:
62        err = abs((target-actual)/target)
63        #err = np.clip(err, 0, 1)
64        pylab.loglog(x, err, '-', label=label)
65    else:
66        limits = np.min(target), np.max(target)
67        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
68
69def compare(x, precision, target):
70    r"""
71    Compare the different computation methods using the given precision.
72    """
73    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits')
74    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision)
75    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision)
76    pylab.xlabel("qr (1/Ang)")
77    if SHOW_DIFF:
78        pylab.ylabel("relative error")
79    else:
80        pylab.ylabel(FUNCTION)
81        pylab.semilogx(x, target,  '-', label="true value")
82    if LINEAR_X:
83        pylab.xscale('linear')
84
85def main():
86    r"""
87    Compare accuracy of different methods for computing $3 j_1(x)/x$.
88    :return:
89    """
90    if LINEAR_X:
91        qr = np.linspace(1,1000,2000)
92    else:
93        qr = np.logspace(-3,5,400)
94    target = np.asarray(mp_fn(qr), 'double')
95    pylab.subplot(121)
96    compare(qr, 'single', target)
97    pylab.legend(loc='best')
98    pylab.subplot(122)
99    compare(qr, 'double', target)
100    pylab.legend(loc='best')
101    pylab.suptitle(FUNCTION)
102
103if __name__ == "__main__":
104    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
105    main()
106    pylab.show()
Note: See TracBrowser for help on using the repository browser.