Changeset 2a7e20e in sasmodels

Timestamp:

Jan 11, 2018 11:02:39 AM (6 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:

e077231

Parents:

2ab331f

Message:

update developer docs with current interpretation of orientation; describe the scripts in the explore directory

Files:

• doc/developer/overview.rst (modified) (3 diffs)
• explore/precision.py (modified) (3 diffs)
• sasmodels/compare.py (modified) (1 diff)
• sasmodels/kernel_iq.c (modified) (5 diffs)

doc/developer/overview.rst

@@ -185 +185 @@
 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
@@ -199 +207 @@
 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
@@ -214 +251 @@
 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

@@ -345 +345 @@
 )
 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,
@@ -609 +609 @@
     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,
@@ -620 +620 @@
       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)

sasmodels/compare.py

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

sasmodels/kernel_iq.c

@@ -479 +479 @@
 #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
@@ -515 +515 @@
     #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];
@@ -526 +526 @@
   // 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; \
@@ -532 +532 @@
       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; \
@@ -542 +542 @@
     } while (0)
   #endif
-#endif // !spherosymmetric projection
+#endif // done defining APPLY_PROJECTION
 
 // ** define looping macros **