Changeset a6f9577 in sasmodels

Timestamp:

Mar 21, 2016 7:40:25 AM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

9f4409a

Parents:

39cc3be

Message:

restore support for simple polydispersity (incomplete)

Location:

sasmodels

Files:

• generate.py (modified) (5 diffs)
• kernel_iq.c (modified) (5 diffs)
• kernelcl.py (modified) (1 diff)

sasmodels/generate.py

@@ -593 +593 @@
       width (e.g., radius +/- 10%) rather than absolute distribution
       width (e.g., theta +/- 6 degrees).
+    * *theta_par* is the index of the polar angle polydispersion parameter
+      or -1 if no such parameter exists
     """
     par_set = {}
@@ -600 +602 @@
     par_set['pd'] = [p.name for p in pars if p.type in ('volume', 'orientation')]
     par_set['pd_relative'] = [p.name for p in pars if p.type == 'volume']
+    if 'theta' in par_set['2d']:
+        # find id of theta in parameter set (or whatever variable is
+        # used for spherical normalization during polydispersity...
+        par_set['theta_par'] = [k for k,p in enumerate(pars) if p.name=='theta'][0]
+    else:
+        par_set['theta_par'] = -1
     return par_set
 
@@ -649 +657 @@
         model_info['demo'] = dict((p.name, p.default) for p in pars)
     model_info['has_2d'] = par_type['orientation'] or par_type['magnetic']
+    # Don't use more polydisperse parameters than are available in the model
+    # Note: we can do polydispersity on arbitrary parameters, so it is not
+    # clear that this is a good idea; it does however make the poly_details
+    # code easier to write, so we will leave it in for now.
+    model_info['max_pd'] = min(len(par_type['pd']), MAX_PD)
 
 def mono_details(model_info):
@@ -659 +672 @@
     details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first
     details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
-    details[2*max_pd:3*max_pd] = range(2, max_pd+2)  # pd_offset: skip scale and background
+    details[2*max_pd:3*max_pd] = range(max_pd)       # pd_offset: consecutive 1.0 weights
     details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
     details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
-    details[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset
+    details[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset: skip scale and background
     details[p+1*npars:p+2*npars] = [0]*npars         # no coordination
     #details[p+npars] = 1 # par_coord[0] is coordinated with the first par?
     details[p+2*npars:p+3*npars] = 0 # fast coord with 0
     details[c]   = 1     # fast_coord_count: one fast index
-    details[c+1] = -1  # theta_var: None
-    details[c+2] = 0   # fast_theta: False
+    details[c+1] = -1    # theta_par: None
     return details
 
-def poly_details(model_info, weights, constraints=None):
-    if constraints is not None:
-        # Need to find the independently varying pars and sort them
-        # Need to build a coordination list for the dependent variables
-        # Need to generate a constraints function which takes values
-        # and weights, returning par blocks
-        raise ValueError("Can't handle constraints yet")
-
+def poly_details(model_info, weights):
+    pars = model_info['parameters'][2:]  # skip scale and background
     max_pd = model_info['max_pd']
-    npars = len(model_info['parameters']) - 2
+    npars = len(pars) # scale and background already removed
     p = 5*max_pd
-    c = p + 3*npars
+    constants_offset = p + 3*npars
 
     # Decreasing list of polydispersity lengths
     # Note: the reversing view, x[::-1], does not require a copy
-    idx = np.argsort([len(w) for w in weights])[::-1]
-    details = np.zeros(c + 3, 'int32')
-    details[0*max_pd:1*max_pd] = idx[:max_pd]        # pd_par: arbitrary order; use first
-    details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
-    details[2*max_pd:3*max_pd] = range(2, max_pd+2)  # pd_offset: skip scale and background
-    details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
-    details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
-    details[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset
-    details[p+1*npars:p+2*npars] = [0]*npars         # no coordination
-    #details[p+npars] = 1 # par_coord[0] is coordinated with the first par?
-    details[p+2*npars:p+3*npars] = 0 # fast coord with 0
-    details[c]   = 1     # fast_coord_count: one fast index
-    details[c+1] = -1  # theta_var: None
-    details[c+2] = 0   # fast_theta: False
+    pd_length = np.array([len(w) for w in weights])
+    pd_offset = np.cumsum(np.hstack((0, pd_length)))
+    pd_isvol = np.array([p.type=='volume' for p in pars])
+    idx = np.argsort(pd_length)[::-1][:max_pd]
+    pd_stride = np.cumprod(np.hstack((1, np.maximum(pd_length[idx][:-1],1))))
+    par_offsets = np.cumsum(np.hstack((2, np.maximum(pd_length, 1))))[:-1]
+    theta_par = model_info['theta_par']
+    if theta_par >= 0 and pd_length[theta_par] <= 1:
+        theta_par = -1
+
+    details = np.empty(constants_offset + 3, 'int32')
+    details[0*max_pd:1*max_pd] = idx             # pd_par
+    details[1*max_pd:2*max_pd] = pd_length[idx]
+    details[2*max_pd:3*max_pd] = pd_offset[idx]
+    details[3*max_pd:4*max_pd] = pd_stride
+    details[4*max_pd:5*max_pd] = pd_isvol[idx]
+    details[p+0*npars:p+1*npars] = par_offsets
+    details[p+1*npars:p+2*npars] = 0  # no coordination for most
+    details[p+2*npars:p+3*npars] = 0  # no fast coord with 0
+    coord_offset = p+1*pars
+    for k,parameter_num in enumerate(idx):
+        details[coord_offset+parameter_num] = 2**k
+    details[constants_offset]   = 1   # fast_coord_count: one fast index
+    details[constants_offset+1] = theta_par
     return details
 
-
-
-    raise ValueError("polydispersity not supported")
+def constrained_poly_details(model_info, weights, constraints):
+    # Need to find the independently varying pars and sort them
+    # Need to build a coordination list for the dependent variables
+    # Need to generate a constraints function which takes values
+    # and weights, returning par blocks
+    raise NotImplementedError("Can't handle constraints yet")
 
 
@@ -777 +796 @@
         name = " ".join(w.capitalize() for w in kernel_id.split('_'))
     model_info = dict(
-        max_pd=MAX_PD,
         id=kernel_id,  # string used to load the kernel
         filename=abspath(kernel_module.__file__),

sasmodels/kernel_iq.c

@@ -18 +18 @@
 
 typedef struct {
-    int32_t pd_par[MAX_PD];     // index of the nth polydispersity variable
+    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
     int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
     int32_t pd_offset[MAX_PD];  // offset of pd weights in the par & weight vector
@@ -25 +25 @@
     int32_t par_offset[NPARS];  // offset of par values in the par & weight vector
     int32_t par_coord[NPARS];   // polydispersity coordination bitvector
-    int32_t fast_coord_index[NPARS]; // index of the fast coordination parameters
+    int32_t fast_coord_pars[NPARS]; // ids of the fast coordination parameters
     int32_t fast_coord_count;   // number of parameters coordinated with pd 1
     int32_t theta_var;          // id of spherical correction variable
-    int32_t fast_theta;         // true if spherical correction depends on pd 1
 } ProblemDetails;
 
@@ -55 +54 @@
   double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
 
-  local int offset[NPARS-2];
+  local int offset[NPARS];  // NPARS excludes scale/background
 
 #if 1 // defined(USE_SHORTCUT_OPTIMIZATION)
@@ -155 +154 @@
         spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var]));
       }
-      if (!problem->fast_theta) {
+      if (problem->theta_var == problem->pd_par[0]) {
         weight *= spherical_correction;
       }
@@ -167 +166 @@
       if (problem->pd_isvol[0]) vol_weight *= wi;
       for (int k=0; k < problem->fast_coord_count; k++) {
-        pvec[problem->fast_coord_index[k]]
-            = pars[offset[problem->fast_coord_index[k]]++];
-      }
-      if (problem->fast_theta) {
+        pvec[problem->fast_coord_pars[k]]
+            = pars[offset[problem->fast_coord_pars[k]]++];
+      }
+      if (problem->theta_var ==problem->pd_par[0]) {
         weight *= fabs(cos(M_PI_180*pvec[problem->theta_var]));
       }

sasmodels/kernelcl.py

@@ -448 +448 @@
         self._need_release = [self.loops_b, self.res_b, self.q_input]
 
-    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
+    def __call__(self, details, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64