source: sasmodels/explore/log_gamma.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.8 KB
Line 
1r"""
2Show numerical precision of $sin(x)/x$.
3"""
4
5import numpy as np
6import scipy.special
7from sympy.mpmath import mp
8#import matplotlib; matplotlib.use('TkAgg')
9import pylab
10
11mp.dps = 150 # number of digits to use in estimating true value
12
13SHOW_DIFF = True  # True if show diff rather than function value
14LINEAR_X = False  # True if q is linearly spaced instead of log spaced
15
16def mp_gamma(vec):
17    """
18    Direct calculation using sympy multiprecision library.
19    """
20    return [_mp_fn(mp.mpf(x)) for x in vec]
21
22def _mp_fn(x):
23    """
24    Helper funciton for mp_j1c
25    """
26    #return mp.gamma(x)
27    return mp.loggamma(x)
28
29def np_gamma(x, dtype):
30    """
31    Direct calculation using scipy.
32    """
33    x = np.asarray(x, dtype)
34    return scipy.special.gammaln(x)
35    #return scipy.special.gamma(x)
36
37def lanczos_gamma(x, dtype):
38    coeff = np.asarray((
39         76.18009172947146,     -86.50532032941677,
40         24.01409824083091,     -1.231739572450155,
41          0.1208650973866179e-2,-0.5395239384953e-5), dtype)
42
43    x = np.asarray(x, dtype)
44    tmp  = x + np.asarray(5.5, dtype)
45    tmp -= (x + np.asarray(0.5, dtype))*np.log(tmp)
46    ser  = np.ones_like(x)*np.asarray(1.000000000190015, dtype)
47    for k,c in enumerate(coeff):
48        ser += c/(x + np.asarray(k+1, dtype))
49    return -tmp+np.log(np.asarray(2.5066282746310005, dtype)*ser/x);
50
51def plotdiff(x, target, actual, label):
52    """
53    Plot the computed value.
54
55    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
56    """
57    if SHOW_DIFF:
58        err = np.clip(abs((target-actual)/target), 0, 1)
59        pylab.loglog(x, err, '-', label=label)
60    else:
61        limits = np.min(target), np.max(target)
62        pylab.loglog(x, np.clip(actual,*limits),  '-', label=label)
63
64def compare(x, precision):
65    r"""
66    Compare the different computation methods using the given precision.
67    """
68    target = np.asarray(mp_gamma(x), 'double')
69    direct = np_gamma(x, precision)
70    approx = lanczos_gamma(x, precision)
71    plotdiff(x, target, direct, 'scipy '+precision)
72    plotdiff(x, target, approx, 'sasmodels '+precision)
73    pylab.xlabel("x (arbitrary)")
74    if SHOW_DIFF:
75        pylab.ylabel("relative error")
76    else:
77        pylab.ylabel("ln(gamma(x))")
78        pylab.loglog(x, target,  '-', label="true value")
79    if LINEAR_X:
80        pylab.xscale('linear')
81
82def main():
83    r"""
84    Compare accuracy of different methods for computing $3 j_1(x)/x$.
85    :return:
86    """
87    if LINEAR_X:
88        qr = np.linspace(1,1000,2000)
89    else:
90        qr = np.logspace(-3,5,400)
91    pylab.subplot(121)
92    compare(qr, 'single')
93    pylab.legend(loc='best')
94    pylab.subplot(122)
95    compare(qr, 'double')
96    pylab.legend(loc='best')
97    pylab.suptitle('ln gamma')
98
99if __name__ == "__main__":
100    main()
101    pylab.show()
Note: See TracBrowser for help on using the repository browser.