Changeset 6aee3ab in sasmodels

Timestamp:

Nov 2, 2017 5:24:41 PM (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:

ce99754

Parents:

9d9fcbd

Message:

make sure that jitter angle defaults to zero if no polydispersity

Location:

sasmodels

Files:

• details.py (modified) (1 diff)
• kernel_iq.c (modified) (12 diffs)
• modelinfo.py (modified) (1 diff)

sasmodels/details.py

@@ -193 +193 @@
     #print("off:",offset)
 
-    # Check that we arn't using too many polydispersity loops
+    # Check that we aren't using too many polydispersity loops
     num_active = np.sum(length > 1)
     max_pd = model_info.parameters.max_pd

sasmodels/kernel_iq.c

@@ -13 +13 @@
 // NOTE: the following macros are defined in generate.py:
 //
-//  MAX_PD : the maximum number of dispersity loops allowed
+//  MAX_PD : the maximum number of dispersity loops allowed for this model,
+//      which will be at most modelinfo.MAX_PD.
 //  NUM_PARS : the number of parameters in the parameter table
 //  NUM_VALUES : the number of values to skip at the start of the
@@ -282 +283 @@
 #endif
 
-  // Storage for the current parameter values.  These will be updated as we
-  // walk the dispersity mesh.
+  // ** Fill in the local values table **
+  // Storage for the current parameter values.
+  // These will be updated as we walk the dispersity mesh.
   ParameterBlock local_values;
-
+  //   values[0] is scale
+  //   values[1] is background
+  #ifdef USE_OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int i=0; i < NUM_PARS; i++) {
+    local_values.vector[i] = values[2+i];
+    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
+  }
+  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
+
+  // ** Precompute magnatism values **
 #if defined(MAGNETIC) && NUM_MAGNETIC>0
   // Location of the sld parameters in the parameter vector.
@@ -302 +316 @@
 #endif // MAGNETIC
 
-  // Fill in the initial variables
-  //   values[0] is scale
-  //   values[1] is background
-  #ifdef USE_OPENMP
-  #pragma omp parallel for
-  #endif
-  for (int i=0; i < NUM_PARS; i++) {
-    local_values.vector[i] = values[2+i];
-//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
-  }
-//if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
-//if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
-
+  // ** Fill in the initial results **
   // If pd_start is zero that means that we are starting a new calculation,
   // and must initialize the result to zero.  Otherwise, we are restarting
@@ -336 +338 @@
       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
     }
-//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 #endif // !USE_OPENCL
 
@@ -390 +392 @@
       PD_OPEN(0,1)
 
-      ... main loop body ...
+      // ... main loop body ...
+      APPLY_PROJECTION    // convert jitter values to spherical coords
+      BUILD_ROTATION      // construct the rotation matrix qxy => qabc
+      for each q
+          FETCH_Q         // set qx,qy from the q input vector
+          APPLY_ROTATION  // convert qx,qy to qa,qb,qc
+          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
 
       ++step;  // increment counter representing position in dispersity mesh
@@ -412 +420 @@
 */
 
-// Define looping variables
-#define PD_INIT(_LOOP) \
-  const int n##_LOOP = details->pd_length[_LOOP]; \
-  const int p##_LOOP = details->pd_par[_LOOP]; \
-  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
-  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
-  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
-
-// Jump into the middle of the dispersity loop
-#define PD_OPEN(_LOOP,_OUTER) \
-  while (i##_LOOP < n##_LOOP) { \
-    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
-    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
-
-// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
-#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
-#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
-
-// Close out the loop
-#define PD_CLOSE(_LOOP) \
-    if (step >= pd_stop) break; \
-    ++i##_LOOP; \
-  } \
-  i##_LOOP = 0;
-
-// The main loop body is essentially:
-//
-//    BUILD_ROTATION      // construct the rotation matrix qxy => qabc
-//    for each q
-//        FETCH_Q         // set qx,qy from the q input vector
-//        APPLY_ROTATION  // convert qx,qy to qa,qb,qc
-//        CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
-//
+
+// ** prepare inner loops **
+
 // Depending on the shape type (radial, axial, triaxial), the variables
-// and calling parameters will be slightly different.  These macros
-// capture the differences in one spot so the rest of the code is easier
-// to read.
+// and calling parameters in the loop body will be slightly different.
+// Macros capture the differences in one spot so the rest of the code
+// is easier to read. The code below both declares variables for the
+// inner loop and defines the macros that use them.
 
 #if defined(CALL_IQ)
@@ -486 +465 @@
   // psi and dpsi are only for IQ_ABC, so they are processed here.
   const double psi = values[details->theta_par+4];
+  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)
@@ -501 +481 @@
   const double theta = values[details->theta_par+2];
   const double phi = values[details->theta_par+3];
+  // Make sure jitter angle defaults to zero if there is no jitter distribution
+  local_values.table.theta = 0.;
+  local_values.table.phi = 0.;
   // The "jitter" angles (dtheta, dphi, dpsi) are stored with the
   // dispersity values and copied to the local parameter table as
@@ -523 +506 @@
 #endif // !spherosymmetric projection
 
+// ** define looping macros **
+
+// Define looping variables
+#define PD_INIT(_LOOP) \
+  const int n##_LOOP = details->pd_length[_LOOP]; \
+  const int p##_LOOP = details->pd_par[_LOOP]; \
+  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
+  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
+
+// Jump into the middle of the dispersity loop
+#define PD_OPEN(_LOOP,_OUTER) \
+  while (i##_LOOP < n##_LOOP) { \
+    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
+    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
+
+// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
+#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
+#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
+
+// Close out the loop
+#define PD_CLOSE(_LOOP) \
+    if (step >= pd_stop) break; \
+    ++i##_LOOP; \
+  } \
+  i##_LOOP = 0;
 
 // ====== construct the loops =======
@@ -537 +546 @@
 int step = pd_start;
 
+// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
+
 // define looping variables
 #if MAX_PD>4
@@ -571 +582 @@
   PD_OPEN(0,1)
 #endif
-
 
 //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");}
@@ -672 +682 @@
 #endif // !USE_OPENCL
 
-// clear the macros in preparation for the next kernel
+// ** clear the macros in preparation for the next kernel **
 #undef PD_INIT
 #undef PD_OPEN

sasmodels/modelinfo.py

@@ -30 +30 @@
     TestCondition = Tuple[ParameterSetUser, TestInput, TestValue]
 
-MAX_PD = 4 #: Maximum number of simultaneously polydisperse parameters
+# If MAX_PD changes, need to change the loop macros in kernel_iq.c
+MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters
 
 # assumptions about common parameters exist throughout the code, such as: