Changeset 31eea1f in sasmodels for explore


Ignore:
Timestamp:
Oct 22, 2017 7:00:50 PM (7 years 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:
f728001
Parents:
2bccb5a
Message:

explore accuracy of different 1D integration schemes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • explore/asymint.py

    r4f611f1 r31eea1f  
    99 
    1010import numpy as np 
     11import mpmath as mp 
    1112from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     13from numpy.polynomial.legendre import leggauss 
    1214from scipy.integrate import dblquad, simps, romb, romberg 
    1315import pylab 
    … …  
    1921from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x 
    2022 
    21 SLD = 3.0 
    22 SLD_SOLVENT = 6.3 
     23def mp_3j1x_x(x): 
     24    return 3*(mp.sin(x)/x - mp.cos(x))/(x*x) 
     25def mp_2J1x_x(x): 
     26    return 2*mp.j1(x)/x 
     27def mp_sinx_x(x): 
     28    return mp.sin(x)/x 
     29 
     30SLD = 3 
     31SLD_SOLVENT = 6 
    2332CONTRAST = SLD - SLD_SOLVENT 
    2433 
     34def 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 
     44def 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 
     55def 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 
     63def 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 
    2572def make_cylinder(radius, length): 
    26     def cylinder(qa, qb, qc): 
     73    def Fq(qa, qb, qc): 
    2774        qab = sqrt(qa**2 + qb**2) 
    2875        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length) 
    2976    volume = pi*radius**2*length 
    30     norm = 1e-4*volume*CONTRAST**2 
    31     return norm, cylinder 
     77    norm = volume*CONTRAST**2/10**4 
     78    return norm, Fq 
     79 
     80def 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 
    3288 
    3389def 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) 
    3792        return sas_3j1x_x(q*radius) 
    3893    volume = 4*pi*radius**3/3 
    39     norm = 1e-4*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 
     97def 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 
     106shape = 'parallelepiped' 
     107#shape = 'triellip' 
     108#shape = 'sphere' 
     109#shape = 'cylinder' 
     110if 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) 
     115elif 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) 
     120elif 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) 
     125elif 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) 
     129else: 
     130    raise ValueError("Unknown shape %r"%shape) 
    45131 
    46132THETA_LOW, THETA_HIGH = 0, pi 
    … …  
    48134SCALE = 1 
    49135 
     136# mathematica code for triaxial_ellipsoid (untested) 
     137_ = """ 
     138R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2] 
     139Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r) 
     140V[a_, b_, c_] := 4/3 pi a b c 
     141Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2 
     142F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]] 
     143I[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}] 
     144I[6/10^3, 63/10, 3, 445, 140, 47] 
     145""" 
     146 
     147 
     148def 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 
    50161 
    51162def kernel(q, theta, phi): 
    … …  
    59170    return NORM*KERNEL(qa, qb, qc)**2 
    60171 
    61  
    62172def scipy_dblquad(q): 
    63173    """ 
    … …  
    66176    """ 
    67177    evals = [0] 
    68     def integrand(theta, phi): 
     178    def integrand(phi, theta): 
    69179        evals[0] += 1 
    70180        Zq = kernel(q, theta=theta, phi=phi) 
    71181        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) 
    75184 
    76185 
    … …  
    80189    complete in a reasonable time.  No idea if it is accurate. 
    81190    """ 
     191    evals = [0] 
    82192    def inner(phi, theta): 
     193        evals[0] += 1 
    83194        return kernel(q, theta=theta, phi=phi) 
    84195    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) 
    87200 
    88201 
    … …  
    96209        return kernel(q, theta=theta, phi=phi) 
    97210    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) 
    103214 
    104215def gauss_quad(q, n=150): 
    … …  
    110221    elif n == 76: 
    111222        z, w = Gauss76Z, Gauss76Wt 
     223    elif n == 150: 
     224        z, w = Gauss150Z, Gauss150Wt 
    112225    else: 
    113         z, w = Gauss150Z, Gauss150Wt 
     226        z, w = leggauss(n) 
    114227    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW 
    115228    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW 
    116229    Atheta, Aphi = np.meshgrid(theta, phi) 
    117230    Aw = w[None, :] * w[:, None] 
    118     sin_theta = np.fmax(abs(sin(Atheta)), 1e-6) 
     231    sin_theta = abs(sin(Atheta)) 
    119232    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_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2 
     235    Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch 
     236    return n**2, Iq 
    122237 
    123238 
    … …  
    134249    Zq *= abs(sin(Atheta)) 
    135250    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)) 
    141255 
    142256def plot(q, n=300): 
    … …  
    161275 
    162276if __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("gauss-20", *gauss_quad(Q, n=20)) 
     284    print("gauss-76", *gauss_quad(Q, n=76)) 
     285    print("gauss-150", *gauss_quad(Q, n=150)) 
     286    print("gauss-500", *gauss_quad(Q, n=500)) 
     287    print("dblquad", *scipy_dblquad(Q)) 
     288    print("semi-romberg-100", *semi_romberg(Q, n=100)) 
     289    print("romberg", *scipy_romberg_2d(Q)) 
    169290    #gridded_integrals(Q, n=2**8+1) 
    170     #gridded_integrals(Q, n=2**10+1) 
     291    gridded_integrals(Q, n=2**10+1) 
    171292    #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) 
Note: See TracChangeset for help on using the changeset viewer.