Changeset fb9a3b6 in sasmodels

Timestamp:

Aug 31, 2017 12:07:59 PM (7 years ago)

Author:

lewis

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

bcdd6c9

Parents:

573ffab

Message:

Allow mixture models to do multiplication as well as addition

Location:

sasmodels

Files:

• core.py (modified) (4 diffs)
• mixture.py (modified) (8 diffs)

sasmodels/core.py

@@ -126 +126 @@
 
 
-def load_model_info(model_name):
+def load_model_info(model_name, force_mixture=False):
     # type: (str) -> modelinfo.ModelInfo
     """
@@ -134 +134 @@
     such as sphere*hardsphere or sphere+cylinder.
 
+    *force_mixture* if true, MixtureModel will be used for combining models.
+    Otherwise either MixtureModel will be used for addition and ProductModel
+    will be used for multiplication
+
     This returns a handle to the module defining the model.  This can be
     used with functions in generate to build the docs or extract model info.
@@ -139 +143 @@
     parts = model_name.split('+')
     if len(parts) > 1:
+        # Always use MixtureModel for addition
         model_info_list = [load_model_info(p) for p in parts]
         return mixture.make_mixture_info(model_info_list)
@@ -144 +149 @@
     parts = model_name.split('*')
     if len(parts) > 1:
+        if force_mixture:
+            # Use MixtureModel for multiplication if forced
+            model_info_list = [load_model_info(p) for p in parts]
+            return mixture.make_mixture_info(model_info_list, operation='*')
         if len(parts) > 2:
             raise ValueError("use P*S to apply structure factor S to model P")
+        # Use ProductModel
         P_info, Q_info = [load_model_info(p) for p in parts]
         return product.make_product_info(P_info, Q_info)

sasmodels/mixture.py

@@ -25 +25 @@
     pass
 
-def make_mixture_info(parts):
+def make_mixture_info(parts, operation='+'):
     # type: (List[ModelInfo]) -> ModelInfo
     """
@@ -46 +46 @@
         # to support vector parameters
         prefix = chr(ord('A')+k) + '_'
-        scale =  Parameter(prefix+'scale', default=1.0,
-                           description="model intensity for " + part.name)
-        combined_pars.append(scale)
+        if operation == '+':
+            # If model is a sum model, each constituent model gets its own scale parameter
+            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)
@@ -63 +65 @@
 
     model_info = ModelInfo()
-    model_info.id = '+'.join(part.id for part in parts)
-    model_info.name = ' + '.join(part.name for part in parts)
+    model_info.id = operation.join(part.id for part in parts)
+    model_info.operation = operation
+    model_info.name = operation.join(part.name for part in parts)
     model_info.filename = None
     model_info.title = 'Mixture model with ' + model_info.name
@@ -116 +119 @@
         self.kernels = kernels
         self.dtype = self.kernels[0].dtype
+        self.operation = model_info.operation
         self.results = []  # type: List[np.ndarray]
 
@@ -124 +128 @@
         # remember the parts for plotting later
         self.results = []  # type: List[np.ndarray]
-        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
+            result = np.array(result).astype(kernel.dtype)
+            # print(kernel.info.name, result)
+            if self.operation == '+':
+                total += result
+            elif self.operation == '*':
+                if np.all(total) == 0.0:
+                    total = result
+                else:
+                    total *= result
             self.results.append(result)
 
@@ -171 +181 @@
 
         self.part_num += 1
-        self.par_index += info.parameters.npars + 1
+        self.par_index += info.parameters.npars
+        if self.model_info.operation == '+':
+            self.par_index += 1 # Account for each constituent model's scale param
         self.mag_index += 3 * len(info.parameters.magnetism_index)
 
@@ -182 +194 @@
         # 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)
+        # Exclude the initial scale and background, so subtract two. If we're
+        # building an addition model, 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.
+        diff = 1 if self.model_info.operation == '+' else 2
+        index = slice(par_index - diff, par_index - diff + info.parameters.npars)
         length = full.length[index]
         offset = full.offset[index]
@@ -196 +209 @@
     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]
+        # Set each constituent model's scale to 1 if this is a multiplication model
+        scale = self.values[par_index] if self.model_info.operation == '+' else 1.0
+        diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model
+        pars = self.values[par_index + diff:par_index + info.parameters.npars + diff]
         nmagnetic = len(info.parameters.magnetism_index)
         if nmagnetic: