source: sasmodels/explore/asymint.py @ 36b3154

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 36b3154 was 31eea1f, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

explore accuracy of different 1D integration schemes

  • Property mode set to 100644
File size: 9.5 KB
Line 
1"""
2Asymmetric shape integration
3"""
4
5from __future__ import print_function, division
6
7import os, sys
8sys.path.insert(0, os.path.dirname(os.path.dirname(__file__)))
9
10import numpy as np
11import mpmath as mp
12from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
13from numpy.polynomial.legendre import leggauss
14from scipy.integrate import dblquad, simps, romb, romberg
15import pylab
16
17from sasmodels.special import square
18from sasmodels.special import Gauss20Wt, Gauss20Z
19from sasmodels.special import Gauss76Wt, Gauss76Z
20from sasmodels.special import Gauss150Wt, Gauss150Z
21from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x
22
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
32CONTRAST = SLD - SLD_SOLVENT
33
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
72def make_cylinder(radius, length):
73    def Fq(qa, qb, qc):
74        qab = sqrt(qa**2 + qb**2)
75        return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length)
76    volume = pi*radius**2*length
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
88
89def make_sphere(radius):
90    def Fq(qa, qb, qc):
91        q = sqrt(qa**2 + qb**2 + qc**2)
92        return sas_3j1x_x(q*radius)
93    volume = 4*pi*radius**3/3
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)
131
132THETA_LOW, THETA_HIGH = 0, pi
133PHI_LOW, PHI_HIGH = 0, 2*pi
134SCALE = 1
135
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
161
162def kernel(q, theta, phi):
163    """
164    S(q) kernel for paracrystal forms.
165    """
166    qab = q*sin(theta)
167    qa = qab*cos(phi)
168    qb = qab*sin(phi)
169    qc = q*cos(theta)
170    return NORM*KERNEL(qa, qb, qc)**2
171
172def scipy_dblquad(q):
173    """
174    Compute the integral using scipy dblquad.  This gets the correct answer
175    eventually, but it is slow.
176    """
177    evals = [0]
178    def integrand(phi, theta):
179        evals[0] += 1
180        Zq = kernel(q, theta=theta, phi=phi)
181        return Zq*sin(theta)
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)
184
185
186def scipy_romberg_2d(q):
187    """
188    Compute the integral using romberg integration.  This function does not
189    complete in a reasonable time.  No idea if it is accurate.
190    """
191    evals = [0]
192    def inner(phi, theta):
193        evals[0] += 1
194        return kernel(q, theta=theta, phi=phi)
195    def outer(theta):
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)
200
201
202def semi_romberg(q, n=100):
203    """
204    Use 1D romberg integration in phi and regular simpsons rule in theta.
205    """
206    evals = [0]
207    def inner(phi, theta):
208        evals[0] += 1
209        return kernel(q, theta=theta, phi=phi)
210    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
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)
214
215def gauss_quad(q, n=150):
216    """
217    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
218    """
219    if n == 20:
220        z, w = Gauss20Z, Gauss20Wt
221    elif n == 76:
222        z, w = Gauss76Z, Gauss76Wt
223    elif n == 150:
224        z, w = Gauss150Z, Gauss150Wt
225    else:
226        z, w = leggauss(n)
227    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
228    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
229    Atheta, Aphi = np.meshgrid(theta, phi)
230    Aw = w[None, :] * w[:, None]
231    sin_theta = abs(sin(Atheta))
232    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
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
237
238
239def gridded_integrals(q, n=300):
240    """
241    Compute the integral on a regular grid using rectangular, trapezoidal,
242    simpsons, and romberg integration.  Romberg integration requires that
243    the grid be of size n = 2**k + 1.
244    """
245    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
246    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
247    Atheta, Aphi = np.meshgrid(theta, phi)
248    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
249    Zq *= abs(sin(Atheta))
250    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
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))
255
256def plot(q, n=300):
257    """
258    Plot the 2D surface that needs to be integrated in order to compute
259    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
260    of points in the grid.
261    """
262    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
263    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
264    Atheta, Aphi = np.meshgrid(theta, phi)
265    Zq = kernel(q=q, theta=Atheta, phi=Aphi)
266    #Zq *= abs(sin(Atheta))
267    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
268    pylab.axis('tight')
269    pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q))
270    pylab.xlabel("theta (degrees)")
271    pylab.ylabel("phi (degrees)")
272    cbar = pylab.colorbar()
273    cbar.set_label('log10 S(q)')
274    pylab.show()
275
276if __name__ == "__main__":
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))
290    #gridded_integrals(Q, n=2**8+1)
291    gridded_integrals(Q, n=2**10+1)
292    #gridded_integrals(Q, n=2**13+1)
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 TracBrowser for help on using the repository browser.