Changeset 31eea1f in sasmodels
 Timestamp:
 Oct 22, 2017 5:00:50 PM (21 months ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 f728001
 Parents:
 2bccb5a
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

explore/asymint.py
r4f611f1 r31eea1f 9 9 10 10 import numpy as np 11 import mpmath as mp 11 12 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 13 from numpy.polynomial.legendre import leggauss 12 14 from scipy.integrate import dblquad, simps, romb, romberg 13 15 import pylab … … 19 21 from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x 20 22 21 SLD = 3.0 22 SLD_SOLVENT = 6.3 23 def mp_3j1x_x(x): 24 return 3*(mp.sin(x)/x  mp.cos(x))/(x*x) 25 def mp_2J1x_x(x): 26 return 2*mp.j1(x)/x 27 def mp_sinx_x(x): 28 return mp.sin(x)/x 29 30 SLD = 3 31 SLD_SOLVENT = 6 23 32 CONTRAST = SLD  SLD_SOLVENT 24 33 34 def make_parallelepiped(a, b, c): 35 def Fq(qa, qb, qc): 36 siA = sas_sinx_x(0.5*a*qa) 37 siB = sas_sinx_x(0.5*b*qb) 38 siC = sas_sinx_x(0.5*c*qc) 39 return siA * siB * siC 40 volume = a*b*c 41 norm = volume*CONTRAST**2/10**4 42 return norm, Fq 43 44 def make_parallelepiped_mp(a, b, c): 45 a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) 46 def Fq(qa, qb, qc): 47 siA = mp_sinx_x(a*qa/2) 48 siB = mp_sinx_x(b*qb/2) 49 siC = mp_sinx_x(c*qc/2) 50 return siA * siB * siC 51 volume = a*b*c 52 norm = (volume*CONTRAST**2)/10000 # mpf since volume=a*b*c is mpf 53 return norm, Fq 54 55 def make_triellip(a, b, c): 56 def Fq(qa, qb, qc): 57 qr = sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) 58 return sas_3j1x_x(qr) 59 volume = 4*pi*a*b*c/3 60 norm = volume*CONTRAST**2/10**4 61 return norm, Fq 62 63 def make_triellip_mp(a, b, c): 64 a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) 65 def Fq(qa, qb, qc): 66 qr = mp.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) 67 return mp_3j1x_x(qr) 68 volume = (4*mp.pi*a*b*c)/3 69 norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf 70 return norm, Fq 71 25 72 def make_cylinder(radius, length): 26 def cylinder(qa, qb, qc):73 def Fq(qa, qb, qc): 27 74 qab = sqrt(qa**2 + qb**2) 28 75 return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length) 29 76 volume = pi*radius**2*length 30 norm = 1e4*volume*CONTRAST**2 31 return norm, cylinder 77 norm = volume*CONTRAST**2/10**4 78 return norm, Fq 79 80 def make_cylinder_mp(radius, length): 81 radius, length = mp.mpf(radius), mp.mpf(length) 82 def Fq(qa, qb, qc): 83 qab = mp.sqrt(qa**2 + qb**2) 84 return mp_2J1x_x(qab*radius) * mp_sinx_x((qc*length)/2) 85 volume = mp.pi*radius**2*length 86 norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf 87 return norm, Fq 32 88 33 89 def make_sphere(radius): 34 def sphere(qa, qb, qc): 35 qab = sqrt(qa**2 + qb**2) 36 q = sqrt(qab**2 + qc**2) 90 def Fq(qa, qb, qc): 91 q = sqrt(qa**2 + qb**2 + qc**2) 37 92 return sas_3j1x_x(q*radius) 38 93 volume = 4*pi*radius**3/3 39 norm = 1e4*volume*CONTRAST**2 40 return norm, sphere 41 42 #NORM, KERNEL = make_cylinder(radius=10., length=100000.) 43 NORM, KERNEL = make_cylinder(radius=10., length=30.) 44 #NORM, KERNEL = make_sphere(radius=50.) 94 norm = volume*CONTRAST**2/10**4 95 return norm, Fq 96 97 def make_sphere_mp(radius): 98 radius = mp.mpf(radius) 99 def Fq(qa, qb, qc): 100 q = mp.sqrt(qa**2 + qb**2 + qc**2) 101 return mp_3j1x_x(q*radius) 102 volume = (4*mp.pi*radius**3)/3 103 norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf 104 return norm, Fq 105 106 shape = 'parallelepiped' 107 #shape = 'triellip' 108 #shape = 'sphere' 109 #shape = 'cylinder' 110 if shape == 'cylinder': 111 #RADIUS, LENGTH = 10, 100000 112 RADIUS, LENGTH = 10, 300 # integer for the sake of mpf 113 NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 114 NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH) 115 elif shape == 'triellip': 116 #A, B, C = 4450, 14000, 47 117 A, B, C = 445, 140, 47 # integer for the sake of mpf 118 NORM, KERNEL = make_triellip(A, B, C) 119 NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C) 120 elif shape == 'parallelepiped': 121 #A, B, C = 4450, 14000, 47 122 A, B, C = 445, 140, 47 # integer for the sake of mpf 123 NORM, KERNEL = make_parallelepiped(A, B, C) 124 NORM_MP, KERNEL_MP = make_parallelepiped_mp(A, B, C) 125 elif shape == 'sphere': 126 RADIUS = 50 # integer for the sake of mpf 127 NORM, KERNEL = make_sphere(radius=RADIUS) 128 NORM_MP, KERNEL_MP = make_sphere_mp(radius=RADIUS) 129 else: 130 raise ValueError("Unknown shape %r"%shape) 45 131 46 132 THETA_LOW, THETA_HIGH = 0, pi … … 48 134 SCALE = 1 49 135 136 # mathematica code for triaxial_ellipsoid (untested) 137 _ = """ 138 R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2] 139 Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r) 140 V[a_, b_, c_] := 4/3 pi a b c 141 Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent  sld)^2 142 F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]] 143 I[q_, sld_, solvent_, a_, b_, c_] := Norm[sld, solvent, a, b, c]/(4 pi) Integrate[F[q, theta, phi, a, b, c]^2 Sin[theta], {phi, 0, 2 pi}, {theta, 0, pi}] 144 I[6/10^3, 63/10, 3, 445, 140, 47] 145 """ 146 147 148 def mp_quad(q, shape): 149 evals = [0] 150 def integrand(theta, phi): 151 evals[0] += 1 152 qab = q*mp.sin(theta) 153 qa = qab*mp.cos(phi) 154 qb = qab*mp.sin(phi) 155 qc = q*mp.cos(theta) 156 Zq = KERNEL_MP(qa, qb, qc)**2 157 return Zq*mp.sin(theta) 158 ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi)) 159 Iq = NORM_MP*ans/(4*mp.pi) 160 return evals[0], Iq 50 161 51 162 def kernel(q, theta, phi): … … 59 170 return NORM*KERNEL(qa, qb, qc)**2 60 171 61 62 172 def scipy_dblquad(q): 63 173 """ … … 66 176 """ 67 177 evals = [0] 68 def integrand( theta, phi):178 def integrand(phi, theta): 69 179 evals[0] += 1 70 180 Zq = kernel(q, theta=theta, phi=phi) 71 181 return Zq*sin(theta) 72 ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]*SCALE/(4*pi) 73 print("dblquad evals =", evals[0]) 74 return ans 182 ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0] 183 return evals[0], ans*SCALE/(4*pi) 75 184 76 185 … … 80 189 complete in a reasonable time. No idea if it is accurate. 81 190 """ 191 evals = [0] 82 192 def inner(phi, theta): 193 evals[0] += 1 83 194 return kernel(q, theta=theta, phi=phi) 84 195 def outer(theta): 85 return romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))*sin(theta) 86 return romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE/(4*pi) 196 Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,)) 197 return Zq*sin(theta) 198 ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100) 199 return evals[0], ans*SCALE/(4*pi) 87 200 88 201 … … 96 209 return kernel(q, theta=theta, phi=phi) 97 210 theta = np.linspace(THETA_LOW, THETA_HIGH, n) 98 f_phi = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) 99 for t in theta] 100 ans = simps(sin(theta)*np.array(f_phi), dx=theta[1]theta[0]) 101 print("semi romberg evals =", evals[0]) 102 return ans*SCALE/(4*pi) 211 Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta] 212 ans = simps(np.array(Zq)*sin(theta), dx=theta[1]theta[0]) 213 return evals[0], ans*SCALE/(4*pi) 103 214 104 215 def gauss_quad(q, n=150): … … 110 221 elif n == 76: 111 222 z, w = Gauss76Z, Gauss76Wt 223 elif n == 150: 224 z, w = Gauss150Z, Gauss150Wt 112 225 else: 113 z, w = Gauss150Z, Gauss150Wt226 z, w = leggauss(n) 114 227 theta = (THETA_HIGHTHETA_LOW)*(z + 1)/2 + THETA_LOW 115 228 phi = (PHI_HIGHPHI_LOW)*(z + 1)/2 + PHI_LOW 116 229 Atheta, Aphi = np.meshgrid(theta, phi) 117 230 Aw = w[None, :] * w[:, None] 118 sin_theta = np.fmax(abs(sin(Atheta)), 1e6)231 sin_theta = abs(sin(Atheta)) 119 232 Zq = kernel(q=q, theta=Atheta, phi=Aphi) 120 print("gauss %d evals ="%n, n**2) 121 return np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) 233 # change from [1,1] x [1,1] range to [0, pi] x [0, 2 pi] range 234 dxdy_stretch = (THETA_HIGHTHETA_LOW)/2 * (PHI_HIGHPHI_LOW)/2 235 Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch 236 return n**2, Iq 122 237 123 238 … … 134 249 Zq *= abs(sin(Atheta)) 135 250 dx, dy = theta[1]theta[0], phi[1]phi[0] 136 print("rect", n, np.sum(Zq)*dx*dy*SCALE/(4*pi)) 137 print("trapz", n, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 138 print("simpson", n, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 139 print("romb", n, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 140 print("gridded %d evals ="%n, n**2) 251 print("rect%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi)) 252 print("trapz%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 253 print("simpson%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 254 print("romb%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 141 255 142 256 def plot(q, n=300): … … 161 275 162 276 if __name__ == "__main__": 163 Q = 0.8 164 print("gauss", 20, gauss_quad(Q, n=20)) 165 print("gauss", 76, gauss_quad(Q, n=76)) 166 print("gauss", 150, gauss_quad(Q, n=150)) 167 #print("dblquad", scipy_dblquad(Q)) 168 #print("semi romberg", semi_romberg(Q)) 277 Qstr = '0.005' 278 #Qstr = '0.8' 279 #Qstr = '0.0003' 280 Q = float(Qstr) 281 if shape == 'sphere': 282 print("exact", NORM*sas_3j1x_x(Q*RADIUS)**2) 283 print("gauss20", *gauss_quad(Q, n=20)) 284 print("gauss76", *gauss_quad(Q, n=76)) 285 print("gauss150", *gauss_quad(Q, n=150)) 286 print("gauss500", *gauss_quad(Q, n=500)) 287 print("dblquad", *scipy_dblquad(Q)) 288 print("semiromberg100", *semi_romberg(Q, n=100)) 289 print("romberg", *scipy_romberg_2d(Q)) 169 290 #gridded_integrals(Q, n=2**8+1) 170 #gridded_integrals(Q, n=2**10+1)291 gridded_integrals(Q, n=2**10+1) 171 292 #gridded_integrals(Q, n=2**13+1) 172 #print("romberg", scipy_romberg(Q)) 173 plot(Q, n=400) 293 #gridded_integrals(Q, n=2**15+1) 294 with mp.workprec(100): 295 print("mpmath", *mp_quad(mp.mpf(Qstr), shape)) 296 #plot(Q, n=200) 
sasmodels/compare.py
re3571cb r31eea1f 1196 1196 elif arg.startswith('accuracy='): opts['accuracy'] = arg[10:] 1197 1197 elif arg.startswith('cutoff='): opts['cutoff'] = arg[8:] 1198 elif arg.startswith('random='): opts['seed'] = int(arg[8:])1199 1198 elif arg.startswith('title='): opts['title'] = arg[7:] 1200 1199 elif arg.startswith('data='): opts['datafile'] = arg[6:] 1201 1200 elif arg.startswith('engine='): opts['engine'] = arg[8:] 1202 1201 elif arg.startswith('neval='): opts['count'] = arg[7:] 1202 elif arg.startswith('random='): 1203 opts['seed'] = int(arg[8:]) 1204 opts['sets'] = 0 1205 elif arg == 'random': 1206 opts['seed'] = np.random.randint(1000000) 1207 opts['sets'] = 0 1203 1208 elif arg.startswith('sphere'): 1204 1209 opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150 1205 1210 opts['is2d'] = True 1206 elif arg == 'random': opts['seed'] = np.random.randint(1000000)1207 1211 elif arg == 'preset': opts['seed'] = 1 1208 1212 elif arg == 'mono': opts['mono'] = True … … 1309 1313 # Set the integration parameters to the half sphere 1310 1314 opts['values'].extend([ 1311 'theta=90',1315 #'theta=90', 1312 1316 'theta_pd=%g'%(90/np.sqrt(3)), 1313 1317 'theta_pd_n=%d'%steps, 1314 1318 'theta_pd_type=rectangle', 1315 'phi=0',1319 #'phi=0', 1316 1320 'phi_pd=%g'%(180/np.sqrt(3)), 1317 1321 'phi_pd_n=%d'%(2*steps), … … 1320 1324 ]) 1321 1325 if 'psi' in opts['info'][0].parameters: 1322 opts['values'].append('psi=0') 1326 #opts['values'].append('psi=0') 1327 pass 1323 1328 1324 1329 def parse_pars(opts, maxdim=np.inf):
Note: See TracChangeset
for help on using the changeset viewer.