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