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