source: sasmodels/explore/sinc.py @ e6f1410

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

examine numerical precision of J1c, sinc, j1c and log gamma

  • Property mode set to 100644
File size: 2.2 KB
Line 
1r"""
2Show numerical precision of $sin(x)/x$.
3"""
4
5import numpy as np
6from sympy.mpmath import mp
7#import matplotlib; matplotlib.use('TkAgg')
8import pylab
9
10mp.dps = 150 # number of digits to use in estimating true value
11
12SHOW_DIFF = True  # True if show diff rather than function value
13LINEAR_X = False  # True if q is linearly spaced instead of log spaced
14
15def mp_sinc(vec):
16    """
17    Direct calculation using sympy multiprecision library.
18    """
19    return [_mp_sinc(mp.mpf(x)) for x in vec]
20
21def _mp_sinc(x):
22    """
23    Helper funciton for mp_j1c
24    """
25    return mp.sin(x)/x
26
27def np_sinc(x, dtype):
28    """
29    Direct calculation using scipy.
30    """
31    x = np.asarray(x, dtype)
32    return np.sin(x)/x
33    #return np.asarray(np.sin(np.double(x))/np.double(x),dtype)
34
35def plotdiff(x, target, actual, label):
36    """
37    Plot the computed value.
38
39    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
40    """
41    if SHOW_DIFF:
42        err = np.clip(abs((target-actual)/target), 0, 1)
43        pylab.loglog(x, err, '-', label=label)
44    else:
45        limits = np.min(target), np.max(target)
46        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
47
48def compare(x, precision):
49    r"""
50    Compare the different computation methods using the given precision.
51    """
52    target = np.asarray(mp_sinc(x), 'double')
53    direct = np_sinc(x, precision)
54    plotdiff(x, target, direct, 'direct '+precision)
55    pylab.xlabel("qr (1/Ang)")
56    if SHOW_DIFF:
57        pylab.ylabel("relative error")
58    else:
59        pylab.ylabel("sin(x)/x")
60        pylab.semilogx(x, target,  '-', label="true value")
61    if LINEAR_X:
62        pylab.xscale('linear')
63
64def main():
65    r"""
66    Compare accuracy of different methods for computing $3 j_1(x)/x$.
67    :return:
68    """
69    if LINEAR_X:
70        qr = np.linspace(1,1000,2000)
71    else:
72        qr = np.logspace(-3,5,400)
73    pylab.subplot(121)
74    compare(qr, 'single')
75    pylab.legend(loc='best')
76    pylab.subplot(122)
77    compare(qr, 'double')
78    pylab.legend(loc='best')
79    pylab.suptitle('sin(x)/x')
80
81if __name__ == "__main__":
82    print "\n".join(str(x) for x in mp_sinc([1e-6,1e-5,1e-4,1e-3]))
83    main()
84    pylab.show()
Note: See TracBrowser for help on using the repository browser.