Changeset a6d3f46 in sasmodels

Timestamp:

Jan 3, 2018 9:42:20 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Children:

5ab99b7

Parents:

33756c8 (diff), 5fb0634 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge remote-tracking branch 'upstream/ticket-1043'

Files:

• doc/guide/orientation/orientation.rst (modified) (3 diffs)
• doc/guide/plugin.rst (modified) (1 diff)
• explore/jitter.py (modified) (1 diff)
• .gitignore (modified) (2 diffs)
• sasmodels/autoc.py (added)
• sasmodels/kerneldll.py (modified) (1 diff)
• sasmodels/modelinfo.py (modified) (5 diffs)
• sasmodels/models/_cylpy.py (added)
• sasmodels/py2c.py (added)

doc/guide/orientation/orientation.rst

@@ -4 +4 @@
 ==================
 
-With two dimensional small angle diffraction data SasView will calculate
+With two dimensional small angle diffraction data sasmodels will calculate
 scattering from oriented particles, applicable for example to shear flow
 or orientation in a magnetic field.
 
 In general we first need to define the reference orientation
-of the particles with respect to the incoming neutron or X-ray beam. This
-is done using three angles: $\theta$ and $\phi$ define the orientation of
-the axis of the particle, angle $\Psi$ is defined as the orientation of
-the major axis of the particle cross section with respect to its starting
-position along the beam direction. The figures below are for an elliptical
-cross section cylinder, but may be applied analogously to other shapes of
-particle.
+of the particle's $a$-$b$-$c$ axes with respect to the incoming
+neutron or X-ray beam. This is done using three angles: $\theta$ and $\phi$
+define the orientation of the $c$-axis of the particle, and angle $\Psi$ is
+defined as the orientation of the major axis of the particle cross section
+with respect to its starting position along the beam direction (or
+equivalently, as rotation about the $c$ axis). There is an unavoidable
+ambiguity when $c$ is aligned with $z$ in that $\phi$ and $\Psi$ both
+serve to rotate the particle about $c$, but this symmetry is destroyed
+when $\theta$ is not a multiple of 180.
+
+The figures below are for an elliptical cross section cylinder, but may
+be applied analogously to other shapes of particle.
 
 .. note::
@@ -29 +34 @@
 
     Definition of angles for oriented elliptical cylinder, where axis_ratio
-    b/a is shown >1, Note that rotation $\theta$, initially in the $x$-$z$
+    b/a is shown >1. Note that rotation $\theta$, initially in the $x$-$z$
     plane, is carried out first, then rotation $\phi$ about the $z$-axis,
     finally rotation $\Psi$ is around the axis of the cylinder. The neutron
-    or X-ray beam is along the $z$ axis.
+    or X-ray beam is along the $-z$ axis.
 
 .. figure::
@@ -40 +45 @@
     with $\Psi$ = 0.
 
-Having established the mean direction of the particle we can then apply
-angular orientation distributions. This is done by a numerical integration
-over a range of angles in a similar way to particle size dispersity.
-In the current version of sasview the orientational dispersity is defined
-with respect to the axes of the particle.
+Having established the mean direction of the particle (the view) we can then
+apply angular orientation distributions (jitter). This is done by a numerical
+integration over a range of angles in a similar way to particle size
+dispersity. The orientation dispersity is defined with respect to the
+$a$-$b$-$c$ axes of the particle, with roll angle $\Psi$ about the $c$-axis,
+yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the
+$a$-axis.
+
+More formally, starting with axes $a$-$b$-$c$ of the particle aligned
+with axes $x$-$y$-$z$ of the laboratory frame, the orientation dispersity
+is applied first, using the
+`Tait-Bryan <https://en.wikipedia.org/wiki/Euler_angles#Conventions_2>`_
+$x$-$y'$-$z''$ convention with angles $\Delta\phi$-$\Delta\theta$-$\Delta\Psi$.
+The reference orientation then follows, using the
+`Euler angles <https://en.wikipedia.org/wiki/Euler_angles#Conventions>`_
+$z$-$y'$-$z''$ with angles $\phi$-$\theta$-$\Psi$.  This is implemented
+using rotation matrices as
+
+.. math::
+
+    R = R_z(\phi)\, R_y(\theta)\, R_z(\Psi)\,
+        R_x(\Delta\phi)\, R_y(\Delta\theta)\, R_z(\Delta\Psi)
+
+To transform detector $(q_x, q_y)$ values into $(q_a, q_b, q_c)$ for the
+shape in its canonical orientation, use
+
+.. math::
+
+    [q_a, q_b, q_c]^T = R^{-1} \, [q_x, q_y, 0]^T
+
+
+The inverse rotation is easily calculated by rotating the opposite directions
+in the reverse order, so
+
+.. math::
+
+    R^{-1} = R_z(-\Delta\Psi)\, R_y(-\Delta\theta)\, R_x(-\Delta\phi)\,
+             R_z(-\Psi)\, R_y(-\theta)\, R_z(-\phi)
+
 
 The $\theta$ and $\phi$ orientation parameters for the cylinder only appear
