Changeset 6ceca44 – SasView

.gitignore

r67cc0ff	r6ceca44
24	24	/example/Fit_*/
25	25	/example/batch_fit.csv
	26	# Note: changes to gauss20, gauss76 and gauss150 are still tracked since they
	27	# are part of the repo and required by some models.
26	28	/sasmodels/models/lib/gauss*.c

.travis.yml

-                      rce8c388
+                      r2d09df1
     env:
     - PY=2.7
+    - DEPLOY=True
   - os: linux
     env:
 …
   on:
     branch: master
+    condition: $DEPLOY = True
 notifications:
   slack:

appveyor.yml

rd810d96	r1258e32
67	67	- cmd: conda install --yes --quiet obvious-ci
68	68	- cmd: conda install --yes --quiet numpy toolchain scipy cython cffi
69		- cmd: conda install --yes --channel conda-forge pyopencl
	69	#- cmd: conda install --yes --channel conda-forge pyopencl
70	70	- cmd: pip install bumps unittest-xml-reporting tinycc
71	71

doc/developer/overview.rst

-                      r3d40839
+                      r2a7e20e
 jitter applied before particle orientation.
+When computing the orientation dispersity integral, the weights for
+the individual points depends on the map projection used to translate jitter
+angles into latitude/longitude.  The choice of projection is set by
+*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for
+equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated
+Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh
+for details.
 For numerical integration within form factors etc. sasmodels is mostly using
 Gaussian quadrature with 20, 76 or 150 points depending on the model. It also
 …
 Useful testing routines include:
+:mod:`asymint` a direct implementation of the surface integral for certain
+models to get a more trusted value for the 1D integral using a
+reimplementation of the 2D kernel in python and mpmath (which computes math
+functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$
+and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including
+the U-substitution for $\theta$.
+:mod:`check1d` uses sasmodels 1D integration and compares that with a
+The *sascomp* utility is used to view and compare models with different
+parameters and calculation engines. The usual case is to simply plot a
+model that you are developing::
+    python sascomp path/to/model.py
+Once the obvious problems are addressed, check the numerical precision
+across a variety of randomly generated inputs::
+    python sascomp -engine=single,double path/to/model.py -sets=10
+You can compare different parameter values for the same or different models.
+For example when looking along the long axis of a cylinder ($\theta=0$),
+dispersity in $\theta$ should be equivalent to dispersity in $\phi$
+when $\phi=90$::
+    python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\
+    phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10
+It turns out that they are not because the equirectangular map projection
+weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical
+to $\Delta\phi$.  Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps
+somewhat.  Postel would help even more in this case, though leading
+to distortions elsewhere.  See :mod:`sasmodels.compare` for many more details.
+*sascomp -ngauss=n* allows you to set the number of quadrature points used
+for the 1D integration for any model.  For example, a carbon nanotube with
+length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$::
+    python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq
+Note: ticket 702 gives some forms for long rods and thin disks that may
+be used in place of the integration, depending on $q$, radius and length; if
+the cylinder model is updated with these corrections then above call will show
+no difference.
+*explore/check1d.py* uses sasmodels 1D integration and compares that with a
 rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to
 $\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle
 …
 similar reasoning.] This should rotate the sample through the entire
 $\theta$-$\phi$ surface according to the pattern that you see in jitter.py when
 you modify it to use 'rectangle' rather than 'gaussian' for its distribution
+without changing the viewing angle. In order to match the 1-D pattern for
 an arbitrary viewing angle on triaxial shapes, we need to integrate
+you use 'rectangle' rather than 'gaussian' for its distribution without
+changing the viewing angle. In order to match the 1-D pattern for an arbitrary
+viewing angle on triaxial shapes, we need to integrate
 over $\theta$, $\phi$ and $\Psi$.
