source: sasmodels/sasmodels/mixture.py @ fb9a3b6

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since fb9a3b6 was fb9a3b6, checked in by lewis, 7 years ago

Allow mixture models to do multiplication as well as addition

  • Property mode set to 100644
File size: 9.0 KB
Line 
1"""
2Mixture model
3-------------
4
5The product model multiplies the structure factor by the form factor,
6modulated by the effective radius of the form.  The resulting model
7has a attributes of both the model description (with parameters, etc.)
8and the module evaluator (with call, release, etc.).
9
10To use it, first load form factor P and structure factor S, then create
11*ProductModel(P, S)*.
12"""
13from __future__ import print_function
14
15from copy import copy
16import numpy as np  # type: ignore
17
18from .modelinfo import Parameter, ParameterTable, ModelInfo
19from .kernel import KernelModel, Kernel
20from .details import make_details
21
22try:
23    from typing import List
24except ImportError:
25    pass
26
27def make_mixture_info(parts, operation='+'):
28    # type: (List[ModelInfo]) -> ModelInfo
29    """
30    Create info block for mixture model.
31    """
32    flatten = []
33    for part in parts:
34        if part.composition and part.composition[0] == 'mixture':
35            flatten.extend(part.composition[1])
36        else:
37            flatten.append(part)
38    parts = flatten
39
40    # Build new parameter list
41    combined_pars = []
42    demo = {}
43    for k, part in enumerate(parts):
44        # Parameter prefix per model, A_, B_, ...
45        # Note that prefix must also be applied to id and length_control
46        # to support vector parameters
47        prefix = chr(ord('A')+k) + '_'
48        if operation == '+':
49            # If model is a sum model, each constituent model gets its own scale parameter
50            scale =  Parameter(prefix+'scale', default=1.0,
51                            description="model intensity for " + part.name)
52            combined_pars.append(scale)
53        for p in part.parameters.kernel_parameters:
54            p = copy(p)
55            p.name = prefix + p.name
56            p.id = prefix + p.id
57            if p.length_control is not None:
58                p.length_control = prefix + p.length_control
59            combined_pars.append(p)
60        demo.update((prefix+k, v) for k, v in part.demo.items()
61                    if k != "background")
62    #print("pars",combined_pars)
63    parameters = ParameterTable(combined_pars)
64    parameters.max_pd = sum(part.parameters.max_pd for part in parts)
65
66    model_info = ModelInfo()
67    model_info.id = operation.join(part.id for part in parts)
68    model_info.operation = operation
69    model_info.name = operation.join(part.name for part in parts)
70    model_info.filename = None
71    model_info.title = 'Mixture model with ' + model_info.name
72    model_info.description = model_info.title
73    model_info.docs = model_info.title
74    model_info.category = "custom"
75    model_info.parameters = parameters
76    #model_info.single = any(part['single'] for part in parts)
77    model_info.structure_factor = False
78    model_info.variant_info = None
79    #model_info.tests = []
80    #model_info.source = []
81    # Iq, Iqxy, form_volume, ER, VR and sesans
82    # Remember the component info blocks so we can build the model
83    model_info.composition = ('mixture', parts)
84    model_info.demo = demo
85    return model_info
86
87
88class MixtureModel(KernelModel):
89    def __init__(self, model_info, parts):
90        # type: (ModelInfo, List[KernelModel]) -> None
91        self.info = model_info
92        self.parts = parts
93
94    def make_kernel(self, q_vectors):
95        # type: (List[np.ndarray]) -> MixtureKernel
96        # Note: may be sending the q_vectors to the n times even though they
97        # are only needed once.  It would mess up modularity quite a bit to
98        # handle this optimally, especially since there are many cases where
99        # separate q vectors are needed (e.g., form in python and structure
100        # in opencl; or both in opencl, but one in single precision and the
101        # other in double precision).
102        kernels = [part.make_kernel(q_vectors) for part in self.parts]
103        return MixtureKernel(self.info, kernels)
104
105    def release(self):
106        # type: () -> None
107        """
108        Free resources associated with the model.
109        """
110        for part in self.parts:
111            part.release()
112
113
114class MixtureKernel(Kernel):
115    def __init__(self, model_info, kernels):
116        # type: (ModelInfo, List[Kernel]) -> None
117        self.dim = kernels[0].dim
118        self.info =  model_info
119        self.kernels = kernels
120        self.dtype = self.kernels[0].dtype
121        self.operation = model_info.operation
122        self.results = []  # type: List[np.ndarray]
123
124    def __call__(self, call_details, values, cutoff, magnetic):
125        # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray
126        scale, background = values[0:2]
127        total = 0.0
128        # remember the parts for plotting later
129        self.results = []  # type: List[np.ndarray]
130        parts = MixtureParts(self.info, self.kernels, call_details, values)
131        for kernel, kernel_details, kernel_values in parts:
132            #print("calling kernel", kernel.info.name)
133            result = kernel(kernel_details, kernel_values, cutoff, magnetic)
134            result = np.array(result).astype(kernel.dtype)
135            # print(kernel.info.name, result)
136            if self.operation == '+':
137                total += result
138            elif self.operation == '*':
139                if np.all(total) == 0.0:
140                    total = result
141                else:
142                    total *= result
143            self.results.append(result)
144
145        return scale*total + background
146
147    def release(self):
148        # type: () -> None
149        for k in self.kernels:
150            k.release()
151
152
153class MixtureParts(object):
154    def __init__(self, model_info, kernels, call_details, values):
155        # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None
156        self.model_info = model_info
157        self.parts = model_info.composition[1]
158        self.kernels = kernels
159        self.call_details = call_details
160        self.values = values
161        self.spin_index = model_info.parameters.npars + 2
162        #call_details.show(values)
163
164    def __iter__(self):
165        # type: () -> PartIterable
166        self.part_num = 0
167        self.par_index = 2
168        self.mag_index = self.spin_index + 3
169        return self
170
171    def next(self):
172        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray]
173        if self.part_num >= len(self.parts):
174            raise StopIteration()
175        info = self.parts[self.part_num]
176        kernel = self.kernels[self.part_num]
177        call_details = self._part_details(info, self.par_index)
178        values = self._part_values(info, self.par_index, self.mag_index)
179        values = values.astype(kernel.dtype)
180        #call_details.show(values)
181
182        self.part_num += 1
183        self.par_index += info.parameters.npars
184        if self.model_info.operation == '+':
185            self.par_index += 1 # Account for each constituent model's scale param
186        self.mag_index += 3 * len(info.parameters.magnetism_index)
187
188        return kernel, call_details, values
189
190    def _part_details(self, info, par_index):
191        # type: (ModelInfo, int) -> CallDetails
192        full = self.call_details
193        # par_index is index into values array of the current parameter,
194        # which includes the initial scale and background parameters.
195        # We want the index into the weight length/offset for each parameter.
196        # Exclude the initial scale and background, so subtract two. If we're
197        # building an addition model, each component has its own scale factor
198        # which we need to skip when constructing the details for the kernel, so
199        # add one, giving a net subtract one.
200        diff = 1 if self.model_info.operation == '+' else 2
201        index = slice(par_index - diff, par_index - diff + info.parameters.npars)
202        length = full.length[index]
203        offset = full.offset[index]
204        # The complete weight vector is being sent to each part so that
205        # offsets don't need to be adjusted.
206        part = make_details(info, length, offset, full.num_weights)
207        return part
208
209    def _part_values(self, info, par_index, mag_index):
210        # type: (ModelInfo, int, int) -> np.ndarray
211        # Set each constituent model's scale to 1 if this is a multiplication model
212        scale = self.values[par_index] if self.model_info.operation == '+' else 1.0
213        diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model
214        pars = self.values[par_index + diff:par_index + info.parameters.npars + diff]
215        nmagnetic = len(info.parameters.magnetism_index)
216        if nmagnetic:
217            spin_state = self.values[self.spin_index:self.spin_index + 3]
218            mag_index = self.values[mag_index:mag_index + 3 * nmagnetic]
219        else:
220            spin_state = []
221            mag_index = []
222        nvalues = self.model_info.parameters.nvalues
223        nweights = self.call_details.num_weights
224        weights = self.values[nvalues:nvalues+2*nweights]
225        zero = self.values.dtype.type(0.)
226        values = [[scale, zero], pars, spin_state, mag_index, weights]
227        # Pad value array to a 32 value boundary
228        spacer = (32 - sum(len(v) for v in values)%32)%32
229        values.append([zero]*spacer)
230        values = np.hstack(values).astype(self.kernels[0].dtype)
231        return values
Note: See TracBrowser for help on using the repository browser.