Changes in / [f475f90:d0462cf] in sasmodels

Location:

sasmodels

Files:

• details.py (modified) (2 diffs)
• kernel_iq.c (modified) (12 diffs)
• modelinfo.py (modified) (1 diff)
• models/core_shell_bicelle_elliptical_belt_rough.c (added)
• models/core_shell_bicelle_elliptical_belt_rough.py (added)
• models/img/core_shell_bicelle_belt_geometry.png (added)
• models/img/core_shell_bicelle_belt_parameters.png (added)
• product.py (modified) (2 diffs)
• sasview_model.py (modified) (7 diffs)

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
@@ -318 +318 @@
     parameter set in the vector.
     """
-    value, dispersity, weight = zip(*mesh)
+    _, dispersity, weight = zip(*mesh)
     #weight = [w if len(w)>0 else [1.] for w in weight]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])

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:

sasmodels/product.py

@@ -101 +101 @@
     model_info.composition = ('product', [p_info, s_info])
     model_info.control = p_info.control
+    model_info.hidden = p_info.hidden
+    if getattr(p_info, 'profile', None) is not None:
+        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
+        def profile(**kwargs):
+            # extract the profile args
+            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
+            return p_info.profile(**kwargs)
+    else:
+        profile = None
+    model_info.profile = profile
+    model_info.profile_axes = p_info.profile_axes
+
     # TODO: delegate random to p_info, s_info
     #model_info.random = lambda: {}
@@ -267 +279 @@
     weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
     npars = model_info.parameters.npars
-    pairs = [(value[offset:offset+length], weight[offset:offset+length])
+    # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the
+    # dispersion mesh code doesn't actually care about the nominal parameter
+    # value, p, so set it to None.
+    pairs = [(None, value[offset:offset+length], weight[offset:offset+length])
              for p, offset, length
              in zip(model_info.parameters.call_parameters[2:2+npars],

sasmodels/sasview_model.py

@@ -205 +205 @@
                                            structure_factor._model_info)
     ConstructedModel = make_model_from_info(model_info)
-    return ConstructedModel(form_factor.multiplicity)    
+    return ConstructedModel(form_factor.multiplicity)
 
 
@@ -673 +673 @@
         #call_details.show()
         #print("pairs", pairs)
+        #for k, p in enumerate(self._model_info.parameters.call_parameters):
+        #    print(k, p.name, *pairs[k])
         #print("params", self.params)
         #print("values", values)
@@ -678 +680 @@
         result = calculator(call_details, values, cutoff=self.cutoff,
                             magnetic=is_magnetic)
+        #print("result", result)
         self._intermediate_results = getattr(calculator, 'results', None)
         calculator.release()
@@ -761 +764 @@
                 return self.multiplicity, [self.multiplicity], [1.0]
             else:
-                # For hidden parameters use the default value.
+                # For hidden parameters use default values.  This sets
+                # scale=1 and background=0 for structure factors
                 default = self._model_info.parameters.defaults.get(par.name, np.NaN)
-                return [default], [1.0]
+                return default, [default], [1.0]
         elif par.polydisperse:
             value = self.params[par.name]
@@ -776 +780 @@
         else:
             value = self.params[par.name]
-            return value, [value if par.relative_pd else 0.0], [1.0]
+            return value, [value], [1.0]
 
 def test_cylinder():
@@ -794 +798 @@
     Model = _make_standard_model('hardsphere')
     model = Model()
+    value2d = model.evalDistribution([0.1, 0.1])
+    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
+    #print("hardsphere", value1d, value2d)
+    if np.isnan(value1d) or np.isnan(value2d):
+        raise ValueError("hardsphere returns nan")
+
+def test_product():
+    # type: () -> float
+    """
+    Test that 2-D hardsphere model runs and doesn't produce NaN.
+    """
+    S = _make_standard_model('hayter_msa')()
+    P = _make_standard_model('cylinder')()
+    model = MultiplicationModel(P, S)
     value = model.evalDistribution([0.1, 0.1])
     if np.isnan(value):
-        raise ValueError("hardsphere returns null")
+        raise ValueError("cylinder*hatyer_msa returns null")
 
 def test_rpa():
@@ -850 +868 @@
 if __name__ == "__main__":
     print("cylinder(0.1,0.1)=%g"%test_cylinder())
+    #test_product()
+    #test_structure_factor()
+    #print("rpa:", test_rpa())
     #test_empty_distribution()