Changeset d4c33d6 in sasmodels

Timestamp:

Apr 12, 2017 10:50:29 AM (7 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:

597878a

Parents:

535fee6

Message:

send cos-weighted theta values to kernel rather than computing them inside

Files:

• explore/jitter.py (modified) (2 diffs)
• sasmodels/details.py (modified) (7 diffs)
• sasmodels/kernel_iq.c (modified) (5 diffs)
• sasmodels/kernel_iq.cl (modified) (6 diffs)
• sasmodels/kernelpy.py (modified) (3 diffs)
• sasmodels/weights.py (modified) (1 diff)

explore/jitter.py

@@ -191 +191 @@
 def apply_jitter(jitter, points):
     dtheta, dphi, dpsi = jitter
-    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
     return points
 
@@ -257 +257 @@
         ax.cla()
         draw_beam(ax, (0, 0))
-        if 0:
-            draw_jitter(ax, view, jitter)
-        else:
-            draw_jitter(ax, view, (0,0,0))
-            draw_mesh(ax, view, jitter)
+        draw_jitter(ax, view, jitter)
+        #draw_jitter(ax, view, (0,0,0))
+        draw_mesh(ax, view, jitter)
         plt.gcf().canvas.draw()
 

sasmodels/details.py

@@ -15 +15 @@
 
 import numpy as np  # type: ignore
-from numpy import pi, cos, sin
+from numpy import pi, cos, sin, radians
 
 try:
@@ -29 +29 @@
 
 try:
-    from typing import List
+    from typing import List, Tuple, Sequence
 except ImportError:
     pass
 else:
     from .modelinfo import ModelInfo
+    from .modelinfo import ParameterTable
 
 
@@ -53 +54 @@
     coordinates, the total circumference decreases as latitude varies from
     pi r^2 at the equator to 0 at the pole, and the weight associated
-    with a range of phi values needs to be scaled by this circumference.
+    with a range of latitude values needs to be scaled by this circumference.
     This scale factor needs to be updated each time the theta value
     changes.  *theta_par* indicates which of the values in the parameter
@@ -231 +232 @@
     nvalues = kernel.info.parameters.nvalues
     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs]
-    values, weights = zip(*pairs[2:npars+2]) if npars else ((),())
+    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ())
+    weights = correct_theta_weights(kernel.info.parameters, values, weights)
     length = np.array([len(w) for w in weights])
     offset = np.cumsum(np.hstack((0, length)))
@@ -244 +246 @@
     return call_details, data, is_magnetic
 
+def correct_theta_weights(parameters, values, weights):
+    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray]
+    """
+    If there is a theta parameter, update the weights of that parameter so that
+    the cosine weighting required for polar integration is preserved.  Avoid
+    evaluation strictly at the pole, which would otherwise send the weight to
+    zero.
+    """
+    if parameters.theta_offset >= 0:
+        index = parameters.theta_offset+len(parameters.COMMON)
+        theta = values[index]
+        theta_weight = np.minimum(cos(radians(theta)), 1e-6)
+        # copy the weights list so we can update it
+        weights = list(weights)
+        weights[index] = theta_weight*weights[index]
+        weights = tuple(weights)
+    return weights
+
 
 def convert_magnetism(parameters, values):
+    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
     """
     Convert magnetism values from polar to rectangular coordinates.
@@ -255 +276 @@
     scale = mag[:,0]
     if np.any(scale):
-        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
+        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
         cos_theta = cos(theta)
         mag[:, 0] = scale*cos_theta*cos(phi)  # mx
@@ -269 +290 @@
     """
     Create a mesh grid of dispersion parameters and weights.
+
+    pars is a list of pairs (values, weights), where the values are the
+    individual parameter values at which to evaluate the polydispersity
+    distribution and weights are the weights associated with each value.
+
+    Only the volume parameters should be included in this list.  Orientation
+    parameters do not affect the calculation of effective radius or volume
+    ratio.
 
     Returns [p1,p2,...],w where pj is a vector of values for parameter j

sasmodels/kernel_iq.c

@@ -25 +25 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -173 +173 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const int fast_theta = (theta_par == p0);
-  const int slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -220 +209 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -244 +227 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
         #ifdef USE_OPENMP
@@ -304 +284 @@
 #endif // !MAGNETIC
 //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
-          result[q_index] += weight * scattering;
+          result[q_index] += weight0 * scattering;
         }
       }

sasmodels/kernel_iq.cl

@@ -25 +25 @@
     int32_t num_weights;        // total length of the weights vector
     int32_t num_active;         // number of non-trivial pd loops
-    int32_t theta_par;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable (not used)
 } ProblemDetails;
 
@@ -169 +169 @@
 
 
-#if MAX_PD>0
-  const int theta_par = details->theta_par;
-  const bool fast_theta = (theta_par == p0);
-  const bool slow_theta = (theta_par >= 0 && !fast_theta);
-  double spherical_correction = 1.0;
-#else
-  // Note: if not polydisperse the weights cancel and we don't need the
-  // spherical correction.
-  const double spherical_correction = 1.0;
-#endif
-
   int step = pd_start;
 
@@ -217 +206 @@
 #endif
 #if MAX_PD>0
-  if (slow_theta) { // Theta is not in inner loop
-    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
-  }
   while(i0 < n0) {
     local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
-    if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
-    }
 #else
     const double weight0 = 1.0;
@@ -232 +215 @@
 
 //if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); }
-//if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
 
     #ifdef INVALID
@@ -241 +223 @@
       // Note: weight==0 must always be excluded
       if (weight0 > cutoff) {
-        // spherical correction is set at a minimum of 1e-6, otherwise there
-        // would be problems looking at models with theta=90.
-        const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values.table);
+        pd_norm += weight0 * CALL_VOLUME(local_values.table);
 
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
@@ -296 +275 @@
         const double scattering = CALL_IQ(q, q_index, local_values.table);
 #endif // !MAGNETIC
-        this_result += weight * scattering;
+        this_result += weight0 * scattering;
       }
     }

sasmodels/kernelpy.py

@@ -197 +197 @@
 
     pd_norm = 0.0
-    spherical_correction = 1.0
     partial_weight = np.NaN
     weight = np.NaN
 
     p0_par = call_details.pd_par[0]
-    p0_is_theta = (p0_par == call_details.theta_par)
     p0_length = call_details.pd_length[0]
     p0_index = p0_length
@@ -219 +217 @@
             parameters[pd_par] = pd_value[pd_offset+pd_index]
             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
-            if call_details.theta_par >= 0:
-                cor = sin(pi / 180 * parameters[call_details.theta_par])
-                spherical_correction = max(abs(cor), 1e-6)
             p0_index = loop_index%p0_length
 
         weight = partial_weight * pd_weight[p0_offset + p0_index]
         parameters[p0_par] = pd_value[p0_offset + p0_index]
-        if p0_is_theta:
-            cor = cos(pi/180 * parameters[p0_par])
-            spherical_correction = max(abs(cor), 1e-6)
         p0_index += 1
         if weight > cutoff:
@@ -239 +231 @@
 
             # update value and norm
-            weight *= spherical_correction
             total += weight * Iq
             pd_norm += weight * form_volume()

sasmodels/weights.py

@@ -55 +55 @@
         """
         sigma = self.width * center if relative else self.width
+        if not relative:
+            # For orientation, the jitter is relative to 0 not the angle
+            #center = 0
+            pass
         if sigma == 0 or self.npts < 2:
             if lb <= center <= ub: