Changeset 242b361 in sasmodels


Ignore:
Timestamp:
Mar 18, 2019 4:52:27 PM (4 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
66dbbfb
Parents:
d8eaa3d
Message:

add quadpy.sphere methods to asymmetric integration tests

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/asymint.py

    rc462169 r242b361  
    2727import os, sys 
    2828sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     29import warnings 
    2930 
    3031import numpy as np 
     
    3637 
    3738import sasmodels.special as sp 
    38  
    39 # Need to parse shape early since it determines the kernel function 
    40 # that will be used for the various integrators 
    41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
    4339 
    4440DTYPE = 'd' 
     
    201197    return norm, Fq 
    202198 
    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 A-B-C to C-B-A 
    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) 
     199NORM = 1.0  # type: float 
     200KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray] 
     201NORM_MP = 1  # type: mpf 
     202KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf] 
     203 
     204SHAPES = [ 
     205    'sphere', 
     206    'cylinder', 
     207    'triaxial_ellipsoid', 
     208    'parallelepiped', 
     209    'core_shell_parallelepiped', 
     210    'fcc_paracrystal', 
     211    'bcc_paracrystal', 
     212    'sc_paracrystal', 
     213] 
     214def 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 A-B-C to C-B-A 
     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) 
    257283 
    258284# Note: hardcoded in mp_quad 
     
    273299 
    274300# 2D integration functions 
    275 def mp_quad_2d(q, shape): 
     301def mp_quad_2d(q): 
    276302    evals = [0] 
    277303    def integrand(theta, phi): 
     
    393419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
    394420    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
     421 
     422def 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: [1-5] 
     429        BazantOh: 9, 11, 13 
     430        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     431            29, 31, 33, 35, 37, 39-[1-2] 
     432        FliegeMaier: 4, 9, 16, 25 
     433        Lebedev: 3[a-c], 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: [1-10] 
     436        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     437            U3 11-[1-3], U3 14-1 
     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)) 
    395465 
    396466def plot_2d(q, n=300): 
     
    414484    pylab.show() 
    415485 
    416 def main(Qstr): 
    417     Q = float(Qstr) 
    418     if shape == 'sphere': 
     486def 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': 
    419506        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2) 
     507 
     508    # Methods from quadpy, if quadpy is available 
     509    #  AlbrechtCollatz: [1-5] 
     510    #  BazantOh: 9, 11, 13 
     511    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     512    #     29, 31, 33, 35, 37, 39-[1-2] 
     513    #  FliegeMaier: 4, 9, 16, 25 
     514    #  Lebedev: 3[a-c], 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: [1-10] 
     517    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     518    #     U3 11-[1-3], U3 14-1 
     519    quadpy_method(Q, "AlbrechtCollatz:5") 
     520    quadpy_method(Q, "HeoXu:39-2") 
     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 14-1") 
     526 
    420527    print("gauss-20", *gauss_quad_2d(Q, n=20)) 
    421528    print("gauss-76", *gauss_quad_2d(Q, n=76)) 
     
    427534    print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
    428535    print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
     536 
    429537    #gridded_2d(Q, n=2**8+1) 
    430538    gridded_2d(Q, n=2**10+1) 
    431539    #gridded_2d(Q, n=2**12+1) 
    432540    #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: 
    435547        print("dblquad", *scipy_dblquad_2d(Q)) 
    436548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
    437549        print("romberg", *scipy_romberg_2d(Q)) 
    438550        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))) 
    440552    plot_2d(Q, n=200) 
    441553 
    442554if __name__ == "__main__": 
    443     main(Qstr) 
     555    main() 
Note: See TracChangeset for help on using the changeset viewer.