1 | """ |
---|
2 | Mixture model |
---|
3 | ------------- |
---|
4 | |
---|
5 | The product model multiplies the structure factor by the form factor, |
---|
6 | modulated by the effective radius of the form. The resulting model |
---|
7 | has a attributes of both the model description (with parameters, etc.) |
---|
8 | and the module evaluator (with call, release, etc.). |
---|
9 | |
---|
10 | To use it, first load form factor P and structure factor S, then create |
---|
11 | *ProductModel(P, S)*. |
---|
12 | """ |
---|
13 | from __future__ import print_function |
---|
14 | |
---|
15 | from copy import copy |
---|
16 | import numpy as np # type: ignore |
---|
17 | |
---|
18 | from .modelinfo import Parameter, ParameterTable, ModelInfo |
---|
19 | from .kernel import KernelModel, Kernel |
---|
20 | from .details import make_details |
---|
21 | |
---|
22 | # pylint: disable=unused-import |
---|
23 | try: |
---|
24 | from typing import List |
---|
25 | except ImportError: |
---|
26 | pass |
---|
27 | # pylint: enable=unused-import |
---|
28 | |
---|
29 | def make_mixture_info(parts, operation='+'): |
---|
30 | # type: (List[ModelInfo]) -> ModelInfo |
---|
31 | """ |
---|
32 | Create info block for mixture model. |
---|
33 | """ |
---|
34 | # Build new parameter list |
---|
35 | combined_pars = [] |
---|
36 | |
---|
37 | all_parts = copy(parts) |
---|
38 | is_flat = False |
---|
39 | while not is_flat: |
---|
40 | is_flat = True |
---|
41 | for part in all_parts: |
---|
42 | if part.composition and part.composition[0] == 'mixture' and \ |
---|
43 | len(part.composition[1]) > 1: |
---|
44 | all_parts += part.composition[1] |
---|
45 | all_parts.remove(part) |
---|
46 | is_flat = False |
---|
47 | |
---|
48 | # When creating a mixture model that is a sum of product models (ie (1*2)+(3*4)) |
---|
49 | # the parameters for models 1 & 2 will be prefixed with A & B respectively, |
---|
50 | # but so will the parameters for models 3 & 4. We need to rename models 3 & 4 |
---|
51 | # so that they are prefixed with C & D to avoid overlap of parameter names. |
---|
52 | used_prefixes = [] |
---|
53 | for part in parts: |
---|
54 | i = 0 |
---|
55 | if part.composition and part.composition[0] == 'mixture': |
---|
56 | npars_list = [info.parameters.npars for info in part.composition[1]] |
---|
57 | for npars in npars_list: |
---|
58 | # List of params of one of the constituent models of part |
---|
59 | submodel_pars = part.parameters.kernel_parameters[i:i+npars] |
---|
60 | # Prefix of the constituent model |
---|
61 | prefix = submodel_pars[0].name[0] |
---|
62 | if prefix not in used_prefixes: # Haven't seen this prefix so far |
---|
63 | used_prefixes.append(prefix) |
---|
64 | i += npars |
---|
65 | continue |
---|
66 | while prefix in used_prefixes: |
---|
67 | # This prefix has been already used, so change it to the |
---|
68 | # next letter that hasn't been used |
---|
69 | prefix = chr(ord(prefix) + 1) |
---|
70 | used_prefixes.append(prefix) |
---|
71 | prefix += "_" |
---|
72 | # Update the parameters of this constituent model to use the |
---|
73 | # new prefix |
---|
74 | for par in submodel_pars: |
---|
75 | par.id = prefix + par.id[2:] |
---|
76 | par.name = prefix + par.name[2:] |
---|
77 | if par.length_control is not None: |
---|
78 | par.length_control = prefix + par.length_control[2:] |
---|
79 | i += npars |
---|
80 | |
---|
81 | for part in parts: |
---|
82 | # Parameter prefix per model, A_, B_, ... |
---|
83 | # Note that prefix must also be applied to id and length_control |
---|
84 | # to support vector parameters |
---|
85 | prefix = '' |
---|
86 | if not part.composition: |
---|
87 | # Model isn't a composition model, so it's parameters don't have a |
---|
88 | # a prefix. Add the next available prefix |
---|
89 | prefix = chr(ord('A')+len(used_prefixes)) |
---|
90 | used_prefixes.append(prefix) |
---|
91 | prefix += '_' |
---|
92 | |
---|
93 | if operation == '+': |
---|
94 | # If model is a sum model, each constituent model gets its own scale parameter |
---|
95 | scale_prefix = prefix |
---|
96 | if prefix == '' and getattr(part, "operation", '') == '*': |
---|
97 | # `part` is a composition product model. Find the prefixes of |
---|
98 | # it's parameters to form a new prefix for the scale. |
---|
99 | # For example, a model with A*B*C will have ABC_scale. |
---|
100 | sub_prefixes = [] |
---|
101 | for param in part.parameters.kernel_parameters: |
---|
102 | # Prefix of constituent model |
---|
103 | sub_prefix = param.id.split('_')[0] |
---|
104 | if sub_prefix not in sub_prefixes: |
---|
105 | sub_prefixes.append(sub_prefix) |
---|
106 | # Concatenate sub_prefixes to form prefix for the scale |
---|
107 | scale_prefix = ''.join(sub_prefixes) + '_' |
---|
108 | scale = Parameter(scale_prefix + 'scale', default=1.0, |
---|
109 | description="model intensity for " + part.name) |
---|
110 | combined_pars.append(scale) |
---|
111 | for p in part.parameters.kernel_parameters: |
---|
112 | p = copy(p) |
---|
113 | p.name = prefix + p.name |
---|
114 | p.id = prefix + p.id |
---|
115 | if p.length_control is not None: |
---|
116 | p.length_control = prefix + p.length_control |
---|
117 | combined_pars.append(p) |
---|
118 | parameters = ParameterTable(combined_pars) |
---|
119 | parameters.max_pd = sum(part.parameters.max_pd for part in parts) |
---|
120 | |
---|
121 | def random(): |
---|
122 | combined_pars = {} |
---|
123 | for k, part in enumerate(parts): |
---|
124 | prefix = chr(ord('A')+k) + '_' |
---|
125 | pars = part.random() |
---|
126 | combined_pars.update((prefix+k, v) for k, v in pars.items()) |
---|
127 | return combined_pars |
---|
128 | |
---|
129 | model_info = ModelInfo() |
---|
130 | model_info.id = operation.join(part.id for part in parts) |
---|
131 | model_info.operation = operation |
---|
132 | model_info.name = '(' + operation.join(part.name for part in parts) + ')' |
---|
133 | model_info.filename = None |
---|
134 | model_info.title = 'Mixture model with ' + model_info.name |
---|
135 | model_info.description = model_info.title |
---|
136 | model_info.docs = model_info.title |
---|
137 | model_info.category = "custom" |
---|
138 | model_info.parameters = parameters |
---|
139 | model_info.random = random |
---|
140 | #model_info.single = any(part['single'] for part in parts) |
---|
141 | model_info.structure_factor = False |
---|
142 | model_info.variant_info = None |
---|
143 | #model_info.tests = [] |
---|
144 | #model_info.source = [] |
---|
145 | # Remember the component info blocks so we can build the model |
---|
146 | model_info.composition = ('mixture', parts) |
---|
147 | return model_info |
---|
148 | |
---|
149 | |
---|
150 | class MixtureModel(KernelModel): |
---|
151 | def __init__(self, model_info, parts): |
---|
152 | # type: (ModelInfo, List[KernelModel]) -> None |
---|
153 | self.info = model_info |
---|
154 | self.parts = parts |
---|
155 | self.dtype = parts[0].dtype |
---|
156 | |
---|
157 | def make_kernel(self, q_vectors): |
---|
158 | # type: (List[np.ndarray]) -> MixtureKernel |
---|
159 | # Note: may be sending the q_vectors to the n times even though they |
---|
160 | # are only needed once. It would mess up modularity quite a bit to |
---|
161 | # handle this optimally, especially since there are many cases where |
---|
162 | # separate q vectors are needed (e.g., form in python and structure |
---|
163 | # in opencl; or both in opencl, but one in single precision and the |
---|
164 | # other in double precision). |
---|
165 | kernels = [part.make_kernel(q_vectors) for part in self.parts] |
---|
166 | return MixtureKernel(self.info, kernels) |
---|
167 | |
---|
168 | def release(self): |
---|
169 | # type: () -> None |
---|
170 | """ |
---|
171 | Free resources associated with the model. |
---|
172 | """ |
---|
173 | for part in self.parts: |
---|
174 | part.release() |
---|
175 | |
---|
176 | |
---|
177 | class MixtureKernel(Kernel): |
---|
178 | def __init__(self, model_info, kernels): |
---|
179 | # type: (ModelInfo, List[Kernel]) -> None |
---|
180 | self.dim = kernels[0].dim |
---|
181 | self.info = model_info |
---|
182 | self.kernels = kernels |
---|
183 | self.dtype = self.kernels[0].dtype |
---|
184 | self.operation = model_info.operation |
---|
185 | self.results = [] # type: List[np.ndarray] |
---|
186 | |
---|
187 | def __call__(self, call_details, values, cutoff, magnetic): |
---|
188 | # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray |
---|
189 | scale, background = values[0:2] |
---|
190 | total = 0.0 |
---|
191 | # remember the parts for plotting later |
---|
192 | self.results = [] # type: List[np.ndarray] |
---|
193 | parts = MixtureParts(self.info, self.kernels, call_details, values) |
---|
194 | for kernel, kernel_details, kernel_values in parts: |
---|
195 | #print("calling kernel", kernel.info.name) |
---|
196 | result = kernel(kernel_details, kernel_values, cutoff, magnetic) |
---|
197 | result = np.array(result).astype(kernel.dtype) |
---|
198 | # print(kernel.info.name, result) |
---|
199 | if self.operation == '+': |
---|
200 | total += result |
---|
201 | elif self.operation == '*': |
---|
202 | if np.all(total) == 0.0: |
---|
203 | total = result |
---|
204 | else: |
---|
205 | total *= result |
---|
206 | self.results.append(result) |
---|
207 | |
---|
208 | return scale*total + background |
---|
209 | |
---|
210 | def release(self): |
---|
211 | # type: () -> None |
---|
212 | for k in self.kernels: |
---|
213 | k.release() |
---|
214 | |
---|
215 | |
---|
216 | class MixtureParts(object): |
---|
217 | def __init__(self, model_info, kernels, call_details, values): |
---|
218 | # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None |
---|
219 | self.model_info = model_info |
---|
220 | self.parts = model_info.composition[1] |
---|
221 | self.kernels = kernels |
---|
222 | self.call_details = call_details |
---|
223 | self.values = values |
---|
224 | self.spin_index = model_info.parameters.npars + 2 |
---|
225 | #call_details.show(values) |
---|
226 | |
---|
227 | def __iter__(self): |
---|
228 | # type: () -> PartIterable |
---|
229 | self.part_num = 0 |
---|
230 | self.par_index = 2 |
---|
231 | self.mag_index = self.spin_index + 3 |
---|
232 | return self |
---|
233 | |
---|
234 | def __next__(self): |
---|
235 | # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] |
---|
236 | if self.part_num >= len(self.parts): |
---|
237 | raise StopIteration() |
---|
238 | info = self.parts[self.part_num] |
---|
239 | kernel = self.kernels[self.part_num] |
---|
240 | call_details = self._part_details(info, self.par_index) |
---|
241 | values = self._part_values(info, self.par_index, self.mag_index) |
---|
242 | values = values.astype(kernel.dtype) |
---|
243 | #call_details.show(values) |
---|
244 | |
---|
245 | self.part_num += 1 |
---|
246 | self.par_index += info.parameters.npars |
---|
247 | if self.model_info.operation == '+': |
---|
248 | self.par_index += 1 # Account for each constituent model's scale param |
---|
249 | self.mag_index += 3 * len(info.parameters.magnetism_index) |
---|
250 | |
---|
251 | return kernel, call_details, values |
---|
252 | |
---|
253 | # CRUFT: py2 support |
---|
254 | next = __next__ |
---|
255 | |
---|
256 | def _part_details(self, info, par_index): |
---|
257 | # type: (ModelInfo, int) -> CallDetails |
---|
258 | full = self.call_details |
---|
259 | # par_index is index into values array of the current parameter, |
---|
260 | # which includes the initial scale and background parameters. |
---|
261 | # We want the index into the weight length/offset for each parameter. |
---|
262 | # Exclude the initial scale and background, so subtract two. If we're |
---|
263 | # building an addition model, each component has its own scale factor |
---|
264 | # which we need to skip when constructing the details for the kernel, so |
---|
265 | # add one, giving a net subtract one. |
---|
266 | diff = 1 if self.model_info.operation == '+' else 2 |
---|
267 | index = slice(par_index - diff, par_index - diff + info.parameters.npars) |
---|
268 | length = full.length[index] |
---|
269 | offset = full.offset[index] |
---|
270 | # The complete weight vector is being sent to each part so that |
---|
271 | # offsets don't need to be adjusted. |
---|
272 | part = make_details(info, length, offset, full.num_weights) |
---|
273 | return part |
---|
274 | |
---|
275 | def _part_values(self, info, par_index, mag_index): |
---|
276 | # type: (ModelInfo, int, int) -> np.ndarray |
---|
277 | # Set each constituent model's scale to 1 if this is a multiplication model |
---|
278 | scale = self.values[par_index] if self.model_info.operation == '+' else 1.0 |
---|
279 | diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model |
---|
280 | pars = self.values[par_index + diff:par_index + info.parameters.npars + diff] |
---|
281 | nmagnetic = len(info.parameters.magnetism_index) |
---|
282 | if nmagnetic: |
---|
283 | spin_state = self.values[self.spin_index:self.spin_index + 3] |
---|
284 | mag_index = self.values[mag_index:mag_index + 3 * nmagnetic] |
---|
285 | else: |
---|
286 | spin_state = [] |
---|
287 | mag_index = [] |
---|
288 | nvalues = self.model_info.parameters.nvalues |
---|
289 | nweights = self.call_details.num_weights |
---|
290 | weights = self.values[nvalues:nvalues+2*nweights] |
---|
291 | zero = self.values.dtype.type(0.) |
---|
292 | values = [[scale, zero], pars, spin_state, mag_index, weights] |
---|
293 | # Pad value array to a 32 value boundary |
---|
294 | spacer = (32 - sum(len(v) for v in values)%32)%32 |
---|
295 | values.append([zero]*spacer) |
---|
296 | values = np.hstack(values).astype(self.kernels[0].dtype) |
---|
297 | return values |
---|