Changes in / [510277d:376b0ee] 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 (deleted)
• models/core_shell_bicelle_elliptical_belt_rough.py (deleted)
• models/img/core_shell_bicelle_belt_geometry.png (deleted)
• models/img/core_shell_bicelle_belt_parameters.png (deleted)
• product.py (modified) (2 diffs)
• sasview_model.py (modified) (7 diffs)

sasmodels/details.py

@@ -193 +193 @@
     #print("off:",offset)
 
-    # Check that we aren't using too many polydispersity loops
+    # Check that we arn'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.
     """
-    _, dispersity, weight = zip(*mesh)
+    value, 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 for this model,
-//      which will be at most modelinfo.MAX_PD.
+//  MAX_PD : the maximum number of dispersity loops allowed
 //  NUM_PARS : the number of parameters in the parameter table
 //  NUM_VALUES : the number of values to skip at the start of the
@@ -283 +282 @@
 #endif
 
-  // ** Fill in the local values table **
-  // Storage for the current parameter values.
-  // These will be updated as we walk the dispersity mesh.
+  // 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.
@@ -316 +302 @@
 #endif // MAGNETIC
 
-  // ** Fill in the initial results **
+  // 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);
+
   // 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
@@ -338 +336 @@
       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
 
@@ -392 +390 @@
       PD_OPEN(0,1)
 
-      // ... 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, ...)
+      ... main loop body ...
 
       ++step;  // increment counter representing position in dispersity mesh
@@ -420 +412 @@
 */
 
-
-// ** prepare inner loops **
-
+// 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, ...)
+//
 // Depending on the shape type (radial, axial, triaxial), the variables
-// 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.
+// 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.
 
 #if defined(CALL_IQ)
@@ -465 +486 @@
   // 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)
@@ -481 +501 @@
   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
@@ -506 +523 @@
 #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 =======
@@ -546 +537 @@
 int step = pd_start;
 
-// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
-
 // define looping variables
 #if MAX_PD>4
@@ -582 +571 @@
   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");}
@@ -682 +672 @@
 #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]
 
-# If MAX_PD changes, need to change the loop macros in kernel_iq.c
-MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters
+MAX_PD = 4 #: 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: {}
@@ -279 +267 @@
     weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
     npars = model_info.parameters.npars
-    # 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])
+    pairs = [(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)
@@ -680 +678 @@
         result = calculator(call_details, values, cutoff=self.cutoff,
                             magnetic=is_magnetic)
-        #print("result", result)
         self._intermediate_results = getattr(calculator, 'results', None)
         calculator.release()
@@ -764 +761 @@
                 return self.multiplicity, [self.multiplicity], [1.0]
             else:
-                # For hidden parameters use default values.  This sets
-                # scale=1 and background=0 for structure factors
+                # For hidden parameters use the default value.
                 default = self._model_info.parameters.defaults.get(par.name, np.NaN)
-                return default, [default], [1.0]
+                return [default], [1.0]
         elif par.polydisperse:
             value = self.params[par.name]
@@ -780 +776 @@
         else:
             value = self.params[par.name]
-            return value, [value], [1.0]
+            return value, [value if par.relative_pd else 0.0], [1.0]
 
 def test_cylinder():
@@ -798 +794 @@
     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("cylinder*hatyer_msa returns null")
+        raise ValueError("hardsphere returns null")
 
 def test_rpa():
@@ -868 +850 @@
 if __name__ == "__main__":
     print("cylinder(0.1,0.1)=%g"%test_cylinder())
-    #test_product()
-    #test_structure_factor()
-    #print("rpa:", test_rpa())
     #test_empty_distribution()