| 1 | r""" |
|---|
| 2 | Show numerical precision of $sin(x)/x$. |
|---|
| 3 | """ |
|---|
| 4 | |
|---|
| 5 | import numpy as np |
|---|
| 6 | import scipy.special |
|---|
| 7 | from sympy.mpmath import mp |
|---|
| 8 | #import matplotlib; matplotlib.use('TkAgg') |
|---|
| 9 | import pylab |
|---|
| 10 | |
|---|
| 11 | mp.dps = 150 # number of digits to use in estimating true value |
|---|
| 12 | |
|---|
| 13 | SHOW_DIFF = True # True if show diff rather than function value |
|---|
| 14 | LINEAR_X = False # True if q is linearly spaced instead of log spaced |
|---|
| 15 | |
|---|
| 16 | def mp_gamma(vec): |
|---|
| 17 | """ |
|---|
| 18 | Direct calculation using sympy multiprecision library. |
|---|
| 19 | """ |
|---|
| 20 | return [_mp_fn(mp.mpf(x)) for x in vec] |
|---|
| 21 | |
|---|
| 22 | def _mp_fn(x): |
|---|
| 23 | """ |
|---|
| 24 | Helper funciton for mp_j1c |
|---|
| 25 | """ |
|---|
| 26 | #return mp.gamma(x) |
|---|
| 27 | return mp.loggamma(x) |
|---|
| 28 | |
|---|
| 29 | def 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 | |
|---|
| 37 | def 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 | |
|---|
| 51 | def 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 | |
|---|
| 64 | def 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 | |
|---|
| 82 | def 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 | |
|---|
| 99 | if __name__ == "__main__": |
|---|
| 100 | main() |
|---|
| 101 | pylab.show() |
|---|
