| 1 | #!/usr/bin/env python |
|---|
| 2 | r""" |
|---|
| 3 | Show numerical precision of $2 J_1(x)/x$. |
|---|
| 4 | """ |
|---|
| 5 | from __future__ import division, print_function |
|---|
| 6 | |
|---|
| 7 | import sys |
|---|
| 8 | import os |
|---|
| 9 | sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))) |
|---|
| 10 | |
|---|
| 11 | import numpy as np |
|---|
| 12 | from numpy import pi, inf |
|---|
| 13 | import scipy.special |
|---|
| 14 | try: |
|---|
| 15 | from mpmath import mp |
|---|
| 16 | except ImportError: |
|---|
| 17 | # CRUFT: mpmath split out into its own package |
|---|
| 18 | from sympy.mpmath import mp |
|---|
| 19 | #import matplotlib; matplotlib.use('TkAgg') |
|---|
| 20 | import pylab |
|---|
| 21 | |
|---|
| 22 | from sasmodels import core, data, direct_model, modelinfo |
|---|
| 23 | |
|---|
| 24 | class Comparator(object): |
|---|
| 25 | def __init__(self, name, mp_function, np_function, ocl_function, xaxis, limits): |
|---|
| 26 | self.name = name |
|---|
| 27 | self.mp_function = mp_function |
|---|
| 28 | self.np_function = np_function |
|---|
| 29 | self.ocl_function = ocl_function |
|---|
| 30 | self.xaxis = xaxis |
|---|
| 31 | self.limits = limits |
|---|
| 32 | |
|---|
| 33 | def __repr__(self): |
|---|
| 34 | return "Comparator(%s)"%self.name |
|---|
| 35 | |
|---|
| 36 | def call_mpmath(self, vec, bits=500): |
|---|
| 37 | """ |
|---|
| 38 | Direct calculation using mpmath extended precision library. |
|---|
| 39 | """ |
|---|
| 40 | with mp.workprec(bits): |
|---|
| 41 | return [self.mp_function(mp.mpf(x)) for x in vec] |
|---|
| 42 | |
|---|
| 43 | def call_numpy(self, x, dtype): |
|---|
| 44 | """ |
|---|
| 45 | Direct calculation using numpy/scipy. |
|---|
| 46 | """ |
|---|
| 47 | x = np.asarray(x, dtype) |
|---|
| 48 | return self.np_function(x) |
|---|
| 49 | |
|---|
| 50 | def call_ocl(self, x, dtype, platform='ocl'): |
|---|
| 51 | """ |
|---|
| 52 | Calculation using sasmodels ocl libraries. |
|---|
| 53 | """ |
|---|
| 54 | x = np.asarray(x, dtype) |
|---|
| 55 | model = core.build_model(self.ocl_function, dtype=dtype) |
|---|
| 56 | calculator = direct_model.DirectModel(data.empty_data1D(x), model) |
|---|
| 57 | return calculator(background=0) |
|---|
| 58 | |
|---|
| 59 | def run(self, xrange="log", diff="relative"): |
|---|
| 60 | r""" |
|---|
| 61 | Compare accuracy of different methods for computing f. |
|---|
| 62 | |
|---|
| 63 | *xrange* is:: |
|---|
| 64 | |
|---|
| 65 | log: [10^-3,10^5] |
|---|
| 66 | logq: [10^-4, 10^1] |
|---|
| 67 | linear: [1,1000] |
|---|
| 68 | zoom: [1000,1010] |
|---|
| 69 | neg: [-100,100] |
|---|
| 70 | |
|---|
| 71 | *diff* is "relative", "absolute" or "none" |
|---|
| 72 | |
|---|
| 73 | *x_bits* is the precision with which the x values are specified. The |
|---|
| 74 | default 23 should reproduce the equivalent of a single precisio |
|---|
| 75 | """ |
|---|
| 76 | linear = not xrange.startswith("log") |
|---|
| 77 | if xrange == "zoom": |
|---|
| 78 | lin_min, lin_max, lin_steps = 1000, 1010, 2000 |
|---|
| 79 | elif xrange == "neg": |
|---|
| 80 | lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 |
|---|
| 81 | elif xrange == "linear": |
|---|
| 82 | lin_min, lin_max, lin_steps = 1, 1000, 2000 |
|---|
| 83 | elif xrange == "log": |
|---|
| 84 | log_min, log_max, log_steps = -3, 5, 400 |
|---|
| 85 | elif xrange == "logq": |
|---|
| 86 | log_min, log_max, log_steps = -4, 1, 400 |
|---|
| 87 | else: |
|---|
| 88 | raise ValueError("unknown range "+xrange) |
|---|
| 89 | with mp.workprec(500): |
|---|
| 90 | # Note: we make sure that we are comparing apples to apples... |
|---|
| 91 | # The x points are set using single precision so that we are |
|---|
| 92 | # examining the accuracy of the transformation from x to f(x) |
|---|
| 93 | # rather than x to f(nearest(x)) where nearest(x) is the nearest |
|---|
| 94 | # value to x in the given precision. |
|---|
| 95 | if linear: |
|---|
| 96 | lin_min = max(lin_min, self.limits[0]) |
|---|
| 97 | lin_max = min(lin_max, self.limits[1]) |
|---|
| 98 | qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single') |
|---|
| 99 | #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double') |
|---|
| 100 | qr = [mp.mpf(float(v)) for v in qrf] |
|---|
| 101 | #qr = mp.linspace(lin_min, lin_max, lin_steps) |
|---|
| 102 | else: |
|---|
| 103 | log_min = np.log10(max(10**log_min, self.limits[0])) |
|---|
| 104 | log_max = np.log10(min(10**log_max, self.limits[1])) |
|---|
| 105 | qrf = np.logspace(log_min, log_max, log_steps, dtype='single') |
|---|
| 106 | #qrf = np.logspace(log_min, log_max, log_steps, dtype='double') |
|---|
| 107 | qr = [mp.mpf(float(v)) for v in qrf] |
|---|
| 108 | #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] |
|---|
| 109 | |
|---|
| 110 | target = self.call_mpmath(qr, bits=500) |
|---|
| 111 | pylab.subplot(121) |
|---|
| 112 | self.compare(qr, 'single', target, linear, diff) |
|---|
| 113 | pylab.legend(loc='best') |
|---|
| 114 | pylab.subplot(122) |
|---|
| 115 | self.compare(qr, 'double', target, linear, diff) |
|---|
| 116 | pylab.legend(loc='best') |
|---|
| 117 | pylab.suptitle(self.name + " compared to 500-bit mpmath") |
|---|
| 118 | |
|---|
| 119 | def compare(self, x, precision, target, linear=False, diff="relative"): |
|---|
| 120 | r""" |
|---|
| 121 | Compare the different computation methods using the given precision. |
|---|
| 122 | """ |
|---|
| 123 | if precision == 'single': |
|---|
| 124 | #n=11; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff) |
|---|
| 125 | #n=23; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff) |
|---|
| 126 | pass |
|---|
| 127 | elif precision == 'double': |
|---|
| 128 | #n=53; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff) |
|---|
| 129 | #n=83; plotdiff(x, target, self.call_mpmath(x, n), 'mp %d bits'%n, diff=diff) |
|---|
| 130 | pass |
|---|
| 131 | plotdiff(x, target, self.call_numpy(x, precision), 'numpy '+precision, diff=diff) |
|---|
| 132 | plotdiff(x, target, self.call_ocl(x, precision, 0), 'OpenCL '+precision, diff=diff) |
|---|
| 133 | pylab.xlabel(self.xaxis) |
|---|
| 134 | if diff == "relative": |
|---|
| 135 | pylab.ylabel("relative error") |
|---|
| 136 | elif diff == "absolute": |
|---|
| 137 | pylab.ylabel("absolute error") |
|---|
| 138 | else: |
|---|
| 139 | pylab.ylabel(self.name) |
|---|
| 140 | pylab.semilogx(x, target, '-', label="true value") |
|---|
| 141 | if linear: |
|---|
| 142 | pylab.xscale('linear') |
|---|
| 143 | |
|---|
| 144 | def plotdiff(x, target, actual, label, diff): |
|---|
| 145 | """ |
|---|
| 146 | Plot the computed value. |
|---|
| 147 | |
|---|
| 148 | Use relative error if SHOW_DIFF, otherwise just plot the value directly. |
|---|
| 149 | """ |
|---|
| 150 | if diff == "relative": |
|---|
| 151 | err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') |
|---|
| 152 | #err = np.clip(err, 0, 1) |
|---|
| 153 | pylab.loglog(x, err, '-', label=label) |
|---|
| 154 | elif diff == "absolute": |
|---|
| 155 | err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd') |
|---|
| 156 | pylab.loglog(x, err, '-', label=label) |
|---|
| 157 | else: |
|---|
| 158 | limits = np.min(target), np.max(target) |
|---|
| 159 | pylab.semilogx(x, np.clip(actual, *limits), '-', label=label) |
|---|
| 160 | |
|---|
| 161 | def make_ocl(function, name, source=[]): |
|---|
| 162 | class Kernel(object): |
|---|
| 163 | pass |
|---|
| 164 | Kernel.__file__ = name+".py" |
|---|
| 165 | Kernel.name = name |
|---|
| 166 | Kernel.parameters = [] |
|---|
| 167 | Kernel.source = source |
|---|
| 168 | Kernel.Iq = function |
|---|
| 169 | model_info = modelinfo.make_model_info(Kernel) |
|---|
| 170 | return model_info |
|---|
| 171 | |
|---|
| 172 | |
|---|
| 173 | # =============== FUNCTION DEFINITIONS ================ |
|---|
| 174 | |
|---|
| 175 | FUNCTIONS = {} |
|---|
| 176 | def add_function(name, mp_function, np_function, ocl_function, |
|---|
| 177 | shortname=None, xaxis="x", limits=(-inf, inf)): |
|---|
| 178 | if shortname is None: |
|---|
| 179 | shortname = name.replace('(x)', '').replace(' ', '') |
|---|
| 180 | FUNCTIONS[shortname] = Comparator(name, mp_function, np_function, ocl_function, xaxis, limits) |
|---|
| 181 | |
|---|
| 182 | add_function( |
|---|
| 183 | name="J0(x)", |
|---|
| 184 | mp_function=mp.j0, |
|---|
| 185 | np_function=scipy.special.j0, |
|---|
| 186 | ocl_function=make_ocl("return sas_J0(q);", "sas_J0", ["lib/polevl.c", "lib/sas_J0.c"]), |
|---|
| 187 | ) |
|---|
| 188 | add_function( |
|---|
| 189 | name="J1(x)", |
|---|
| 190 | mp_function=mp.j1, |
|---|
| 191 | np_function=scipy.special.j1, |
|---|
| 192 | ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 193 | ) |
|---|
| 194 | add_function( |
|---|
| 195 | name="JN(-3, x)", |
|---|
| 196 | mp_function=lambda x: mp.besselj(-3, x), |
|---|
| 197 | np_function=lambda x: scipy.special.jn(-3, x), |
|---|
| 198 | ocl_function=make_ocl("return sas_JN(-3, q);", "sas_JN", |
|---|
| 199 | ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]), |
|---|
| 200 | shortname="J-3", |
|---|
| 201 | ) |
|---|
| 202 | add_function( |
|---|
| 203 | name="JN(3, x)", |
|---|
| 204 | mp_function=lambda x: mp.besselj(3, x), |
|---|
| 205 | np_function=lambda x: scipy.special.jn(3, x), |
|---|
| 206 | ocl_function=make_ocl("return sas_JN(3, q);", "sas_JN", |
|---|
| 207 | ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]), |
|---|
| 208 | shortname="J3", |
|---|
| 209 | ) |
|---|
| 210 | add_function( |
|---|
| 211 | name="JN(2, x)", |
|---|
| 212 | mp_function=lambda x: mp.besselj(2, x), |
|---|
| 213 | np_function=lambda x: scipy.special.jn(2, x), |
|---|
| 214 | ocl_function=make_ocl("return sas_JN(2, q);", "sas_JN", |
|---|
| 215 | ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c"]), |
|---|
| 216 | shortname="J2", |
|---|
| 217 | ) |
|---|
| 218 | add_function( |
|---|
| 219 | name="2 J1(x)/x", |
|---|
| 220 | mp_function=lambda x: 2*mp.j1(x)/x, |
|---|
| 221 | np_function=lambda x: 2*scipy.special.j1(x)/x, |
|---|
| 222 | ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 223 | ) |
|---|
| 224 | add_function( |
|---|
| 225 | name="J1(x)", |
|---|
| 226 | mp_function=mp.j1, |
|---|
| 227 | np_function=scipy.special.j1, |
|---|
| 228 | ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 229 | ) |
|---|
| 230 | add_function( |
|---|
| 231 | name="Si(x)", |
|---|
| 232 | mp_function=mp.si, |
|---|
| 233 | np_function=lambda x: scipy.special.sici(x)[0], |
|---|
| 234 | ocl_function=make_ocl("return sas_Si(q);", "sas_Si", ["lib/sas_Si.c"]), |
|---|
| 235 | ) |
|---|
| 236 | #import fnlib |
|---|
| 237 | #add_function( |
|---|
| 238 | # name="fnlibJ1", |
|---|
| 239 | # mp_function=mp.j1, |
|---|
| 240 | # np_function=fnlib.J1, |
|---|
| 241 | # ocl_function=make_ocl("return sas_J1(q);", "sas_J1", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 242 | #) |
|---|
| 243 | add_function( |
|---|
| 244 | name="sin(x)", |
|---|
| 245 | mp_function=mp.sin, |
|---|
| 246 | np_function=np.sin, |
|---|
| 247 | #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return sn;", "sas_sin"), |
|---|
| 248 | ocl_function=make_ocl("return sin(q);", "sas_sin"), |
|---|
| 249 | ) |
|---|
| 250 | add_function( |
|---|
| 251 | name="sin(x)/x", |
|---|
| 252 | mp_function=lambda x: mp.sin(x)/x if x != 0 else 1, |
|---|
| 253 | ## scipy sinc function is inaccurate and has an implied pi*x term |
|---|
| 254 | #np_function=lambda x: scipy.special.sinc(x/pi), |
|---|
| 255 | ## numpy sin(x)/x needs to check for x=0 |
|---|
| 256 | np_function=lambda x: np.sin(x)/x, |
|---|
| 257 | ocl_function=make_ocl("return sas_sinx_x(q);", "sas_sinc"), |
|---|
| 258 | ) |
|---|
| 259 | add_function( |
|---|
| 260 | name="cos(x)", |
|---|
| 261 | mp_function=mp.cos, |
|---|
| 262 | np_function=np.cos, |
|---|
| 263 | #ocl_function=make_ocl("double sn, cn; SINCOS(q,sn,cn); return cn;", "sas_cos"), |
|---|
| 264 | ocl_function=make_ocl("return cos(q);", "sas_cos"), |
|---|
| 265 | ) |
|---|
| 266 | add_function( |
|---|
| 267 | name="gamma(x)", |
|---|
| 268 | mp_function=mp.gamma, |
|---|
| 269 | np_function=scipy.special.gamma, |
|---|
| 270 | ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]), |
|---|
| 271 | limits=(-3.1, 10), |
|---|
| 272 | ) |
|---|
| 273 | add_function( |
|---|
| 274 | name="erf(x)", |
|---|
| 275 | mp_function=mp.erf, |
|---|
| 276 | np_function=scipy.special.erf, |
|---|
| 277 | ocl_function=make_ocl("return sas_erf(q);", "sas_erf", ["lib/polevl.c", "lib/sas_erf.c"]), |
|---|
| 278 | limits=(-5., 5.), |
|---|
| 279 | ) |
|---|
| 280 | add_function( |
|---|
| 281 | name="erfc(x)", |
|---|
| 282 | mp_function=mp.erfc, |
|---|
| 283 | np_function=scipy.special.erfc, |
|---|
| 284 | ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]), |
|---|
| 285 | limits=(-5., 5.), |
|---|
| 286 | ) |
|---|
| 287 | add_function( |
|---|
| 288 | name="arctan(x)", |
|---|
| 289 | mp_function=mp.atan, |
|---|
| 290 | np_function=np.arctan, |
|---|
| 291 | ocl_function=make_ocl("return atan(q);", "sas_arctan"), |
|---|
| 292 | ) |
|---|
| 293 | add_function( |
|---|
| 294 | name="3 j1(x)/x", |
|---|
| 295 | mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x), |
|---|
| 296 | # Note: no taylor expansion near 0 |
|---|
| 297 | np_function=lambda x: 3*(np.sin(x)/x - np.cos(x))/(x*x), |
|---|
| 298 | ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]), |
|---|
| 299 | ) |
|---|
| 300 | add_function( |
|---|
| 301 | name="(1-cos(x))/x^2", |
|---|
| 302 | mp_function=lambda x: (1 - mp.cos(x))/(x*x), |
|---|
| 303 | np_function=lambda x: (1 - np.cos(x))/(x*x), |
|---|
| 304 | ocl_function=make_ocl("return (1-cos(q))/q/q;", "sas_1mcosx_x2"), |
|---|
| 305 | ) |
|---|
| 306 | add_function( |
|---|
| 307 | name="(1-sin(x)/x)/x", |
|---|
| 308 | mp_function=lambda x: 1/x - mp.sin(x)/(x*x), |
|---|
| 309 | np_function=lambda x: 1/x - np.sin(x)/(x*x), |
|---|
| 310 | ocl_function=make_ocl("return (1-sas_sinx_x(q))/q;", "sas_1msinx_x_x"), |
|---|
| 311 | ) |
|---|
| 312 | add_function( |
|---|
| 313 | name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x", |
|---|
| 314 | mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x, |
|---|
| 315 | np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x, |
|---|
| 316 | ocl_function=make_ocl("return (0.5-sin(q)/q + (1-cos(q))/q/q)/q;", "sas_T2"), |
|---|
| 317 | ) |
|---|
| 318 | add_function( |
|---|
| 319 | name="fmod_2pi", |
|---|
| 320 | mp_function=lambda x: mp.fmod(x, 2*mp.pi), |
|---|
| 321 | np_function=lambda x: np.fmod(x, 2*np.pi), |
|---|
| 322 | ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"), |
|---|
| 323 | ) |
|---|
| 324 | |
|---|
| 325 | RADIUS=3000 |
|---|
| 326 | LENGTH=30 |
|---|
| 327 | THETA=45 |
|---|
| 328 | def mp_cyl(x): |
|---|
| 329 | f = mp.mpf |
|---|
| 330 | theta = f(THETA)*mp.pi/f(180) |
|---|
| 331 | qr = x * f(RADIUS)*mp.sin(theta) |
|---|
| 332 | qh = x * f(LENGTH)/f(2)*mp.cos(theta) |
|---|
| 333 | be = f(2)*mp.j1(qr)/qr |
|---|
| 334 | si = mp.sin(qh)/qh |
|---|
| 335 | background = f(0) |
|---|
| 336 | #background = f(1)/f(1000) |
|---|
| 337 | volume = mp.pi*f(RADIUS)**f(2)*f(LENGTH) |
|---|
| 338 | contrast = f(5) |
|---|
| 339 | units = f(1)/f(10000) |
|---|
| 340 | #return be |
|---|
| 341 | #return si |
|---|
| 342 | return units*(volume*contrast*be*si)**f(2)/volume + background |
|---|
| 343 | def np_cyl(x): |
|---|
| 344 | f = np.float64 if x.dtype == np.float64 else np.float32 |
|---|
| 345 | theta = f(THETA)*f(np.pi)/f(180) |
|---|
| 346 | qr = x * f(RADIUS)*np.sin(theta) |
|---|
| 347 | qh = x * f(LENGTH)/f(2)*np.cos(theta) |
|---|
| 348 | be = f(2)*scipy.special.j1(qr)/qr |
|---|
| 349 | si = np.sin(qh)/qh |
|---|
| 350 | background = f(0) |
|---|
| 351 | #background = f(1)/f(1000) |
|---|
| 352 | volume = f(np.pi)*f(RADIUS)**2*f(LENGTH) |
|---|
| 353 | contrast = f(5) |
|---|
| 354 | units = f(1)/f(10000) |
|---|
| 355 | #return be |
|---|
| 356 | #return si |
|---|
| 357 | return units*(volume*contrast*be*si)**f(2)/volume + background |
|---|
| 358 | ocl_cyl = """\ |
|---|
| 359 | double THETA = %(THETA).15e*M_PI_180; |
|---|
| 360 | double qr = q*%(RADIUS).15e*sin(THETA); |
|---|
| 361 | double qh = q*0.5*%(LENGTH).15e*cos(THETA); |
|---|
| 362 | double be = sas_2J1x_x(qr); |
|---|
| 363 | double si = sas_sinx_x(qh); |
|---|
| 364 | double background = 0; |
|---|
| 365 | //double background = 0.001; |
|---|
| 366 | double volume = M_PI*square(%(RADIUS).15e)*%(LENGTH).15e; |
|---|
| 367 | double contrast = 5.0; |
|---|
| 368 | double units = 1e-4; |
|---|
| 369 | //return be; |
|---|
| 370 | //return si; |
|---|
| 371 | return units*square(volume*contrast*be*si)/volume + background; |
|---|
| 372 | """%{"LENGTH":LENGTH, "RADIUS": RADIUS, "THETA": THETA} |
|---|
| 373 | add_function( |
|---|
| 374 | name="cylinder(r=%g, l=%g, theta=%g)"%(RADIUS, LENGTH, THETA), |
|---|
| 375 | mp_function=mp_cyl, |
|---|
| 376 | np_function=np_cyl, |
|---|
| 377 | ocl_function=make_ocl(ocl_cyl, "ocl_cyl", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 378 | shortname="cylinder", |
|---|
| 379 | xaxis="$q/A^{-1}$", |
|---|
| 380 | ) |
|---|
| 381 | |
|---|
| 382 | lanczos_gamma = """\ |
|---|
| 383 | const double coeff[] = { |
|---|
| 384 | 76.18009172947146, -86.50532032941677, |
|---|
| 385 | 24.01409824083091, -1.231739572450155, |
|---|
| 386 | 0.1208650973866179e-2,-0.5395239384953e-5 |
|---|
| 387 | }; |
|---|
| 388 | const double x = q; |
|---|
| 389 | double tmp = x + 5.5; |
|---|
| 390 | tmp -= (x + 0.5)*log(tmp); |
|---|
| 391 | double ser = 1.000000000190015; |
|---|
| 392 | for (int k=0; k < 6; k++) ser += coeff[k]/(x + k+1); |
|---|
| 393 | return -tmp + log(2.5066282746310005*ser/x); |
|---|
| 394 | """ |
|---|
| 395 | add_function( |
|---|
| 396 | name="log gamma(x)", |
|---|
| 397 | mp_function=mp.loggamma, |
|---|
| 398 | np_function=scipy.special.gammaln, |
|---|
| 399 | ocl_function=make_ocl(lanczos_gamma, "lgamma"), |
|---|
| 400 | ) |
|---|
| 401 | |
|---|
| 402 | # Alternate versions of 3 j1(x)/x, for posterity |
|---|
| 403 | def taylor_3j1x_x(x): |
|---|
| 404 | """ |
|---|
| 405 | Calculation using taylor series. |
|---|
| 406 | """ |
|---|
| 407 | # Generate coefficients using the precision of the target value. |
|---|
| 408 | n = 5 |
|---|
| 409 | cinv = [3991680, -45360, 840, -30, 3] |
|---|
| 410 | three = x.dtype.type(3) |
|---|
| 411 | p = three/np.array(cinv, x.dtype) |
|---|
| 412 | return np.polyval(p[-n:], x*x) |
|---|
| 413 | add_function( |
|---|
| 414 | name="3 j1(x)/x: taylor", |
|---|
| 415 | mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x), |
|---|
| 416 | np_function=taylor_3j1x_x, |
|---|
| 417 | ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]), |
|---|
| 418 | ) |
|---|
| 419 | def trig_3j1x_x(x): |
|---|
| 420 | r""" |
|---|
| 421 | Direct calculation using linear combination of sin/cos. |
|---|
| 422 | |
|---|
| 423 | Use the following trig identity: |
|---|
| 424 | |
|---|
| 425 | .. math:: |
|---|
| 426 | |
|---|
| 427 | a \sin(x) + b \cos(x) = c \sin(x + \phi) |
|---|
| 428 | |
|---|
| 429 | where $c = \surd(a^2+b^2)$ and $\phi = \tan^{-1}(b/a) to calculate the |
|---|
| 430 | numerator $\sin(x) - x\cos(x)$. |
|---|
| 431 | """ |
|---|
| 432 | one = x.dtype.type(1) |
|---|
| 433 | three = x.dtype.type(3) |
|---|
| 434 | c = np.sqrt(one + x*x) |
|---|
| 435 | phi = np.arctan2(-x, one) |
|---|
| 436 | return three*(c*np.sin(x+phi))/(x*x*x) |
|---|
| 437 | add_function( |
|---|
| 438 | name="3 j1(x)/x: trig", |
|---|
| 439 | mp_function=lambda x: 3*(mp.sin(x)/x - mp.cos(x))/(x*x), |
|---|
| 440 | np_function=trig_3j1x_x, |
|---|
| 441 | ocl_function=make_ocl("return sas_3j1x_x(q);", "sas_j1c", ["lib/sas_3j1x_x.c"]), |
|---|
| 442 | ) |
|---|
| 443 | def np_2J1x_x(x): |
|---|
| 444 | """ |
|---|
| 445 | numpy implementation of 2J1(x)/x using single precision algorithm |
|---|
| 446 | """ |
|---|
| 447 | # pylint: disable=bad-continuation |
|---|
| 448 | f = x.dtype.type |
|---|
| 449 | ax = abs(x) |
|---|
| 450 | if ax < f(8.0): |
|---|
| 451 | y = x*x |
|---|
| 452 | ans1 = f(2)*(f(72362614232.0) |
|---|
| 453 | + y*(f(-7895059235.0) |
|---|
| 454 | + y*(f(242396853.1) |
|---|
| 455 | + y*(f(-2972611.439) |
|---|
| 456 | + y*(f(15704.48260) |
|---|
| 457 | + y*(f(-30.16036606))))))) |
|---|
| 458 | ans2 = (f(144725228442.0) |
|---|
| 459 | + y*(f(2300535178.0) |
|---|
| 460 | + y*(f(18583304.74) |
|---|
| 461 | + y*(f(99447.43394) |
|---|
| 462 | + y*(f(376.9991397) |
|---|
| 463 | + y))))) |
|---|
| 464 | return ans1/ans2 |
|---|
| 465 | else: |
|---|
| 466 | y = f(64.0)/(ax*ax) |
|---|
| 467 | xx = ax - f(2.356194491) |
|---|
| 468 | ans1 = (f(1.0) |
|---|
| 469 | + y*(f(0.183105e-2) |
|---|
| 470 | + y*(f(-0.3516396496e-4) |
|---|
| 471 | + y*(f(0.2457520174e-5) |
|---|
| 472 | + y*f(-0.240337019e-6))))) |
|---|
| 473 | ans2 = (f(0.04687499995) |
|---|
| 474 | + y*(f(-0.2002690873e-3) |
|---|
| 475 | + y*(f(0.8449199096e-5) |
|---|
| 476 | + y*(f(-0.88228987e-6) |
|---|
| 477 | + y*f(0.105787412e-6))))) |
|---|
| 478 | sn, cn = np.sin(xx), np.cos(xx) |
|---|
| 479 | ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x |
|---|
| 480 | return -ans if (x < f(0.0)) else ans |
|---|
| 481 | add_function( |
|---|
| 482 | name="2 J1(x)/x:alt", |
|---|
| 483 | mp_function=lambda x: 2*mp.j1(x)/x, |
|---|
| 484 | np_function=lambda x: np.asarray([np_2J1x_x(v) for v in x], x.dtype), |
|---|
| 485 | ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]), |
|---|
| 486 | ) |
|---|
| 487 | |
|---|
| 488 | ALL_FUNCTIONS = set(FUNCTIONS.keys()) |
|---|
| 489 | ALL_FUNCTIONS.discard("loggamma") # OCL version not ready yet |
|---|
| 490 | ALL_FUNCTIONS.discard("3j1/x:taylor") |
|---|
| 491 | ALL_FUNCTIONS.discard("3j1/x:trig") |
|---|
| 492 | ALL_FUNCTIONS.discard("2J1/x:alt") |
|---|
| 493 | |
|---|
| 494 | # =============== MAIN PROGRAM ================ |
|---|
| 495 | |
|---|
| 496 | def usage(): |
|---|
| 497 | names = ", ".join(sorted(ALL_FUNCTIONS)) |
|---|
| 498 | print("""\ |
|---|
| 499 | usage: precision.py [-f/a/r] [-x<range>] name... |
|---|
| 500 | where |
|---|
| 501 | -f indicates that the function value should be plotted, |
|---|
| 502 | -a indicates that the absolute error should be plotted, |
|---|
| 503 | -r indicates that the relative error should be plotted (default), |
|---|
| 504 | -x<range> indicates the steps in x, where <range> is one of the following |
|---|
| 505 | log indicates log stepping in [10^-3, 10^5] (default) |
|---|
| 506 | logq indicates log stepping in [10^-4, 10^1] |
|---|
| 507 | linear indicates linear stepping in [1, 1000] |
|---|
| 508 | zoom indicates linear stepping in [1000, 1010] |
|---|
| 509 | neg indicates linear stepping in [-100.1, 100.1] |
|---|
| 510 | and name is "all [first]" or one of: |
|---|
| 511 | """+names) |
|---|
| 512 | sys.exit(1) |
|---|
| 513 | |
|---|
| 514 | def main(): |
|---|
| 515 | import sys |
|---|
| 516 | diff = "relative" |
|---|
| 517 | xrange = "log" |
|---|
| 518 | options = [v for v in sys.argv[1:] if v.startswith('-')] |
|---|
| 519 | for opt in options: |
|---|
| 520 | if opt == '-f': |
|---|
| 521 | diff = "none" |
|---|
| 522 | elif opt == '-r': |
|---|
| 523 | diff = "relative" |
|---|
| 524 | elif opt == '-a': |
|---|
| 525 | diff = "absolute" |
|---|
| 526 | elif opt.startswith('-x'): |
|---|
| 527 | xrange = opt[2:] |
|---|
| 528 | else: |
|---|
| 529 | usage() |
|---|
| 530 | |
|---|
| 531 | names = [v for v in sys.argv[1:] if not v.startswith('-')] |
|---|
| 532 | if not names: |
|---|
| 533 | usage() |
|---|
| 534 | |
|---|
| 535 | if names[0] == "all": |
|---|
| 536 | cutoff = names[1] if len(names) > 1 else "" |
|---|
| 537 | names = list(sorted(ALL_FUNCTIONS)) |
|---|
| 538 | names = [k for k in names if k >= cutoff] |
|---|
| 539 | if any(k not in FUNCTIONS for k in names): |
|---|
| 540 | usage() |
|---|
| 541 | multiple = len(names) > 1 |
|---|
| 542 | pylab.interactive(multiple) |
|---|
| 543 | for k in names: |
|---|
| 544 | pylab.clf() |
|---|
| 545 | comparator = FUNCTIONS[k] |
|---|
| 546 | comparator.run(xrange=xrange, diff=diff) |
|---|
| 547 | if multiple: |
|---|
| 548 | raw_input() |
|---|
| 549 | if not multiple: |
|---|
| 550 | pylab.show() |
|---|
| 551 | |
|---|
| 552 | if __name__ == "__main__": |
|---|
| 553 | main() |
|---|
