Changeset 242b361 in sasmodels
 Timestamp:
 Mar 18, 2019 4:52:27 PM (4 years ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 66dbbfb
 Parents:
 d8eaa3d
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/asymint.py
rc462169 r242b361 27 27 import os, sys 28 28 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 29 import warnings 29 30 30 31 import numpy as np … … 36 37 37 38 import sasmodels.special as sp 38 39 # Need to parse shape early since it determines the kernel function40 # that will be used for the various integrators41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]43 39 44 40 DTYPE = 'd' … … 201 197 return norm, Fq 202 198 203 if shape == 'sphere': 204 RADIUS = 50 # integer for the sake of mpf 205 NORM, KERNEL = make_sphere(radius=RADIUS) 206 NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 207 elif shape == 'cylinder': 208 #RADIUS, LENGTH = 10, 100000 209 RADIUS, LENGTH = 10, 300 # integer for the sake of mpf 210 NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 211 NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 212 elif shape == 'triaxial_ellipsoid': 213 #A, B, C = 4450, 14000, 47 214 A, B, C = 445, 140, 47 # integer for the sake of mpf 215 NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 216 NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 217 elif shape == 'parallelepiped': 218 #A, B, C = 4450, 14000, 47 219 A, B, C = 445, 140, 47 # integer for the sake of mpf 220 NORM, KERNEL = make_parallelepiped(A, B, C) 221 NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 222 elif shape == 'core_shell_parallelepiped': 223 #A, B, C = 4450, 14000, 47 224 #A, B, C = 445, 140, 47 # integer for the sake of mpf 225 A, B, C = 114, 1380, 6800 226 DA, DB, DC = 21, 58, 2300 227 SLDA, SLDB, SLDC = "5", "0.3", "11.5" 228 ## default parameters from sasmodels 229 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 230 ## swap ABC to CBA 231 #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 232 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 233 #SLD_SOLVENT,CONTRAST = 0, 4 234 if 1: # C shortest 235 B, C = C, B 236 DB, DC = DC, DB 237 SLDB, SLDC = SLDC, SLDB 238 elif 0: # C longest 239 A, C = C, A 240 DA, DC = DC, DA 241 SLDA, SLDC = SLDC, SLDA 242 #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 243 NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 244 NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 245 elif shape == 'paracrystal': 246 LATTICE = 'bcc' 247 #LATTICE = 'fcc' 248 #LATTICE = 'sc' 249 DNN, D_FACTOR = 220, '0.06' # mpmath needs to initialize floats from string 250 RADIUS = 40 # integer for the sake of mpf 251 NORM, KERNEL = make_paracrystal( 252 radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 253 NORM_MP, KERNEL_MP = make_paracrystal( 254 radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 255 else: 256 raise ValueError("Unknown shape %r"%shape) 199 NORM = 1.0 # type: float 200 KERNEL = None # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray] 201 NORM_MP = 1 # type: mpf 202 KERNEL = None # type: CALLABLE[[mpf, mpf, mpf], mpf] 203 204 SHAPES = [ 205 'sphere', 206 'cylinder', 207 'triaxial_ellipsoid', 208 'parallelepiped', 209 'core_shell_parallelepiped', 210 'fcc_paracrystal', 211 'bcc_paracrystal', 212 'sc_paracrystal', 213 ] 214 def build_shape(shape, **pars): 215 global NORM, KERNEL 216 global NORM_MP, KERNEL_MP 217 218 # Note: using integer or string defaults for the sake of mpf 219 if shape == 'sphere': 220 RADIUS = pars.get('radius', 50) 221 NORM, KERNEL = make_sphere(radius=RADIUS) 222 NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 223 elif shape == 'cylinder': 224 #RADIUS, LENGTH = 10, 100000 225 RADIUS = pars.get('radius', 10) 226 LENGTH = pars.get('radius', 300) 227 NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 228 NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 229 elif shape == 'triaxial_ellipsoid': 230 #A, B, C = 4450, 14000, 47 231 A = pars.get('a', 445) 232 B = pars.get('b', 140) 233 C = pars.get('c', 47) 234 NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 235 NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 236 elif shape == 'parallelepiped': 237 #A, B, C = 4450, 14000, 47 238 A = pars.get('a', 445) 239 B = pars.get('b', 140) 240 C = pars.get('c', 47) 241 NORM, KERNEL = make_parallelepiped(A, B, C) 242 NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 243 elif shape == 'core_shell_parallelepiped': 244 #A, B, C = 4450, 14000, 47 245 #A, B, C = 445, 140, 47 # integer for the sake of mpf 246 A = pars.get('a', 114) 247 B = pars.get('b', 1380) 248 C = pars.get('c', 6800) 249 DA = pars.get('da', 21) 250 DB = pars.get('db', 58) 251 DC = pars.get('dc', 2300) 252 SLDA = pars.get('slda', "5") 253 SLDB = pars.get('sldb', "0.3") 254 SLDC = pars.get('sldc', "11.5") 255 ## default parameters from sasmodels 256 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 257 ## swap ABC to CBA 258 #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 259 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 260 #SLD_SOLVENT,CONTRAST = 0, 4 261 if 1: # C shortest 262 B, C = C, B 263 DB, DC = DC, DB 264 SLDB, SLDC = SLDC, SLDB 265 elif 0: # C longest 266 A, C = C, A 267 DA, DC = DC, DA 268 SLDA, SLDC = SLDC, SLDA 269 #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 270 NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 271 NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 272 elif shape.endswith('paracrystal'): 273 LATTICE, _ = shape.split('_') 274 DNN = pars.get('dnn', 220) 275 D_FACTOR = pars.get('d_factor', '0.06') 276 RADIUS = pars.get('radius', 40) 277 NORM, KERNEL = make_paracrystal( 278 radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 279 NORM_MP, KERNEL_MP = make_paracrystal( 280 radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 281 else: 282 raise ValueError("Unknown shape %r"%shape) 257 283 258 284 # Note: hardcoded in mp_quad … … 273 299 274 300 # 2D integration functions 275 def mp_quad_2d(q , shape):301 def mp_quad_2d(q): 276 302 evals = [0] 277 303 def integrand(theta, phi): … … 393 419 print("simpson%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 394 420 print("romb%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 421 422 def quadpy_method(q, rule): 423 """ 424 Use *rule*="name:index" where name and index are chosen from below. 425 426 Available rule names and the corresponding indices:: 427 428 AlbrechtCollatz: [15] 429 BazantOh: 9, 11, 13 430 HeoXu: 13, 15, 17, 19[12], 21[16], 23[13], 25[12], 27[13], 431 29, 31, 33, 35, 37, 39[12] 432 FliegeMaier: 4, 9, 16, 25 433 Lebedev: 3[ac], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 434 41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 435 McLaren: [110] 436 Stroud: U3 31, U3 5[15], U3 7[12], U3 81, U3 9[13], 437 U3 11[13], U3 141 438 """ 439 try: 440 import quadpy 441 except ImportError: 442 warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests") 443 return 444 445 from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu, 446 FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical) 447 RULES = { 448 'AlbrechtCollatz': AlbrechtCollatz, 449 'BazantOh': BazantOh, 450 'HeoXu': HeoXu, 451 'FliegeMaier': FliegeMaier, 452 'Lebedev': Lebedev, 453 'McLaren': McLaren, 454 'Stroud': Stroud, 455 } 456 int_index = 'AlbrechtCollatz', 'McLaren' 457 458 rule_name, rule_index = rule.split(':') 459 index = int(rule_index) if rule_name in int_index else rule_index 460 rule_obj = RULES[rule_name](index) 461 fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal) 462 Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi) 463 print("%s degree=%d points=%s => %g" 464 % (rule, rule_obj.degree, len(rule_obj.points), Iq)) 395 465 396 466 def plot_2d(q, n=300): … … 414 484 pylab.show() 415 485 416 def main(Qstr): 417 Q = float(Qstr) 418 if shape == 'sphere': 486 def main(): 487 import argparse 488 489 parser = argparse.ArgumentParser( 490 description="asymmetric integration explorer", 491 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 492 ) 493 parser.add_argument('s', 'shape', choices=SHAPES, 494 default='parallelepiped', 495 help='oriented shape') 496 parser.add_argument('q', 'q_value', type=str, default='0.005', 497 help='Q value to evaluate') 498 parser.add_argument('pars', type=str, nargs='*', default=[], 499 help='p=val for p in shape parameters') 500 opts = parser.parse_args() 501 pars = {k: v for par in opts.pars for k, v in [par.split('=')]} 502 build_shape(opts.shape, **pars) 503 504 Q = float(opts.q_value) 505 if opts.shape == 'sphere': 419 506 print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2) 507 508 # Methods from quadpy, if quadpy is available 509 # AlbrechtCollatz: [15] 510 # BazantOh: 9, 11, 13 511 # HeoXu: 13, 15, 17, 19[12], 21[16], 23[13], 25[12], 27[13], 512 # 29, 31, 33, 35, 37, 39[12] 513 # FliegeMaier: 4, 9, 16, 25 514 # Lebedev: 3[ac], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 515 # 41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 516 # McLaren: [110] 517 # Stroud: U3 31, U3 5[15], U3 7[12], U3 81, U3 9[13], 518 # U3 11[13], U3 141 519 quadpy_method(Q, "AlbrechtCollatz:5") 520 quadpy_method(Q, "HeoXu:392") 521 quadpy_method(Q, "FliegeMaier:25") 522 quadpy_method(Q, "Lebedev:19") 523 quadpy_method(Q, "Lebedev:131") 524 quadpy_method(Q, "McLaren:10") 525 quadpy_method(Q, "Stroud:U3 141") 526 420 527 print("gauss20", *gauss_quad_2d(Q, n=20)) 421 528 print("gauss76", *gauss_quad_2d(Q, n=76)) … … 427 534 print("gauss76 usub", *gauss_quad_usub(Q, n=76)) 428 535 print("gauss150 usub", *gauss_quad_usub(Q, n=150)) 536 429 537 #gridded_2d(Q, n=2**8+1) 430 538 gridded_2d(Q, n=2**10+1) 431 539 #gridded_2d(Q, n=2**12+1) 432 540 #gridded_2d(Q, n=2**15+1) 433 if shape not in ('paracrystal', 'core_shell_parallelepiped'): 434 # adaptive forms on models for which the calculations are fast enough 541 # adaptive forms on models for which the calculations are fast enough 542 SLOW_SHAPES = { 543 'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 544 'core_shell_parallelepiped', 545 } 546 if opts.shape not in SLOW_SHAPES: 435 547 print("dblquad", *scipy_dblquad_2d(Q)) 436 548 print("semiromberg100", *semi_romberg_2d(Q, n=100)) 437 549 print("romberg", *scipy_romberg_2d(Q)) 438 550 with mp.workprec(100): 439 print("mpmath", *mp_quad_2d(mp.mpf( Qstr), shape))551 print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value))) 440 552 plot_2d(Q, n=200) 441 553 442 554 if __name__ == "__main__": 443 main( Qstr)555 main()
Note: See TracChangeset
for help on using the changeset viewer.