-when fitting 2d data. On introducing "Orientational Distribution" in
-the angles, "distribution of theta" and "distribution of phi" parameters will
+when fitting 2d data. On introducing "Orientational Distribution" in the
+angles, "distribution of theta" and "distribution of phi" parameters will
 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$
-of the cylinder, the $b$ and $a$ axes of the cylinder cross section. (When
-$\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the
-instrument.) The third orientation distribution, in $\Psi$, is about the $c$
-axis of the particle. Some experimentation may be required to understand the
-2d patterns fully. A number of different shapes of distribution are
-available, as described for polydispersity, see :ref:`polydispersityhelp` .
+of the cylinder, which correspond to the $b$ and $a$ axes of the cylinder
+cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and
+$X$ axes of the instrument.) The third orientation distribution, in $\Psi$,
+is about the $c$ axis of the particle. Some experimentation may be required
+to understand the 2d patterns fully. A number of different shapes of
+distribution are available, as described for size dispersity, see
+:ref:`polydispersityhelp`.
 
-Earlier versions of SasView had numerical integration issues in some
-circumstances when distributions passed through 90 degrees. The distributions
-in particle coordinates are more robust, but should still be approached with
-care for large ranges of angle.
+Given that the angular dispersion distribution is defined in cartesian space,
+over a cube defined by
+
+.. math::
+
+    [-\Delta \theta, \Delta \theta] \times
+    [-\Delta \phi, \Delta \phi] \times
+    [-\Delta \Psi, \Delta \Psi]
+
+but the orientation is defined over a sphere, we are left with a
+`map projection <https://en.wikipedia.org/wiki/List_of_map_projections>`_
+problem, with different tradeoffs depending on how values in $\Delta\theta$
+and $\Delta\phi$ are translated into latitude/longitude on the sphere.
+
+Sasmodels is using the
+`equirectangular projection <https://en.wikipedia.org/wiki/Equirectangular_projection>`_.
+In this projection, square patches in angular dispersity become wedge-shaped
+patches on the sphere. To correct for the changing point density, there is a
+scale factor of $\sin(\Delta\theta)$ that applies to each point in the
+integral. This is not enough, though. Consider a shape which is tumbling
+freely around the $b$ axis, with $\Delta\theta$ uniform in $[-180, 180]$. At
+$\pm 90$, all points in $\Delta\phi$ map to the pole, so the jitter will have
+a distinct angular preference. If the spin axis is along the beam (which
+will be the case for $\theta=90$ and $\Psi=90$) the scattering pattern
+should be circularly symmetric, but it will go to zero at $q_x = 0$ due to the
+$\sin(\Delta\theta)$ correction. This problem does not appear for a shape
+that is tumbling freely around the $a$ axis, with $\Delta\phi$ uniform in
+$[-180, 180]$, so swap the $a$ and $b$ axes so $\Delta\theta < \Delta\phi$
+and adjust $\Psi$ by 90. This works with the current sasmodels shapes due to
+symmetry.
+
+Alternative projections were considered.
+The `sinusoidal projection <https://en.wikipedia.org/wiki/Sinusoidal_projection>`_
+works by scaling $\Delta\phi$ as $\Delta\theta$ increases, and dropping those
+points outside $[-180, 180]$. The distortions are a little less for middle
+ranges of $\Delta\theta$, but they are still severe for large $\Delta\theta$
+and the model is much harder to explain.
+The `azimuthal equidistance projection <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>`_
+also improves on the equirectangular projection by extending the range of
+reasonable values for the $\Delta\theta$ range, with $\Delta\phi$ forming a
+wedge that cuts to the opposite side of the sphere rather than cutting to the
+pole. This projection has the nice property that distance from the center are
+preserved, and that $\Delta\theta$ and $\Delta\phi$ act the same.
+The `azimuthal equal area projection <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>`_
+is like the azimuthal equidistance projection, but it preserves area instead
+of distance. It also has the same behaviour for $\Delta\theta$ and $\Delta\phi$.
+The `Guyou projection <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>`_
+has an excellent balance with reasonable distortion in both $\Delta\theta$
+and $\Delta\phi$, as well as preserving small patches. However, it requires
+considerably more computational overhead, and we have not yet derived the
+formula for the distortion correction, measuring the degree of stretch at
+the point $(\Delta\theta, \Delta\phi)$ on the map.
 
 .. note::
-    Note that the form factors for oriented particles are also performing
-    numerical integrations over one or more variables, so care should be taken,
-    especially with very large particles or more extreme aspect ratios. In such 
-    cases results may not be accurate, particularly at very high Q, unless the model
-    has been specifically coded to use limiting forms of the scattering equations.
-    
-    For best numerical results keep the $\theta$ distribution narrower than the $\phi$ 
-    distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may 
-    need to reorder the sizes of the three axes to acheive the desired result. 
-    This is due to the issues of mapping a rectangular distribution onto the 
-    surface of a sphere.
+    Note that the form factors for oriented particles are performing
+    numerical integrations over one or more variables, so care should be
+    taken, especially with very large particles or more extreme aspect
+    ratios. In such cases results may not be accurate, particularly at very
+    high Q, unless the model has been specifically coded to use limiting
+    forms of the scattering equations.
 
