Changeset ac98886 in sasmodels for sasmodels/mixture.py

Timestamp:

Aug 14, 2016 11:15:15 PM (8 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:

725ee36

Parents:

bde38b5

Message:

restore support for mixture models

File:

• sasmodels/mixture.py (modified) (8 diffs)

sasmodels/mixture.py

@@ -11 +11 @@
 *ProductModel(P, S)*.
 """
+from __future__ import print_function
+
 from copy import copy
 import numpy as np  # type: ignore
@@ -16 +18 @@
 from .modelinfo import Parameter, ParameterTable, ModelInfo
 from .kernel import KernelModel, Kernel
+from .details import make_details
 
 try:
     from typing import List
-    from .details import CallDetails
 except ImportError:
     pass
@@ -43 +45 @@
         # to support vector parameters
         prefix = chr(ord('A')+k) + '_'
-        combined_pars.append(Parameter(prefix+'scale'))
+        scale =  Parameter(prefix+'scale', default=1.0,
+                           description="model intensity for " + part.name)
+        combined_pars.append(scale)
         for p in part.parameters.kernel_parameters:
             p = copy(p)
@@ -51 +55 @@
                 p.length_control = prefix + p.length_control
             combined_pars.append(p)
+    #print("pars",combined_pars)
     parameters = ParameterTable(combined_pars)
+    parameters.max_pd = sum(part.parameters.max_pd for part in parts)
 
     model_info = ModelInfo()
@@ -70 +76 @@
     # Remember the component info blocks so we can build the model
     model_info.composition = ('mixture', parts)
+    model_info.demo = {}
+    return model_info
 
 
@@ -78 +86 @@
         self.parts = parts
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         # type: (List[np.ndarray]) -> MixtureKernel
         # Note: may be sending the q_vectors to the n times even though they
@@ -104 +112 @@
         self.info =  model_info
         self.kernels = kernels
-
-    def __call__(self, call_details, value, weight, cutoff):
+        self.dtype = self.kernels[0].dtype
+
+    def __call__(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray
-        scale, background = value[0:2]
+        scale, background = values[0:2]
         total = 0.0
         # remember the parts for plotting later
         self.results = []
-        for kernel, kernel_details in zip(self.kernels, call_details.parts):
-            part_result = kernel(kernel_details, value, weight, cutoff)
-            total += part_result
-            self.results.append(part_result)
+        offset = 2 # skip scale & background
+        parts = MixtureParts(self.info, self.kernels, call_details, values)
+        for kernel, kernel_details, kernel_values in parts:
+            #print("calling kernel", kernel.info.name)
+            result = kernel(kernel_details, kernel_values, cutoff, magnetic)
+            #print(kernel.info.name, result)
+            total += result
+            self.results.append(result)
 
         return scale*total + background
@@ -123 +136 @@
             k.release()
 
+
+class MixtureParts(object):
+    def __init__(self, model_info, kernels, call_details, values):
+        self.model_info = model_info
+        self.parts = model_info.composition[1]
+        self.kernels = kernels
+        self.call_details = call_details
+        self.values = values
+        self.spin_index = model_info.parameters.npars + 2
+        #call_details.show(values)
+
+    def __iter__(self):
+        # type: () -> PartIterable
+        self.part_num = 0
+        self.par_index = 2
+        self.mag_index = self.spin_index + 3
+        return self
+
+    def next(self):
+        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray]
+        if self.part_num >= len(self.parts):
+            raise StopIteration()
+        info = self.parts[self.part_num]
+        kernel = self.kernels[self.part_num]
+        call_details = self._part_details(info, self.par_index)
+        values = self._part_values(info, self.par_index, self.mag_index)
+        values = values.astype(kernel.dtype)
+        #call_details.show(values)
+
+        self.part_num += 1
+        self.par_index += info.parameters.npars + 1
+        self.mag_index += 3 * len(info.parameters.magnetism_index)
+
+        return kernel, call_details, values
+
+    def _part_details(self, info, par_index):
+        # type: (ModelInfo, int) -> CallDetails
+        full = self.call_details
+        # par_index is index into values array of the current parameter,
+        # which includes the initial scale and background parameters.
+        # We want the index into the weight length/offset for each parameter.
+        # Exclude the initial scale and background, so subtract two, but each
+        # component has its own scale factor which we need to skip when
+        # constructing the details for the kernel, so add one, giving a
+        # net subtract one.
+        index = slice(par_index - 1, par_index - 1 + info.parameters.npars)
+        length = full.length[index]
+        offset = full.offset[index]
+        # The complete weight vector is being sent to each part so that
+        # offsets don't need to be adjusted.
+        part = make_details(info, length, offset, full.num_weights)
+        return part
+
+    def _part_values(self, info, par_index, mag_index):
+        # type: (ModelInfo, int, int) -> np.ndarray
+        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1])
+        scale = self.values[par_index]
+        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1]
+        nmagnetic = len(info.parameters.magnetism_index)
+        if nmagnetic:
+            spin_state = self.values[self.spin_index:self.spin_index + 3]
+            mag_index = self.values[mag_index:mag_index + 3 * nmagnetic]
+        else:
+            spin_state = []
+            mag_index = []
+        nvalues = self.model_info.parameters.nvalues
+        nweights = self.call_details.num_weights
+        weights = self.values[nvalues:nvalues+2*nweights]
+        zero = self.values.dtype.type(0.)
+        values = [[scale, zero], pars, spin_state, mag_index, weights]
+        # Pad value array to a 32 value boundary
+        spacer = (32 - sum(len(v) for v in values)%32)%32
+        values.append([zero]*spacer)
+        values = np.hstack(values).astype(self.kernels[0].dtype)
+        return values