+When computing the dispersity integral, weights are scaled by
+$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer
+together as $\delta \theta$ increases.
+[This will probably change so that instead of adjusting the weights, we will
+adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in
+$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in
+*kernel_iq.c* selects an alternative algorithm.]
+The integrated dispersion is computed at a set of $(qx, qy)$ points $(q
+\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for
+each $q$ used in the 1-D integration. The individual $q$ points should be
+equivalent to asymint rect-n when the viewing angle is set to
+$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d
+models are consistent with 1d models.
+:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to
+compute the pattern of the $q_x$-$q_y$ grid.
+The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to
+compare results for two sets of parameters or processor types, for example
+these two oriented cylinders here should be equivalent.
+:mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10`
+*sascomp -sphere=n* uses the same rectangular distribution as check1d to
+compute the pattern of the $q_x$-$q_y$ grid.  This ought to produce a
+spherically symmetric pattern.
+*explore/precision.py* investigates the accuracy of individual functions
+on the different execution platforms.  This includes the basic special
+functions as well as more complex expressions used in scattering.  In many
+cases the OpenCL function in sasmodels will use a polynomial approximation
+over part of the range to improve accuracy, as shown in::
+    python explore/precision.py 3j1/x
+*explore/realspace.py* allows you to set up a Monte Carlo simulation of your
+model by sampling random points within and computing the 1D and 2D scattering
+patterns.  This was used to check the core shell parallelepiped example.  This
+is similar to the general sas calculator in sasview, though it uses different
+code.
+*explore/jitter.py* is for exploring different options for handling
+orientation and orientation dispersity.  It uses *explore/guyou.py* to
+generate the Guyou projection.
+*explore/asymint.py* is a direct implementation of the 1D integration for
+a number of different models.  It calculates the integral for a particular
+$q$ using several different integration schemes, including mpmath with 100
+bits of precision (33 digits), so you can use it to check the target values
+for the 1D model tests.  This is not a general purpose tool; you will need to
+edit the file to change the model parameters. It does not currently
+apply the $u=cos(\theta)$ substitution that many models are using
+internally so the 76-point gaussian quadrature may not match the value
+produced by the eqivalent model from sasmodels.
+*explore/symint.py* is for rotationally symmetric models (just cylinder for
+now), so it can compute an entire curve rather than a single $q$ point.  It
+includes code to compare the long cylinder approximation to cylinder.
+*explore/rpa.py* is for checking different implementations of the RPA model
+against calculations for specific blends.  This is a work in (suspended)
+progress.
+There are a few modules left over in *explore* that can probably be removed.
+*angular_pd.py* was an early version of *jitter.py*.  *sc.py* and *sc.c*
+was used to explore different calculations for paracrystal models; it has
+been absorbed into *asymint.py*. *transform_angles.py* translates orientation
+parameters from the SasView 3.x definition to sasmodels.
+*explore/angles.py* generates code for the view and jitter transformations.
+Keep this around since it may be needed if we add new projections.
 Testing

explore/precision.py

-                      ra1c5758
+                      r2a7e20e
+)
 add_function(
     name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x",
+    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x",
     mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x,
     np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x,
 …
     names = ", ".join(sorted(ALL_FUNCTIONS))
     print("""\
 usage: precision.py [-f/a/r] [-x<range>] name...
+usage: precision.py [-f/a/r] [-x<range>] "name" ...
 where
     -f indicates that the function value should be plotted,
 …
       zoom indicates linear stepping in [1000, 1010]
       neg indicates linear stepping in [-100.1, 100.1]
 and name is "all [first]" or one of:
+and name is "all" or one of:
     """+names)
     sys.exit(1)

explore/realspace.py

-                      ra1c32c2
+                      r8fb2a94
     """
     return Rx(phi)*Ry(theta)*Rz(psi)
+def apply_view(points, view):
+    """
+    Apply the view transform (theta, phi, psi) to a set of points.
+    Points are stored in a 3 x n numpy array.
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T
+def invert_view(qx, qy, view):
+    """
+    Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector
+    pixel (qx, qy).
+    View angles are in degrees.
+    """
+    theta, phi, psi = view
+    q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten()))
+    return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q))
 class Shape:
 …
         values = self.value.repeat(points.shape[0])
         return values, self._adjust(points)
+NUMBA = False
+if NUMBA:
+    from numba import njit
+    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])")
+    def _Iqxy(values, x, y, z, qa, qb, qc):
+        Iq = np.zeros_like(qa)
+        for j in range(len(Iq)):
+            total = 0. + 0j
+            for k in range(len(Iq)):
+                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k]))
+            Iq[j] = abs(total)**2
+        return Iq
+def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)):
+    qx, qy = np.broadcast_arrays(qx, qy)
+    qa, qb, qc = invert_view(qx, qy, view)
+    rho, volume = np.broadcast_arrays(rho, volume)
+    values = rho*volume
+    x, y, z = points.T
+    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)
+    if NUMBA:
+        Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten())
+    else:
+        Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2
+              for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)]
+    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume)
 def _calc_Pr_nonuniform(r, rho, points):
 …
             print("processing %d of %d"%(k, len(rho)-1))
     Pr = extended_Pr[1:-1]
+    return Pr / Pr.max()
+def calc_Pr(r, rho, points):
+    # P(r) with uniform steps in r is 3x faster; check if we are uniform
+    # before continuing
+    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
+        return _calc_Pr_nonuniform(r, rho, points)
+    return Pr
+def _calc_Pr_uniform(r, rho, points):
     # Make Pr a little be bigger than necessary so that only distances
     # min < d < max end up in Pr
     n_max = len(r)
+    dr, n_max = r[0], len(r)
     extended_Pr = np.zeros(n_max+1, 'd')
     t0 = time.clock()
     t_next = t0 + 3
-    row_zero = np.zeros(len(rho), 'i')
     for k, rho_k in enumerate(rho[:-1]):
         distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1))
         weights = rho_k * rho[k+1:]
         index = np.minimum(np.asarray(distances/r[0], 'i'), n_max)
+        index = np.minimum(np.asarray(distances/dr, 'i'), n_max)
         # Note: indices may be duplicated, so "Pr[index] += w" will not work!!
         extended_Pr += np.bincount(index, weights, n_max+1)
 …
     # we are only accumulating the upper triangular distances.
     #Pr = Pr * 2 / len(rho)**2
     return Pr / Pr.max()
+    return Pr
     # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even
 …
+}
 """
+if NUMBA:
+    @njit("f8[:](f8[:], f8[:], f8[:,:])")
+    def _calc_Pr_uniform_jit(r, rho, points):
+        dr = r[0]
+        n_max = len(r)
+        Pr = np.zeros_like(r)
+        for j in range(len(rho) - 1):
+            x, y, z = points[j, 0], points[j, 1], points[j, 2]
+            for k in range(j+1, len(rho)):
+                distance = sqrt((x - points[k, 0])**2
+                                + (y - points[k, 1])**2
+                                + (z - points[k, 2])**2)
+                index = int(distance/dr)
+                if index < n_max:
+                    Pr[index] += rho[j] * rho[k]
+        # Make Pr independent of sampling density.  The factor of 2 comes because
+        # we are only accumulating the upper triangular distances.
+        #Pr = Pr * 2 / len(rho)**2
+        return Pr
+def calc_Pr(r, rho, points):
+    # P(r) with uniform steps in r is 3x faster; check if we are uniform
+    # before continuing
+    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01:
+        Pr = _calc_Pr_nonuniform(r, rho, points)
+    else:
+        if NUMBA:
+            Pr = _calc_Pr_uniform_jit(r, rho, points)
+        else:
+            Pr = _calc_Pr_uniform(r, rho, points)
+    return Pr / Pr.max()
 def j0(x):
 …
         plt.legend()
+def plot_calc_2d(qx, qy, Iqxy, theory=None):
+    import matplotlib.pyplot as plt
+    qx, qy = bin_edges(qx), bin_edges(qy)
+    #qx, qy = np.meshgrid(qx, qy)
+    if theory is not None:
+        plt.subplot(121)
+    plt.pcolormesh(qx, qy, np.log10(Iqxy))
+    plt.xlabel('qx (1/A)')
+    plt.ylabel('qy (1/A)')
+    if theory is not None:
+        plt.subplot(122)
+        plt.pcolormesh(qx, qy, np.log10(theory))
+        plt.xlabel('qx (1/A)')
 def plot_points(rho, points):
     import mpl_toolkits.mplot3d
 …
         pass
     n = len(points)
+    index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None)
+    #print("len points", n)
+    index = np.random.choice(n, size=500) if n > 500 else slice(None, None)
     ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index])
     #low, high = points.min(axis=0), points.max(axis=0)
     #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]])
     ax.autoscale(True)
+def check_shape(shape, fn=None):
+    rho_solvent = 0
+    q = np.logspace(-3, 0, 200)
+    r = shape.r_bins(q, r_step=0.01)
+    sampling_density = 6*5000 / shape.volume()
+    rho, points = shape.sample(sampling_density)
+    t0 = time.time()
+    Pr = calc_Pr(r, rho-rho_solvent, points)
+    print("calc Pr time", time.time() - t0)
+    Iq = calc_Iq(q, r, Pr)
+    theory = (q, fn(q)) if fn is not None else None
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc(r, Pr, q, Iq, theory=theory)
+    pylab.show()
+def check_shape_2d(shape, fn=None, view=(0, 0, 0)):
+    rho_solvent = 0
+    nq, qmax = 100, 1.0
+    qx = np.linspace(0.0, qmax, nq)
+    qy = np.linspace(0.0, qmax, nq)
+    Qx, Qy = np.meshgrid(qx, qy)
+    sampling_density = 50000 / shape.volume()
+    #t0 = time.time()
+    rho, points = shape.sample(sampling_density)
+    #print("sample time", time.time() - t0)
+    t0 = time.time()
+    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view)
+    print("calc time", time.time() - t0)
+    theory = fn(Qx, Qy) if fn is not None else None
+    Iqxy += 0.001 * Iqxy.max()
+    if theory is not None:
+        theory += 0.001 * theory.max()
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc_2d(qx, qy, Iqxy, theory=theory)
+    pylab.show()
+def sas_sinx_x(x):
+    with np.errstate(all='ignore'):
+        retvalue = sin(x)/x
+    retvalue[x == 0.] = 1.
+    return retvalue
 def sas_2J1x_x(x):
 …
     Iq = Iq/Iq[0]
     return Iq
+def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+    qab = np.sqrt(qa**2 + qb**2)
+    Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2)
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 def sphere_Iq(q, radius):
 …
     return Iq/Iq[0]
+def check_shape(shape, fn=None):
+    rho_solvent = 0
+    q = np.logspace(-3, 0, 200)
+    r = shape.r_bins(q, r_step=0.01)
+    sampling_density = 15000 / shape.volume()
+    rho, points = shape.sample(sampling_density)
+    Pr = calc_Pr(r, rho-rho_solvent, points)
+    Iq = calc_Iq(q, r, Pr)
+    theory = (q, fn(q)) if fn is not None else None
+    import pylab
+    #plot_points(rho, points); pylab.figure()
+    plot_calc(r, Pr, q, Iq, theory=theory)
+    pylab.show()
+def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)):
+    qa, qb, qc = invert_view(qx, qy, view)
+    sld_solvent = 0
+    overlapping = False
+    dr0 = sld_core - sld_solvent
+    drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent
+    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
+    siA = a*sas_sinx_x(a*qa/2)
+    siB = b*sas_sinx_x(b*qb/2)
+    siC = c*sas_sinx_x(c*qc/2)
+    siAt = tA*sas_sinx_x(tA*qa/2)
+    siBt = tB*sas_sinx_x(tB*qb/2)
+    siCt = tC*sas_sinx_x(tC*qc/2)
+    Fq = (dr0*siA*siB*siC
+          + drA*(siAt-siA)*siB*siC
+          + drB*siA*(siBt-siB)*siC
+          + drC*siA*siB*(siCt-siC))
+    Iq = Fq**2
+    return Iq.reshape(qx.shape)
 def check_cylinder(radius=25, length=125, rho=2.):
 …
     fn = lambda q: cylinder_Iq(q, radius, length)
     check_shape(shape, fn)
+def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)):
+    shape = EllipticalCylinder(radius, radius, length, rho)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
+def check_cylinder_2d_lattice(radius=25, length=125, rho=2.,
+                              view=(0, 0, 0)):
+    nx, dx = 1, 2*radius
+    ny, dy = 30, 2*radius
+    nz, dz = 30, length
+    dx, dy, dz = 2*dx, 2*dy, 2*dz
+    def center(*args):
+        sigma = 0.333
+        space = 2
+        return [(space*n+np.random.randn()*sigma)*x for n, x in args]
+    shapes = [EllipticalCylinder(radius, radius, length, rho,
+                                 #center=(ix*dx, iy*dy, iz*dz)
+                                 orientation=np.random.randn(3)*0,
+                                 center=center((ix, dx), (iy, dy), (iz, dz))
+                                )
+              for ix in range(nx)
+              for iy in range(ny)
+              for iz in range(nz)]
+    shape = Composite(shapes)
+    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view)
+    check_shape_2d(shape, fn, view=view)
 def check_sphere(radius=125, rho=2):
 …
     side_c2 = copy(side_c).shift(0, 0, -c-dc)
     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2))
+    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
+    check_shape(shape, fn)
+    def fn(q):
+        return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core)
+    #check_shape(shape, fn)
+    view = (20, 30, 40)
+    def fn_xy(qx, qy):
+        return csbox_Iqxy(qx, qy, a, b, c, da, db, dc,
+                          slda, sldb, sldc, sld_core, view=view)
+    check_shape_2d(shape, fn_xy, view=view)
 if __name__ == "__main__":
     check_cylinder(radius=10, length=40)
+    #check_cylinder_2d(radius=10, length=40, view=(90,30,0))
+    #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0))
     #check_sphere()
     #check_csbox()

sasmodels/compare.py

ra261a83	r2a7e20e
92	92	-accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
93	93	-neval=1 sets the number of evals for more accurate timing
	94	-ngauss=0 overrides the number of points in the 1-D gaussian quadrature
94	95
95	96	=== precision options ===

sasmodels/kernel_iq.c

-                      r108e70e
+                      r924a119
 static double
 qac_apply(
     QACRotation rotation,
+    QACRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qc_out)
+{
     const double dqc = rotation.R31*qx + rotation.R32*qy;
+    const double dqc = rotation->R31*qx + rotation->R32*qy;
     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2
     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy);
 …
 static double
 qabc_apply(
     QABCRotation rotation,
+    QABCRotation *rotation,
     double qx, double qy,
     double *qa_out, double *qb_out, double *qc_out)
+{
     *qa_out = rotation.R11*qx + rotation.R12*qy;
     *qb_out = rotation.R21*qx + rotation.R22*qy;
     *qc_out = rotation.R31*qx + rotation.R32*qy;
+    *qa_out = rotation->R11*qx + rotation->R12*qy;
+    *qb_out = rotation->R21*qx + rotation->R22*qy;
+    *qc_out = rotation->R31*qx + rotation->R32*qy;
+}
 …
   // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code.
   #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi);
   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc)
+  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc)
   #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table)
 …
   local_values.table.psi = 0.;
   #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
+  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
 #elif defined(CALL_IQ_XY)
 …
 #endif
 // Doing jitter projection code outside the previous if block so that we don't
 // need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This
 // will become more important if we implement more projections, or more
 // complicated projections.
 #if defined(CALL_IQ) || defined(CALL_IQ_A)
+// Define APPLY_PROJECTION depending on model symmetries. We do this outside
+// the previous if block so that we don't need to repeat the identical
+// logic in the IQ_AC and IQ_ABC branches.  This will become more important
+// if we implement more projections, or more complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
 #elif defined(CALL_IQ_XY)
+#elif defined(CALL_IQ_XY) // pass orientation to the model
   // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc
   // Need to plug the values for the orientation angles back into parameter
 …
     #define APPLY_PROJECTION() const double weight=weight0
   #endif
 #else // !spherosymmetric projection
+#else // apply jitter and view before calling the model
   // Grab the "view" angles (theta, phi, psi) from the initial parameter table.
   const double theta = values[details->theta_par+2];
 …
   // we go through the mesh.
   double dtheta, dphi, weight;
   #if PROJECTION == 1
+  #if PROJECTION == 1 // equirectangular
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
 …
       weight = fabs(cos(dtheta*M_PI_180)) * weight0; \
     } while (0)
   #elif PROJECTION == 2
+  #elif PROJECTION == 2 // sinusoidal
     #define APPLY_PROJECTION() do { \
       dtheta = local_values.table.theta; \
 …
     } while (0)
   #endif
 #endif // !spherosymmetric projection
+#endif // done defining APPLY_PROJECTION
 // ** define looping macros **

sasmodels/models/core_shell_parallelepiped.c

-                      r108e70e
+                      re077231
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
     // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
     // Code rewritten (PAK)
+    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
+    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
     const double half_q = 0.5*q;
 …
     const double drC = crim_sld-solvent_sld;
-    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B.
     const double tA = length_a + 2.0*thick_rim_a;
     const double tB = length_b + 2.0*thick_rim_b;

sasmodels/models/core_shell_parallelepiped.py

-                      r10ee838
+                      r97be877
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
+* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
+* **Currently Under review by:** Paul Butler
+* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017
+* Cross-checked against hollow rectangular prism and rectangular prism for
+  equal thickness overlapping sides, and by Monte Carlo sampling of points
+  within the shape for non-uniform, non-overlapping sides.
 """

sasmodels/models/lib/gauss150.c

-                      r74768cb
+                      r99b84ec
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+// Created by Andrew Jackson on 4/23/07
  #ifdef GAUSS_N
  # undef GAUSS_N
 …
  #define GAUSS_W Gauss150Wt
-// Note: using array size 152 so that it is a multiple of 4
+// Gaussians
+// Note: using array size 152 rather than 150 so that it is a multiple of 4.
+// Some OpenCL devices prefer that vectors start and end on nice boundaries.
 constant double Gauss150Z[152]={
         -0.9998723404457334,
 …
 .9993274305065947,
 .9998723404457334,
 .,
 .
+., // zero padding is ignored
+.  // zero padding is ignored
 };
 …
 .0007624720924706,
 .0003276086705538,
 .,
 .
+., // zero padding is ignored
+.  // zero padding is ignored
 };

sasmodels/models/lib/gauss20.c

-                      r74768cb
+                      r99b84ec
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+// Created by Andrew Jackson on 4/23/07
  #ifdef GAUSS_N
  # undef GAUSS_N

sasmodels/models/lib/gauss76.c

-                      r74768cb
+                      r99b84ec
+/*
+ *  GaussWeights.c
+ *  SANSAnalysis
+ *
+ *  Created by Andrew Jackson on 4/23/07.
+ *  Copyright 2007 __MyCompanyName__. All rights reserved.
+ *
+ */
+// Created by Andrew Jackson on 4/23/07
  #ifdef GAUSS_N
  # undef GAUSS_N

sasmodels/kerneldll.py

-                      r2d81cfe
+                      r1ddb794
         subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT)
     except subprocess.CalledProcessError as exc:
+        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output))
+        raise RuntimeError("compile failed.\n%s\n%s"
+                           % (command_str, exc.output.decode()))
     if not os.path.exists(output):
         raise RuntimeError("compile failed.  File is in %r"%source)

sasmodels/modelinfo.py

-                      r108e70e
+                      r5ab99b7
 from os.path import abspath, basename, splitext
 import inspect
+import logging
 import numpy as np  # type: ignore
+from . import autoc
 # Optional typing
 …
     TestCondition = Tuple[ParameterSetUser, TestInput, TestValue]
 # pylint: enable=unused-import
+logger = logging.getLogger(__name__)
 # If MAX_PD changes, need to change the loop macros in kernel_iq.c
 …
     info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore
     # Default single and opencl to True for C models.  Python models have callable Iq.
-    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
-    info.single = getattr(kernel_module, 'single', not callable(info.Iq))
     info.random = getattr(kernel_module, 'random', None)
 …
     info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore
+    info.lineno = {}
+    _find_source_lines(info, kernel_module)
+    if getattr(kernel_module, 'py2c', False):
+        try:
+            autoc.convert(info, kernel_module)
+        except Exception as exc:
+            logger.warn(str(exc) + " while converting %s from C to python"%name)
+    # Needs to come after autoc.convert since the Iq symbol may have been
+    # converted from python to C
+    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
+    info.single = getattr(kernel_module, 'single', not callable(info.Iq))
     if callable(info.Iq) and parameters.has_2d:
         raise ValueError("oriented python models not supported")
-    info.lineno = {}
-    _find_source_lines(info, kernel_module)
     return info

SasView

Changeset 6ceca44 in sasmodels

Legend:

.gitignore

.travis.yml

appveyor.yml

doc/developer/overview.rst

explore/precision.py

explore/realspace.py

sasmodels/compare.py

sasmodels/kernel_iq.c

sasmodels/models/core_shell_parallelepiped.c

sasmodels/models/core_shell_parallelepiped.py

sasmodels/models/lib/gauss150.c

sasmodels/models/lib/gauss20.c

sasmodels/models/lib/gauss76.c

sasmodels/kerneldll.py

sasmodels/modelinfo.py

Download in other formats: