From 5ac6f8334ac97e19b0dd59adf5a300293996117a Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 13 Mar 2024 12:11:21 -0400 Subject: [PATCH] add precision checks for the two yukawa expressions [tweak cutoff values] --- explore/precision.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/explore/precision.py b/explore/precision.py index 97e047f6..735328af 100755 --- a/explore/precision.py +++ b/explore/precision.py @@ -392,9 +392,9 @@ def add_function(name, mp_function, np_function, ocl_function, name="(1-cos(x))/x^2", mp_function=lambda x: (1 - mp.cos(x))/(x*x), #np_function=lambda x: (1 - np.cos(x))/(x*x), - np_function=lambda x: ( - (x>0.5)*(1 - np.cos(x))/(x*x) - +(x<=0.5)*(1.0/2.0 + x*x*(1./-24. + x*x*(1./720. + x*x*(1./-40320. + x*x*(1./3628800. + x*x*(1./-479001600. + x*x*1./87178291200.)))))) + np_function=lambda x, cutoff=0.5: ( + (x>cutoff)*(1 - np.cos(x))/(x*x) + +(x<=cutoff)*(1.0/2.0 + x*x*(1./-24. + x*x*(1./720. + x*x*(1./-40320. + x*x*(1./3628800. + x*x*(1./-479001600. + x*x*1./87178291200.)))))) ), ocl_function=make_ocl("return (1-cos(q))/q/q;", "sas_1mcosx_x2"), ) @@ -402,27 +402,28 @@ def add_function(name, mp_function, np_function, ocl_function, # Note: this is just j1(x)/x name="TY_pc_t1", mp_function=lambda x: (mp.sin(x)/x - mp.cos(x))/(x*x), - np_function=lambda x: ( - (x>0.5)*(np.sin(x)/x - np.cos(x))/(x*x) - +(x<=0.5)*(1.0/3.0 + x*x*(1./-30. + x*x*(1./840. + x*x*(1./-45360. + x*x*(1./3991680. + x*x*(1./-518918400. + x*x*(1./93405312000.))))))) + np_function=lambda x, cutoff=0.5: ( + (x>cutoff)*(np.sin(x)/x - np.cos(x))/(x*x) + +(x<=cutoff)*(1.0/3.0 + x*x*(1./-30. + x*x*(1./840. + x*x*(1./-45360. + x*x*(1./3991680. + x*x*(1./-518918400. + x*x*(1./93405312000.))))))) ), - ocl_function=make_ocl("return sas_3j1x_x(q)/3.;", "sas_j1x_x", ["lib/sas_3j1x_x.c"]), + #ocl_function=make_ocl("return sas_3j1x_x(q)/3.;", "sas_j1x_x", ["lib/sas_3j1x_x.c"]), + ocl_function=make_ocl("return (sin(q)/q - cos(q))/(q*q);", "j1x_x"), ) add_function( name="TY_pc_t2", mp_function=lambda x: (2.*x*mp.sin(x) - (x*x - 2.)*mp.cos(x) - 2.)/x**4, - np_function=lambda x: ( - (x>1.0)*(2.*x*np.sin(x) - (x*x - 2.)*np.cos(x) - 2.)/x**4 - +(x<=1.0)*(1./4. + x*x*(1./-36 + x*x*(1./960. + x*x*(1./-50400. + x*x*(1./4354560. + x*x*(1./-558835200. + x*x*(1./99632332800.))))))) + np_function=lambda x, cutoff=0.5: ( + (x>cutoff)*(2.*x*np.sin(x) - (x*x - 2.)*np.cos(x) - 2.)/x**4 + +(x<=cutoff)*(1./4. + x*x*(1./-36 + x*x*(1./960. + x*x*(1./-50400. + x*x*(1./4354560. + x*x*(1./-558835200. + x*x*(1./99632332800.))))))) ), ocl_function=make_ocl("return (2.*q*sin(q) - (q*q - 2.)*cos(q) - 2.) / pown(q,4);","TY_pc_t2"), ) add_function( name="TY_pc_t3", mp_function=lambda x: ( 4*x*(x*x - 6)*mp.sin(x) - (x**4 - 12*x*x + 24)*mp.cos(x) + 24 ) / x**6, - np_function=lambda x: ( - (x>1.0)*(4.*x*(x*x - 6.)*np.sin(x) - ((x*x - 12.)*x*x + 24.)*np.cos(x) + 24.) / x**6 - +(x<=1.0)*(1/6 + x*x*(1./-48. + x*x*(1./1200. + x*x*(1./-60480. + x*x*(1./5080320. + x*x*(1./-638668800. + x*x*(1./112086374400.))))))) + np_function=lambda x, cutoff=0.9: ( + (x>cutoff)*(4.*x*(x*x - 6.)*np.sin(x) - ((x*x - 12.)*x*x + 24.)*np.cos(x) + 24.) / x**6 + +(x<=cutoff)*(1/6 + x*x*(1./-48. + x*x*(1./1200. + x*x*(1./-60480. + x*x*(1./5080320. + x*x*(1./-638668800. + x*x*(1./112086374400.))))))) ), ocl_function=make_ocl("return (4*q*(q*q - 6)*sin(q) - ((q*q - 12)*q*q + 24)*cos(q) + 24 ) / pown(q,6);","ty_pc_t2"), )