source: sasmodels/explore/asymint.py @ d0dc9a3

#!/usr/bin/env python
"""
Asymmetric shape integration

Usage:

    explore/asymint.py [MODEL] [q-value]

Computes the numerical integral over theta and phi of the given model at a
single point q using different algorithms or the same algorithm with different
precision.  It also displays a 2-D image of the theta-phi surface that is
being integrated.

The available models are:

    triaxial_ellipsoid, parallelpiped, paracrystal, cylinder, sphere

Cylinder and sphere are included as simple checks on the integration
algorithms. Cylinder is better investigated using 1-D integration methods in
explore/symint.py.  Sphere has an easily computed analytic value which is
identical for all theta-phi for a given q, so it is useful for checking
that the normalization constants are correct for the different algorithms.
"""

from __future__ import print_function, division

import os, sys
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))

import numpy as np
import mpmath as mp
from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10
from numpy.polynomial.legendre import leggauss
from scipy.integrate import dblquad, simps, romb, romberg
import pylab

import sasmodels.special as sp

# Need to parse shape early since it determines the kernel function
# that will be used for the various integrators
shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]

class MPenv:
    sqrt = staticmethod(mp.sqrt)
    exp = staticmethod(mp.exp)
    expm1 = staticmethod(mp.expm1)
    cos = staticmethod(mp.cos)
    sin = staticmethod(mp.sin)
    tan = staticmethod(mp.tan)
    @staticmethod
    def sas_3j1x_x(x):
        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
    @staticmethod
    def sas_2J1x_x(x):
        return 2*mp.j1(x)/x
    @staticmethod
    def sas_sinx_x(x):
        return mp.sin(x)/x
    pi = mp.pi
    mpf = staticmethod(mp.mpf)

class NPenv:
    sqrt = staticmethod(np.sqrt)
    exp = staticmethod(np.exp)
    expm1 = staticmethod(np.expm1)
    cos = staticmethod(np.cos)
    sin = staticmethod(np.sin)
    tan = staticmethod(np.tan)
    sas_3j1x_x = staticmethod(sp.sas_3j1x_x)
    sas_2J1x_x = staticmethod(sp.sas_2J1x_x)
    sas_sinx_x = staticmethod(sp.sas_sinx_x)
    pi = np.pi
    mpf = staticmethod(float)

SLD = 3
SLD_SOLVENT = 6
CONTRAST = SLD - SLD_SOLVENT

# Carefully code models so that mpmath will use full precision.  That means:
#    * wrap inputs in env.mpf
#    * don't use floating point constants, only integers
#    * for division, make sure the numerator or denominator is env.mpf
#    * use env.pi, env.sas_sinx_x, etc. for functions
def make_parallelepiped(a, b, c, env=NPenv):
    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
    def Fq(qa, qb, qc):
        siA = env.sas_sinx_x(a*qa/2)
        siB = env.sas_sinx_x(b*qb/2)
        siC = env.sas_sinx_x(c*qc/2)
        return siA * siB * siC
    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
    volume = a*b*c
    norm = CONTRAST**2*volume/10000
    return norm, Fq

def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):
    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)
    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)
    drV0 = CONTRAST*a*b*c
    dra, drb, drc = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT
    Aa, Ab, Ac = b*c, a*c, a*b
    Ta, Tb, Tc = a + 2*da, b + 2*db, c + 2*dc
    drVa, drVb, drVc = dra*a*Aa, drb*b*Ab, drc*c*Ac
    drVTa, drVTb, drVTc = dra*Ta*Aa, drb*Tb*Ab, drc*Tc*Ac
    def Fq(qa, qb, qc):
        siA = env.sas_sinx_x(a*qa/2)
        siB = env.sas_sinx_x(b*qb/2)
        siC = env.sas_sinx_x(c*qc/2)
        siAt = env.sas_sinx_x(Ta*qa/2)
        siBt = env.sas_sinx_x(Tb*qb/2)
        siCt = env.sas_sinx_x(Tc*qc/2)
        return (drV0*siA*siB*siC
            + (drVTa*siAt-drVa*siA)*siB*siC
            + siA*(drVTb*siBt-drVb*siB)*siC
            + siA*siB*(drVTc*siCt-drVc*siC))
    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)
    volume = a*b*c + 2*da*Aa + 2*db*Ab + 2*dc*Ac
    norm = 1/(volume*10000)
    return norm, Fq