-Users can experiment with the values of *Npts* and *Nsigs*, the number of steps 
-used in the integration and the range spanned in number of standard deviations. 
-The standard deviation is entered in units of degrees. For a "rectangular" 
-distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations. 
-The new "uniform" distribution avoids this by letting you directly specify the 
+    For best numerical results keep the $\theta$ distribution narrower than
+    the $\phi$ distribution. Thus for asymmetric particles, such as
+    elliptical_cylinder, you may need to reorder the sizes of the three axes
+    to acheive the desired result. This is due to the issues of mapping a
+    rectanglar distribution onto the surface of a sphere.
+
+Users can experiment with the values of *Npts* and *Nsigs*, the number of steps
+used in the integration and the range spanned in number of standard deviations.
+The standard deviation is entered in units of degrees. For a "rectangular"
+distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.
+The new "uniform" distribution avoids this by letting you directly specify the
 half width.
 
-The angular distributions will be truncated outside of the range -180 to +180 
-degrees, so beware of using saying a broad Gaussian distribution with large value
-of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would 
-give a good integration,as well as becoming rather meaningless. (At some point 
-in the future the actual polydispersity arrays may be made available to the user 
-for inspection.)
+The angular distributions may be truncated outside of the range -180 to +180
+degrees, so beware of using saying a broad Gaussian distribution with large
+value of *Nsigs*, as the array of *Npts* may be truncated to many fewer
+points than would give a good integration,as well as becoming rather
+meaningless. (At some point in the future the actual dispersion arrays may be
+made available to the user for inspection.)
 
 Some more detailed technical notes are provided in the developer section of
 this manual :ref:`orientation_developer` .
 
+This definition of orientation is new to SasView 4.2.  In earlier versions,
+the orientation distribution appeared as a distribution of view angles.
+This led to strange effects when $c$ was aligned with $z$, where changes
+to the $\phi$ angle served only to rotate the shape about $c$, rather than
+having a consistent interpretation as the pitch of the shape relative to
+the flow field defining the reference orientation.  Prior to SasView 4.1,
+the reference orientation was defined using a Tait-Bryan convention, making
+it difficult to control.  Now, rotation in $\theta$ modifies the spacings
+in the refraction pattern, and rotation in $\phi$ rotates it in the detector
+plane.
+
+
 *Document History*
 
 | 2017-11-06 Richard Heenan
+| 2017-12-20 Paul Kienzle

doc/guide/plugin.rst

@@ -538 +538 @@
 If the scattering is dependent on the orientation of the shape, then you
 will need to include *orientation* parameters *theta*, *phi* and *psi*
-at the end of the parameter table.  Shape orientation uses *a*, *b* and *c*
-axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate
-system, with *z* along the beam and *x*-*y* in the detector plane, with *x*
-horizontal and *y* vertical.  The *psi* parameter rotates the shape
-about its *c* axis, the *theta* parameter then rotates the *c* axis toward
-the *x* axis of the detector, then *phi* rotates the shape in the detector
-plane.  (Prior to these rotations, orientation dispersity will be applied
-as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate
-system.)  A particular *qx*, *qy* point on the detector, then corresponds
-to *qa*, *qb*, *qc* with respect to the shape.
-
-The oriented C model is called as *Iqabc(qa, qb, qc, par1, par2, ...)* where
+at the end of the parameter table.  As described in the section
+:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will
+be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its
+canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the
+laboratory frame and beam travelling along $-z$.
+
+The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
 *par1*, etc. are the parameters to the model.  If the shape is rotationally
 symmetric about *c* then *psi* is not needed, and the model is called

explore/jitter.py

@@ -165 +165 @@
     # constants in kernel_iq.c
     'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
+    'azimuthal_equal_area',
 ]
 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',

.gitignore

@@ -8 +8 @@
 *.so
 *.obj
+*.o
 /doc/_build/
 /doc/api/
@@ -19 +20 @@
 /.pydevproject
 /.idea
+.vscode
 /sasmodels.egg-info/
 /example/Fit_*/

sasmodels/kerneldll.py

@@ -185 +185 @@
         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

@@ -12 +12 @@
 from os.path import abspath, basename, splitext
 import inspect
+import logging
 
 import numpy as np  # type: ignore
+
+from . import autoc
 
 # Optional typing
@@ -32 +35 @@
     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
@@ -789 +794 @@
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
+    info.c_code = getattr(kernel_module, 'c_code', None)
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
@@ -804 +810 @@
     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)
 
@@ -814 +818 @@
     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