source: sasmodels/explore/J1c.py @ 0a6da3c

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0a6da3c was 0a6da3c, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

explore numerical precision of cylinder form

  • Property mode set to 100644
File size: 5.0 KB
Line 
1r"""
2Show numerical precision of $2 J_1(x)/x$.
3"""
4
5import numpy as np
6from sympy.mpmath import mp
7#import matplotlib; matplotlib.use('TkAgg')
8import pylab
9
10
11SHOW_DIFF = True  # True if show diff rather than function value
12LINEAR_X = False  # True if q is linearly spaced instead of log spaced
13
14def mp_J1c(vec, bits=500):
15    """
16    Direct calculation using sympy multiprecision library.
17    """
18    with mp.workprec(bits):
19        return [_mp_J1c(mp.mpf(x)) for x in vec]
20
21def _mp_J1c(x):
22    """
23    Helper funciton for mp_j1c
24    """
25    return mp.mpf(2)*mp.j1(x)/x
26
27def np_J1c(x, dtype):
28    """
29    Direct calculation using scipy.
30    """
31    from scipy.special import j1 as J1
32    x = np.asarray(x, dtype)
33    return np.asarray(2, dtype)*J1(x)/x
34
35def cephes_J1c(x, dtype, n):
36    """
37    Calculation using pade approximant.
38    """
39    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
40    x = np.asarray(x, dtype)
41    ans = np.empty_like(x)
42    ax = abs(x)
43
44    # Branch a
45    a_idx = ax < f(8.0)
46    a_xsq = x[a_idx]**2
47    a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606)))
48    a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0)))
49    a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq)
50    a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq)
51    ans[a_idx] = f(2.0)*a_ans1/a_ans2
52
53    # Branch b
54    b_idx = ~a_idx
55    b_ax = ax[b_idx]
56    b_x = x[b_idx]
57
58    b_y = f(64.0)/(b_ax**2)
59    b_xx = b_ax - f(2.356194491)
60    b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6)))
61    b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)))
62    b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y)
63    b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y)
64    b_sn, b_cn = np.sin(b_xx), np.cos(b_xx)
65    ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x
66
67    return ans
68
69def div_J1c(x, dtype):
70    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
71    x = np.asarray(x, dtype)
72    return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype)
73
74def _J1(x, f):
75    ax = abs(x)
76    if ax < f(8.0):
77        y = x*x
78        ans1 = x*(f(72362614232.0)
79                  + y*(f(-7895059235.0)
80                  + y*(f(242396853.1)
81                  + y*(f(-2972611.439)
82                  + y*(f(15704.48260)
83                  + y*(f(-30.16036606)))))))
84        ans2 = (f(144725228442.0)
85                  + y*(f(2300535178.0)
86                  + y*(f(18583304.74)
87                  + y*(f(99447.43394)
88                  + y*(f(376.9991397)
89                  + y)))))
90        return ans1/ans2
91    else:
92        y = f(64.0)/(ax*ax)
93        xx = ax - f(2.356194491)
94        ans1 = (f(1.0)
95                  + y*(f(0.183105e-2)
96                  + y*(f(-0.3516396496e-4)
97                  + y*(f(0.2457520174e-5)
98                  + y*f(-0.240337019e-6)))))
99        ans2 = (f(0.04687499995)
100                  + y*(f(-0.2002690873e-3)
101                  + y*(f(0.8449199096e-5)
102                  + y*(f(-0.88228987e-6)
103                  + y*f(0.105787412e-6)))))
104        sn, cn = np.sin(xx), np.cos(xx)
105        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2)
106        return -ans if (x < f(0.0)) else ans
107
108def plotdiff(x, target, actual, label):
109    """
110    Plot the computed value.
111
112    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
113    """
114    if SHOW_DIFF:
115        err = np.clip(abs((target-actual)/target), 0, 1)
116        pylab.loglog(x, err, '-', label=label)
117    else:
118        limits = np.min(target), np.max(target)
119        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
120
121def compare(x, precision):
122    r"""
123    Compare the different computation methods using the given precision.
124    """
125    target = np.asarray(mp_J1c(x), 'double')
126    #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits')
127    plotdiff(x, target, np_J1c(x, precision), 'direct '+precision)
128    plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision)
129    #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision)
130    #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision)
131    pylab.xlabel("qr (1/Ang)")
132    if SHOW_DIFF:
133        pylab.ylabel("relative error")
134    else:
135        pylab.ylabel("2 J1(x)/x")
136        pylab.semilogx(x, target,  '-', label="true value")
137    if LINEAR_X:
138        pylab.xscale('linear')
139
140def main():
141    r"""
142    Compare accuracy of different methods for computing $3 j_1(x)/x$.
143    :return:
144    """
145    if LINEAR_X:
146        qr = np.linspace(1,1000,2000)
147    else:
148        qr = np.logspace(-3,5,400)
149    pylab.subplot(121)
150    compare(qr, 'single')
151    pylab.legend(loc='best')
152    pylab.subplot(122)
153    compare(qr, 'double')
154    pylab.legend(loc='best')
155    pylab.suptitle('2 J1(x)/x')
156
157if __name__ == "__main__":
158    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
159    main()
160    pylab.show()
Note: See TracBrowser for help on using the repository browser.