def make_triaxial_ellipsoid(a, b, c, env=NPenv):
    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
    def Fq(qa, qb, qc):
        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
        return env.sas_3j1x_x(qr)
    Fq.__doc__ = "triaxial ellipsoid minor=%g, major=%g polar=%g"%(a, b, c)
    volume = 4*env.pi*a*b*c/3
    norm = CONTRAST**2*volume/10000
    return norm, Fq

def make_cylinder(radius, length, env=NPenv):
    radius, length = env.mpf(radius), env.mpf(length)
    def Fq(qa, qb, qc):
        qab = env.sqrt(qa**2 + qb**2)
        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2)
    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
    volume = env.pi*radius**2*length
    norm = CONTRAST**2*volume/10000
    return norm, Fq

def make_sphere(radius, env=NPenv):
    radius = env.mpf(radius)
    def Fq(qa, qb, qc):
        q = env.sqrt(qa**2 + qb**2 + qc**2)
        return env.sas_3j1x_x(q*radius)
    Fq.__doc__ = "sphere radius=%g"%(radius, )
    volume = 4*pi*radius**3
    norm = CONTRAST**2*volume/10000
    return norm, Fq

def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv):
    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
    def sc(qa, qb, qc):
        return qa, qb, qc
    def bcc(qa, qb, qc):
        a1 = (+qa + qb + qc)/2
        a2 = (-qa - qb + qc)/2
        a3 = (-qa + qb - qc)/2
        return a1, a2, a3
    def fcc(qa, qb, qc):
        a1 = ( 0 + qb + qc)/2
        a2 = (-qa + 0 + qc)/2
        a3 = (-qa + qb + 0)/2
        return a1, a2, a3
    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice]
    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
    def Fq(qa, qb, qc):
        a1, a2, a3 = lattice_fn(qa, qb, qc)
        # Note: paper says that different directions can have different
        # distoration factors.  Easy enough to add to the code.
        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2
        exp_arg = env.exp(arg)
        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)]
        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2])

        q = env.sqrt(qa**2 + qb**2 + qc**2)
        Fq = env.sas_3j1x_x(q*radius)
        # the caller computes F(q)**2, but we need it to compute S(q)*F(q)**2
        return env.sqrt(Sq)*Fq
    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius)
    def sphere_volume(r): return 4*env.pi*r**3/3
    Vf = {
        'sc': sphere_volume(radius/dnn),
        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn),
        'fcc': 4*sphere_volume(1/env.sqrt(2)*radius/dnn),
    }[lattice]
    volume = sphere_volume(radius)
    norm = CONTRAST**2*volume/10000*Vf
    return norm, Fq

if shape == 'sphere':
    RADIUS = 50  # integer for the sake of mpf
    NORM, KERNEL = make_sphere(radius=RADIUS)
    NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
elif shape == 'cylinder':
    #RADIUS, LENGTH = 10, 100000
    RADIUS, LENGTH = 10, 300  # integer for the sake of mpf
    NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
    NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
elif shape == 'triaxial_ellipsoid':
    #A, B, C = 4450, 14000, 47
    A, B, C = 445, 140, 47  # integer for the sake of mpf
    NORM, KERNEL = make_triaxial_ellipsoid(A, B, C)
    NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv)
elif shape == 'parallelepiped':
    #A, B, C = 4450, 14000, 47
    A, B, C = 445, 140, 47  # integer for the sake of mpf
    NORM, KERNEL = make_parallelepiped(A, B, C)
    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
elif shape == 'core_shell_parallelepiped':
    #A, B, C = 4450, 14000, 47
    #A, B, C = 445, 140, 47  # integer for the sake of mpf
    A, B, C = 6800, 114, 1380
    DA, DB, DC = 2300, 21, 58
    SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
    if 1: # C shortest
        B, C = C, B
        DB, DC = DC, DB
        SLDB, SLDC = SLDC, SLDB
    elif 0: # C longest
        A, C = C, A
        DA, DC = DC, DA
        SLDA, SLDC = SLDC, SLDA
    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
elif shape == 'paracrystal':
    LATTICE = 'bcc'
    #LATTICE = 'fcc'
    #LATTICE = 'sc'
    DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string
    RADIUS = 40  # integer for the sake of mpf
    NORM, KERNEL = make_paracrystal(
        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
    NORM_MP, KERNEL_MP = make_paracrystal(
        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
else:
    raise ValueError("Unknown shape %r"%shape)

# Note: hardcoded in mp_quad
THETA_LOW, THETA_HIGH = 0, pi
PHI_LOW, PHI_HIGH = 0, 2*pi
SCALE = 1

# mathematica code for triaxial_ellipsoid (untested)
_ = """
R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2]
Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r)
V[a_, b_, c_] := 4/3 pi a b c
Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2
F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]]
I[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}]
I[6/10^3, 63/10, 3, 445, 140, 47]
"""

# 2D integration functions
def mp_quad_2d(q, shape):
    evals = [0]
    def integrand(theta, phi):
        evals[0] += 1
        qab = q*mp.sin(theta)
        qa = qab*mp.cos(phi)
        qb = qab*mp.sin(phi)
        qc = q*mp.cos(theta)
        Zq = KERNEL_MP(qa, qb, qc)**2
        return Zq*mp.sin(theta)
    ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi))
    Iq = NORM_MP*ans/(4*mp.pi)
    return evals[0], Iq

