source: sasmodels/explore/asymint.py @ f691fac

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since f691fac was 66dbbfb, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

prettier output for explore/asymint.py
  • Property mode set to 100755
File size: 20.1 KB
Line 
1#!/usr/bin/env python
2"""
3Asymmetric shape integration
4
5Usage:
6
7    explore/asymint.py [MODEL] [q-value]
8
9Computes the numerical integral over theta and phi of the given model at a
10single point q using different algorithms or the same algorithm with different
11precision.  It also displays a 2-D image of the theta-phi surface that is
12being integrated.
13
14The available models are:
15
16    triaxial_ellipsoid, parallelpiped, paracrystal, cylinder, sphere
17
18Cylinder and sphere are included as simple checks on the integration
19algorithms. Cylinder is better investigated using 1-D integration methods in
20explore/symint.py.  Sphere has an easily computed analytic value which is
21identical for all theta-phi for a given q, so it is useful for checking
22that the normalization constants are correct for the different algorithms.
23"""
24
25from __future__ import print_function, division
26
27import os, sys
28sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
29import warnings
30
31import numpy as np
32import mpmath as mp
33from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos
34from numpy.polynomial.legendre import leggauss
35from scipy.integrate import dblquad, simps, romb, romberg
36import pylab
37
38import sasmodels.special as sp
39
40DTYPE = 'd'
41
42class MPenv:
43    sqrt = staticmethod(mp.sqrt)
44    exp = staticmethod(mp.exp)
45    expm1 = staticmethod(mp.expm1)
46    cos = staticmethod(mp.cos)
47    sin = staticmethod(mp.sin)
48    tan = staticmethod(mp.tan)
49    @staticmethod
50    def sas_3j1x_x(x):
51        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
52    @staticmethod
53    def sas_2J1x_x(x):
54        return 2*mp.j1(x)/x
55    @staticmethod
56    def sas_sinx_x(x):
57        return mp.sin(x)/x
58    pi = mp.pi
59    mpf = staticmethod(mp.mpf)
60
61class NPenv:
62    sqrt = staticmethod(np.sqrt)
63    exp = staticmethod(np.exp)
64    expm1 = staticmethod(np.expm1)
65    cos = staticmethod(np.cos)
66    sin = staticmethod(np.sin)
67    tan = staticmethod(np.tan)
68    sas_3j1x_x = staticmethod(sp.sas_3j1x_x)
69    sas_2J1x_x = staticmethod(sp.sas_2J1x_x)
70    sas_sinx_x = staticmethod(sp.sas_sinx_x)
71    pi = np.pi
72    #mpf = staticmethod(float)
73    mpf = staticmethod(lambda x: np.array(x, DTYPE))
74
75SLD = 3
76SLD_SOLVENT = 6
77CONTRAST = SLD - SLD_SOLVENT
78
79# Carefully code models so that mpmath will use full precision.  That means:
80#    * wrap inputs in env.mpf
81#    * don't use floating point constants, only integers
82#    * for division, make sure the numerator or denominator is env.mpf
83#    * use env.pi, env.sas_sinx_x, etc. for functions
84def make_parallelepiped(a, b, c, env=NPenv):
85    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
86    def Fq(qa, qb, qc):
87        siA = env.sas_sinx_x(a*qa/2)
88        siB = env.sas_sinx_x(b*qb/2)
89        siC = env.sas_sinx_x(c*qc/2)
90        return siA * siB * siC
91    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
92    volume = a*b*c
93    norm = CONTRAST**2*volume/10000
94    return norm, Fq
95
96def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):
97    overlapping = False
98    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
99    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)
100    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)
101    dr0 = CONTRAST
102    drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT
103    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
104    def Fq(qa, qb, qc):
105        siA = a*env.sas_sinx_x(a*qa/2)
106        siB = b*env.sas_sinx_x(b*qb/2)
107        siC = c*env.sas_sinx_x(c*qc/2)
108        siAt = tA*env.sas_sinx_x(tA*qa/2)
109        siBt = tB*env.sas_sinx_x(tB*qb/2)
110        siCt = tC*env.sas_sinx_x(tC*qc/2)
111        if overlapping:
112            return (dr0*siA*siB*siC
113                    + drA*(siAt-siA)*siB*siC
114                    + drB*siAt*(siBt-siB)*siC
115                    + drC*siAt*siBt*(siCt-siC))
116        else:
117            return (dr0*siA*siB*siC
118                    + drA*(siAt-siA)*siB*siC
119                    + drB*siA*(siBt-siB)*siC
120                    + drC*siA*siB*(siCt-siC))
121    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)
122    if overlapping:
123        volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc
124    else:
125        volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc
126    norm = 1/(volume*10000)
127    return norm, Fq
128
129def make_triaxial_ellipsoid(a, b, c, env=NPenv):
130    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
131    def Fq(qa, qb, qc):
132        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
133        return env.sas_3j1x_x(qr)
134    Fq.__doc__ = "triaxial ellipsoid minor=%g, major=%g polar=%g"%(a, b, c)
135    volume = 4*env.pi*a*b*c/3
136    norm = CONTRAST**2*volume/10000
137    return norm, Fq
138
139def make_cylinder(radius, length, env=NPenv):
140    radius, length = env.mpf(radius), env.mpf(length)
141    def Fq(qa, qb, qc):
142        qab = env.sqrt(qa**2 + qb**2)
143        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2)
144    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
145    volume = env.pi*radius**2*length
146    norm = CONTRAST**2*volume/10000
147    return norm, Fq
148
149def make_sphere(radius, env=NPenv):
150    radius = env.mpf(radius)
151    def Fq(qa, qb, qc):
152        q = env.sqrt(qa**2 + qb**2 + qc**2)
153        return env.sas_3j1x_x(q*radius)
154    Fq.__doc__ = "sphere radius=%g"%(radius, )
155    volume = 4*pi*radius**3
156    norm = CONTRAST**2*volume/10000
157    return norm, Fq
158
159def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv):
160    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
161    def sc(qa, qb, qc):
162        return qa, qb, qc
163    def bcc(qa, qb, qc):
164        a1 = (+qa + qb + qc)/2
165        a2 = (-qa - qb + qc)/2
166        a3 = (-qa + qb - qc)/2
167        return a1, a2, a3
168    def fcc(qa, qb, qc):
169        a1 = ( 0 + qb + qc)/2
170        a2 = (-qa + 0 + qc)/2
171        a3 = (-qa + qb + 0)/2
172        return a1, a2, a3
173    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice]
174    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
175    def Fq(qa, qb, qc):
176        a1, a2, a3 = lattice_fn(qa, qb, qc)
177        # Note: paper says that different directions can have different
178        # distoration factors.  Easy enough to add to the code.
179        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2
180        exp_arg = env.exp(arg)
181        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)]
182        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2])
183
184        q = env.sqrt(qa**2 + qb**2 + qc**2)
185        Fq = env.sas_3j1x_x(q*radius)
186        # the caller computes F(q)**2, but we need it to compute S(q)*F(q)**2
187        return env.sqrt(Sq)*Fq
188    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius)
189    def sphere_volume(r): return 4*env.pi*r**3/3
190    Vf = {
191        'sc': sphere_volume(radius/dnn),
192        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn),
193        'fcc': 4*sphere_volume(1/env.sqrt(2)*radius/dnn),
194    }[lattice]
195    volume = sphere_volume(radius)
196    norm = CONTRAST**2*volume/10000*Vf
197    return norm, Fq
198
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)
283
284# Note: hardcoded in mp_quad
285THETA_LOW, THETA_HIGH = 0, pi
286PHI_LOW, PHI_HIGH = 0, 2*pi
287SCALE = 1
288
289# mathematica code for triaxial_ellipsoid (untested)
290_ = """
291R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2]
292Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r)
293V[a_, b_, c_] := 4/3 pi a b c
294Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2
295F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]]
296I[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}]
297I[6/10^3, 63/10, 3, 445, 140, 47]
298"""
299
300# 2D integration functions
301def mp_quad_2d(q):
302    evals = [0]
303    def integrand(theta, phi):
304        evals[0] += 1
305        qab = q*mp.sin(theta)
306        qa = qab*mp.cos(phi)
307        qb = qab*mp.sin(phi)
308        qc = q*mp.cos(theta)
309        Zq = KERNEL_MP(qa, qb, qc)**2
310        return Zq*mp.sin(theta)
311    ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi))
312    Iq = NORM_MP*ans/(4*mp.pi)
313    return evals[0], Iq
314
315def kernel_2d(q, theta, phi):
316    """
317    S(q) kernel for paracrystal forms.
318    """
319    qab = q*sin(theta)
320    qa = qab*cos(phi)
321    qb = qab*sin(phi)
322    qc = q*cos(theta)
323    return NORM*KERNEL(qa, qb, qc)**2
324
325def scipy_dblquad_2d(q):
326    """
327    Compute the integral using scipy dblquad.  This gets the correct answer
328    eventually, but it is slow.
329    """
330    evals = [0]
331    def integrand(phi, theta):
332        evals[0] += 1
333        Zq = kernel_2d(q, theta=theta, phi=phi)
334        return Zq*sin(theta)
335    ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]
336    return evals[0], ans*SCALE/(4*pi)
337
338
339def scipy_romberg_2d(q):
340    """
341    Compute the integral using romberg integration.  This function does not
342    complete in a reasonable time.  No idea if it is accurate.
343    """
344    evals = [0]
345    def inner(phi, theta):
346        evals[0] += 1
347        return kernel_2d(q, theta=theta, phi=phi)
348    def outer(theta):
349        Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))
350        return Zq*sin(theta)
351    ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)
352    return evals[0], ans*SCALE/(4*pi)
353
354
355def semi_romberg_2d(q, n=100):
356    """
357    Use 1D romberg integration in phi and regular simpsons rule in theta.
358    """
359    evals = [0]
360    def inner(phi, theta):
361        evals[0] += 1
362        return kernel_2d(q, theta=theta, phi=phi)
363    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
364    Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta]
365    ans = simps(np.array(Zq)*sin(theta), dx=theta[1]-theta[0])
366    return evals[0], ans*SCALE/(4*pi)
367
368def gauss_quad_2d(q, n=150):
369    """
370    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
371    """
372    z, w = leggauss(n)
373    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
374    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
375    Atheta, Aphi = np.meshgrid(theta, phi)
376    Aw = w[None, :] * w[:, None]
377    sin_theta = abs(sin(Atheta))
378    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
379    # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range
380    dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2
381    Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch
382    return n**2, Iq
383
384def gauss_quad_usub(q, n=150, dtype=DTYPE):
385    """
386    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
387
388    Use *u = sin theta* substitution, and restrict integration over a single
389    quadrant for shapes that are mirror symmetric about AB, AC and BC planes.
390
391    Note that this doesn't work for fcc/bcc paracrystals, which instead step
392    over the entire 4 pi surface uniformly in theta-phi.
393    """
394    z, w = leggauss(n)
395    cos_theta = 0.5 * (z + 1)
396    theta = arccos(cos_theta)
397    phi = pi/2*(0.5 * (z + 1))
398    Atheta, Aphi = np.meshgrid(theta, phi)
399    Aw = w[None, :] * w[:, None]
400    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)]
401    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
402    Iq = np.sum(Zq*Aw)*0.25
403    return n**2, Iq
404
405def gridded_2d(q, n=300):
406    """
407    Compute the integral on a regular grid using rectangular, trapezoidal,
408    simpsons, and romberg integration.  Romberg integration requires that
409    the grid be of size n = 2**k + 1.
410    """
411    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
412    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
413    Atheta, Aphi = np.meshgrid(theta, phi)
414    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
415    Zq *= abs(sin(Atheta))
416    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
417    print("rect-%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi))
418    print("trapz-%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
420    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 => %.15g"
464          % (rule, rule_obj.degree, len(rule_obj.points), Iq))
465
466def plot_2d(q, n=300):
467    """
468    Plot the 2D surface that needs to be integrated in order to compute
469    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
470    of points in the grid.
471    """
472    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
473    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
474    Atheta, Aphi = np.meshgrid(theta, phi)
475    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
476    #Zq *= abs(sin(Atheta))
477    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
478    pylab.axis('tight')
479    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q))
480    pylab.xlabel("theta (degrees)")
481    pylab.ylabel("phi (degrees)")
482    cbar = pylab.colorbar()
483    cbar.set_label('log10 S(q)')
484    pylab.show()
485
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':
506        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
527    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20))
528    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76))
529    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150))
530    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500))
531    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025))
532    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049))
533    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20))
534    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76))
535    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150))
536
537    #gridded_2d(Q, n=2**8+1)
538    gridded_2d(Q, n=2**10+1)
539    #gridded_2d(Q, n=2**12+1)
540    #gridded_2d(Q, n=2**15+1)
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:
547        print("dblquad", *scipy_dblquad_2d(Q))
548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
549        print("romberg", *scipy_romberg_2d(Q))
550        with mp.workprec(100):
551            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value)))
552    plot_2d(Q, n=200)
553
554if __name__ == "__main__":
555    main()
Note: See TracBrowser for help on using the repository browser.