def kernel_2d(q, theta, phi):
    """
    S(q) kernel for paracrystal forms.
    """
    qab = q*sin(theta)
    qa = qab*cos(phi)
    qb = qab*sin(phi)
    qc = q*cos(theta)
    return NORM*KERNEL(qa, qb, qc)**2

def scipy_dblquad_2d(q):
    """
    Compute the integral using scipy dblquad.  This gets the correct answer
    eventually, but it is slow.
    """
    evals = [0]
    def integrand(phi, theta):
        evals[0] += 1
        Zq = kernel_2d(q, theta=theta, phi=phi)
        return Zq*sin(theta)
    ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]
    return evals[0], ans*SCALE/(4*pi)


def scipy_romberg_2d(q):
    """
    Compute the integral using romberg integration.  This function does not
    complete in a reasonable time.  No idea if it is accurate.
    """
    evals = [0]
    def inner(phi, theta):
        evals[0] += 1
        return kernel_2d(q, theta=theta, phi=phi)
    def outer(theta):
        Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))
        return Zq*sin(theta)
    ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)
    return evals[0], ans*SCALE/(4*pi)


def semi_romberg_2d(q, n=100):
    """
    Use 1D romberg integration in phi and regular simpsons rule in theta.
    """
    evals = [0]
    def inner(phi, theta):
        evals[0] += 1
        return kernel_2d(q, theta=theta, phi=phi)
    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
    Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta]
    ans = simps(np.array(Zq)*sin(theta), dx=theta[1]-theta[0])
    return evals[0], ans*SCALE/(4*pi)

def gauss_quad_2d(q, n=150):
    """
    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
    """
    z, w = leggauss(n)
    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
    Atheta, Aphi = np.meshgrid(theta, phi)
    Aw = w[None, :] * w[:, None]
    sin_theta = abs(sin(Atheta))
    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
    # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range
    dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2
    Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch
    return n**2, Iq

def gridded_2d(q, n=300):
    """
    Compute the integral on a regular grid using rectangular, trapezoidal,
    simpsons, and romberg integration.  Romberg integration requires that
    the grid be of size n = 2**k + 1.
    """
    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
    Atheta, Aphi = np.meshgrid(theta, phi)
    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
    Zq *= abs(sin(Atheta))
    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
    print("rect-%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi))
    print("trapz-%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))

def plot_2d(q, n=300):
    """
    Plot the 2D surface that needs to be integrated in order to compute
    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
    of points in the grid.
    """
    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
    Atheta, Aphi = np.meshgrid(theta, phi)
    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
    #Zq *= abs(sin(Atheta))
    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
    pylab.axis('tight')
    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q))
    pylab.xlabel("theta (degrees)")
    pylab.ylabel("phi (degrees)")
    cbar = pylab.colorbar()
    cbar.set_label('log10 S(q)')
    pylab.show()

def main(Qstr):
    Q = float(Qstr)
    if shape == 'sphere':
        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2)
    print("gauss-20", *gauss_quad_2d(Q, n=20))
    print("gauss-76", *gauss_quad_2d(Q, n=76))
    print("gauss-150", *gauss_quad_2d(Q, n=150))
    print("gauss-500", *gauss_quad_2d(Q, n=500))
    print("gauss-1025", *gauss_quad_2d(Q, n=1025))
    print("gauss-2049", *gauss_quad_2d(Q, n=2049))
    #gridded_2d(Q, n=2**8+1)
    gridded_2d(Q, n=2**10+1)
    #gridded_2d(Q, n=2**12+1)
    #gridded_2d(Q, n=2**15+1)
    if shape not in ('paracrystal', 'core_shell_parallelepiped'):
        # adaptive forms on models for which the calculations are fast enough
        print("dblquad", *scipy_dblquad_2d(Q))
        print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
        print("romberg", *scipy_romberg_2d(Q))
        with mp.workprec(100):
            print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape))
    plot_2d(Q, n=200)

if __name__ == "__main__":
    main(Qstr)