1 | """ |
---|
2 | Product 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 | *make_product_info(P, S)*. |
---|
12 | """ |
---|
13 | from __future__ import print_function, division |
---|
14 | |
---|
15 | from collections import OrderedDict |
---|
16 | |
---|
17 | from copy import copy |
---|
18 | import numpy as np # type: ignore |
---|
19 | |
---|
20 | from .modelinfo import ParameterTable, ModelInfo, parse_parameter |
---|
21 | from .kernel import KernelModel, Kernel |
---|
22 | from .details import make_details |
---|
23 | |
---|
24 | # pylint: disable=unused-import |
---|
25 | try: |
---|
26 | from typing import Tuple, Callable, Union |
---|
27 | except ImportError: |
---|
28 | pass |
---|
29 | else: |
---|
30 | from .modelinfo import ParameterSet, Parameter |
---|
31 | # pylint: enable=unused-import |
---|
32 | |
---|
33 | # TODO: make shape averages available to constraints |
---|
34 | #ESTIMATED_PARAMETERS = [ |
---|
35 | # ["mean_radius_effective", "A", 0.0, [0, np.inf], "", "mean effective radius"], |
---|
36 | # ["mean_volume", "A", 0.0, [0, np.inf], "", "mean volume"], |
---|
37 | # ["mean_volume_ratio", "", 1.0, [0, np.inf], "", "mean form: mean shell volume ratio"], |
---|
38 | #] |
---|
39 | STRUCTURE_MODE_ID = "structure_factor_mode" |
---|
40 | RADIUS_MODE_ID = "radius_effective_mode" |
---|
41 | RADIUS_ID = "radius_effective" |
---|
42 | VOLFRAC_ID = "volfraction" |
---|
43 | def make_extra_pars(p_info): |
---|
44 | # type: (ModelInfo) -> List[Parameter] |
---|
45 | """ |
---|
46 | Create parameters for structure factor and effective radius modes. |
---|
47 | """ |
---|
48 | pars = [] |
---|
49 | if p_info.have_Fq: |
---|
50 | par = parse_parameter( |
---|
51 | STRUCTURE_MODE_ID, |
---|
52 | "", |
---|
53 | 0, |
---|
54 | [["P*S", "P*(1+beta*(S-1))"]], |
---|
55 | "", |
---|
56 | "Structure factor calculation") |
---|
57 | pars.append(par) |
---|
58 | if p_info.radius_effective_modes is not None: |
---|
59 | par = parse_parameter( |
---|
60 | RADIUS_MODE_ID, |
---|
61 | "", |
---|
62 | 1, |
---|
63 | [["unconstrained"] + p_info.radius_effective_modes], |
---|
64 | "", |
---|
65 | "Effective radius calculation") |
---|
66 | pars.append(par) |
---|
67 | return pars |
---|
68 | |
---|
69 | def make_product_info(p_info, s_info): |
---|
70 | # type: (ModelInfo, ModelInfo) -> ModelInfo |
---|
71 | """ |
---|
72 | Create info block for product model. |
---|
73 | """ |
---|
74 | # Make sure effective radius is the first parameter and |
---|
75 | # make sure volume fraction is the second parameter of the |
---|
76 | # structure factor calculator. Structure factors should not |
---|
77 | # have any magnetic parameters |
---|
78 | if not len(s_info.parameters.kernel_parameters) >= 2: |
---|
79 | raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID)) |
---|
80 | if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID: |
---|
81 | raise TypeError("S needs {} as first parameter".format(RADIUS_ID)) |
---|
82 | if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID: |
---|
83 | raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID)) |
---|
84 | if not s_info.parameters.magnetism_index == []: |
---|
85 | raise TypeError("S should not have SLD parameters") |
---|
86 | p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters |
---|
87 | s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters |
---|
88 | |
---|
89 | # Create list of parameters for the combined model. If there |
---|
90 | # are any names in P that overlap with those in S, modify the name in S |
---|
91 | # to distinguish it. |
---|
92 | p_set = set(p.id for p in p_pars.kernel_parameters) |
---|
93 | s_list = [(_tag_parameter(par) if par.id in p_set else par) |
---|
94 | for par in s_pars.kernel_parameters] |
---|
95 | # Check if still a collision after renaming. This could happen if for |
---|
96 | # example S has volfrac and P has both volfrac and volfrac_S. |
---|
97 | if any(p.id in p_set for p in s_list): |
---|
98 | raise TypeError("name collision: P has P.name and P.name_S while S has S.name") |
---|
99 | |
---|
100 | # make sure effective radius is not a polydisperse parameter in product |
---|
101 | s_list[0] = copy(s_list[0]) |
---|
102 | s_list[0].polydisperse = False |
---|
103 | |
---|
104 | translate_name = dict((old.id, new.id) for old, new |
---|
105 | in zip(s_pars.kernel_parameters, s_list)) |
---|
106 | combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) |
---|
107 | parameters = ParameterTable(combined_pars) |
---|
108 | parameters.max_pd = p_pars.max_pd + s_pars.max_pd |
---|
109 | def random(): |
---|
110 | """Random set of model parameters for product model""" |
---|
111 | combined_pars = p_info.random() |
---|
112 | s_names = set(par.id for par in s_pars.kernel_parameters) |
---|
113 | combined_pars.update((translate_name[k], v) |
---|
114 | for k, v in s_info.random().items() |
---|
115 | if k in s_names) |
---|
116 | return combined_pars |
---|
117 | |
---|
118 | model_info = ModelInfo() |
---|
119 | model_info.id = '@'.join((p_id, s_id)) |
---|
120 | model_info.name = '@'.join((p_name, s_name)) |
---|
121 | model_info.filename = None |
---|
122 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
---|
123 | model_info.description = model_info.title |
---|
124 | model_info.docs = model_info.title |
---|
125 | model_info.category = "custom" |
---|
126 | model_info.parameters = parameters |
---|
127 | model_info.random = random |
---|
128 | #model_info.single = p_info.single and s_info.single |
---|
129 | model_info.structure_factor = False |
---|
130 | model_info.variant_info = None |
---|
131 | #model_info.tests = [] |
---|
132 | #model_info.source = [] |
---|
133 | # Remember the component info blocks so we can build the model |
---|
134 | model_info.composition = ('product', [p_info, s_info]) |
---|
135 | model_info.hidden = p_info.hidden |
---|
136 | if getattr(p_info, 'profile', None) is not None: |
---|
137 | profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) |
---|
138 | def profile(**kwargs): |
---|
139 | """Return SLD profile of the form factor as a function of radius.""" |
---|
140 | # extract the profile args |
---|
141 | kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars) |
---|
142 | return p_info.profile(**kwargs) |
---|
143 | else: |
---|
144 | profile = None |
---|
145 | model_info.profile = profile |
---|
146 | model_info.profile_axes = p_info.profile_axes |
---|
147 | |
---|
148 | # TODO: delegate random to p_info, s_info |
---|
149 | #model_info.random = lambda: {} |
---|
150 | |
---|
151 | ## Show the parameter table |
---|
152 | #from .compare import get_pars, parlist |
---|
153 | #print("==== %s ====="%model_info.name) |
---|
154 | #values = get_pars(model_info) |
---|
155 | #print(parlist(model_info, values, is2d=True)) |
---|
156 | return model_info |
---|
157 | |
---|
158 | def _tag_parameter(par): |
---|
159 | """ |
---|
160 | Tag the parameter name with _S to indicate that the parameter comes from |
---|
161 | the structure factor parameter set. This is only necessary if the |
---|
162 | form factor model includes a parameter of the same name as a parameter |
---|
163 | in the structure factor. |
---|
164 | """ |
---|
165 | par = copy(par) |
---|
166 | # Protect against a vector parameter in S by appending the vector length |
---|
167 | # to the renamed parameter. Note: haven't tested this since no existing |
---|
168 | # structure factor models contain vector parameters. |
---|
169 | vector_length = par.name[len(par.id):] |
---|
170 | par.id = par.id + "_S" |
---|
171 | par.name = par.id + vector_length |
---|
172 | return par |
---|
173 | |
---|
174 | def _intermediates( |
---|
175 | F1, # type: np.ndarray |
---|
176 | F2, # type: np.ndarray |
---|
177 | S, # type: np.ndarray |
---|
178 | scale, # type: float |
---|
179 | radius_effective, # type: float |
---|
180 | beta_mode, # type: bool |
---|
181 | ): |
---|
182 | # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] |
---|
183 | """ |
---|
184 | Returns intermediate results for beta approximation-enabled product. |
---|
185 | The result may be an array or a float. |
---|
186 | """ |
---|
187 | # CRUFT: remove effective_radius once SasView 5.0 is released. |
---|
188 | if beta_mode: |
---|
189 | # TODO: 1. include calculated Q vector |
---|
190 | # TODO: 2. consider implications if there are intermediate results in P(Q) |
---|
191 | parts = OrderedDict(( |
---|
192 | ("P(Q)", scale*F2), |
---|
193 | ("S(Q)", S), |
---|
194 | ("beta(Q)", F1**2 / F2), |
---|
195 | ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), |
---|
196 | ("effective_radius", radius_effective), |
---|
197 | ("radius_effective", radius_effective), |
---|
198 | # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), |
---|
199 | )) |
---|
200 | else: |
---|
201 | parts = OrderedDict(( |
---|
202 | ("P(Q)", scale*F2), |
---|
203 | ("S(Q)", S), |
---|
204 | ("effective_radius", radius_effective), |
---|
205 | ("radius_effective", radius_effective), |
---|
206 | )) |
---|
207 | return parts |
---|
208 | |
---|
209 | class ProductModel(KernelModel): |
---|
210 | """ |
---|
211 | Model definition for product model. |
---|
212 | """ |
---|
213 | def __init__(self, model_info, P, S): |
---|
214 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
---|
215 | #: Combined info plock for the product model |
---|
216 | self.info = model_info |
---|
217 | #: Form factor modelling individual particles. |
---|
218 | self.P = P |
---|
219 | #: Structure factor modelling interaction between particles. |
---|
220 | self.S = S |
---|
221 | |
---|
222 | #: Model precision. This is not really relevant, since it is the |
---|
223 | #: individual P and S models that control the effective dtype, |
---|
224 | #: converting the q-vectors to the correct type when the kernels |
---|
225 | #: for each are created. Ideally this should be set to the more |
---|
226 | #: precise type to avoid loss of precision, but precision in q is |
---|
227 | #: not critical (single is good enough for our purposes), so it just |
---|
228 | #: uses the precision of the form factor. |
---|
229 | self.dtype = P.dtype # type: np.dtype |
---|
230 | |
---|
231 | def make_kernel(self, q_vectors): |
---|
232 | # type: (List[np.ndarray]) -> Kernel |
---|
233 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
234 | # are only needed once. It would mess up modularity quite a bit to |
---|
235 | # handle this optimally, especially since there are many cases where |
---|
236 | # separate q vectors are needed (e.g., form in python and structure |
---|
237 | # in opencl; or both in opencl, but one in single precision and the |
---|
238 | # other in double precision). |
---|
239 | |
---|
240 | p_kernel = self.P.make_kernel(q_vectors) |
---|
241 | s_kernel = self.S.make_kernel(q_vectors) |
---|
242 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
243 | make_kernel.__doc__ = KernelModel.make_kernel.__doc__ |
---|
244 | |
---|
245 | def release(self): |
---|
246 | # type: (None) -> None |
---|
247 | """ |
---|
248 | Free resources associated with the model. |
---|
249 | """ |
---|
250 | self.P.release() |
---|
251 | self.S.release() |
---|
252 | |
---|
253 | |
---|
254 | class ProductKernel(Kernel): |
---|
255 | """ |
---|
256 | Instantiated kernel for product model. |
---|
257 | """ |
---|
258 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
259 | # type: (ModelInfo, Kernel, Kernel) -> None |
---|
260 | self.info = model_info |
---|
261 | self.p_kernel = p_kernel |
---|
262 | self.s_kernel = s_kernel |
---|
263 | self.dtype = p_kernel.dtype |
---|
264 | self.results = [] # type: List[np.ndarray] |
---|
265 | |
---|
266 | def Iq(self, call_details, values, cutoff, magnetic): |
---|
267 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
---|
268 | |
---|
269 | p_info, s_info = self.info.composition[1] |
---|
270 | p_npars = p_info.parameters.npars |
---|
271 | p_length = call_details.length[:p_npars] |
---|
272 | p_offset = call_details.offset[:p_npars] |
---|
273 | s_npars = s_info.parameters.npars |
---|
274 | s_length = call_details.length[p_npars:p_npars+s_npars] |
---|
275 | s_offset = call_details.offset[p_npars:p_npars+s_npars] |
---|
276 | |
---|
277 | # Beta mode parameter is the first parameter after P and S parameters |
---|
278 | have_beta_mode = p_info.have_Fq |
---|
279 | beta_mode_offset = 2+p_npars+s_npars |
---|
280 | beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False |
---|
281 | if beta_mode and self.p_kernel.dim == '2d': |
---|
282 | raise NotImplementedError("beta not yet supported for 2D") |
---|
283 | |
---|
284 | # R_eff type parameter is the second parameter after P and S parameters |
---|
285 | # unless the model doesn't support beta mode, in which case it is first |
---|
286 | have_radius_type = p_info.radius_effective_modes is not None |
---|
287 | #print(p_npars,s_npars) |
---|
288 | radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0) |
---|
289 | #print(values[radius_type_offset]) |
---|
290 | radius_type = int(values[radius_type_offset]) if have_radius_type else 0 |
---|
291 | |
---|
292 | # Retrieve the volume fraction, which is the second of the |
---|
293 | # 'S' parameters in the parameter list, or 2+np in 0-origin, |
---|
294 | # as well as the scale and background. |
---|
295 | volfrac = values[3+p_npars] |
---|
296 | scale, background = values[0], values[1] |
---|
297 | |
---|
298 | # if there are magnetic parameters, they will only be on the |
---|
299 | # form factor P, not the structure factor S. |
---|
300 | nmagnetic = len(self.info.parameters.magnetism_index) |
---|
301 | if nmagnetic: |
---|
302 | spin_index = self.info.parameters.npars + 2 |
---|
303 | magnetism = values[spin_index: spin_index+3+3*nmagnetic] |
---|
304 | else: |
---|
305 | magnetism = [] |
---|
306 | nvalues = self.info.parameters.nvalues |
---|
307 | nweights = call_details.num_weights |
---|
308 | weights = values[nvalues:nvalues + 2*nweights] |
---|
309 | |
---|
310 | # Construct the calling parameters for P. |
---|
311 | p_details = make_details(p_info, p_length, p_offset, nweights) |
---|
312 | p_values = [ |
---|
313 | [1., 0.], # scale=1, background=0, |
---|
314 | values[2:2+p_npars], |
---|
315 | magnetism, |
---|
316 | weights] |
---|
317 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
---|
318 | p_values.append([0.]*spacer) |
---|
319 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
---|
320 | |
---|
321 | # Construct the calling parameters for S. |
---|
322 | if radius_type > 0: |
---|
323 | # If R_eff comes from form factor, make sure it is monodisperse. |
---|
324 | # weight is set to 1 later, after the value array is created |
---|
325 | s_length[0] = 1 |
---|
326 | s_details = make_details(s_info, s_length, s_offset, nweights) |
---|
327 | s_values = [ |
---|
328 | [1., 0.], # scale=1, background=0, |
---|
329 | values[2+p_npars:2+p_npars+s_npars], |
---|
330 | weights, |
---|
331 | ] |
---|
332 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
---|
333 | s_values.append([0.]*spacer) |
---|
334 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
---|
335 | |
---|
336 | # Call the form factor kernel to compute <F> and <F^2>. |
---|
337 | # If the model doesn't support Fq the returned <F> will be None. |
---|
338 | F1, F2, radius_effective, shell_volume, volume_ratio = self.p_kernel.Fq( |
---|
339 | p_details, p_values, cutoff, magnetic, radius_type) |
---|
340 | |
---|
341 | # Call the structure factor kernel to compute S. |
---|
342 | # Plug R_eff from the form factor into structure factor parameters |
---|
343 | # and scale volume fraction by form:shell volume ratio. These changes |
---|
344 | # needs to be both in the initial value slot as well as the |
---|
345 | # polydispersity distribution slot in the values array due to |
---|
346 | # implementation details in kernel_iq.c. |
---|
347 | #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g" |
---|
348 | # % (radius_type, radius_effective, volfrac, volume_ratio)) |
---|
349 | if radius_type > 0: |
---|
350 | # set the value to the model R_eff and set the weight to 1 |
---|
351 | s_values[2] = s_values[2+s_npars+s_offset[0]] = radius_effective |
---|
352 | s_values[2+s_npars+s_offset[0]+nweights] = 1.0 |
---|
353 | s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio |
---|
354 | S = self.s_kernel.Iq(s_details, s_values, cutoff, False) |
---|
355 | |
---|
356 | # Determine overall scale factor. Hollow shapes are weighted by |
---|
357 | # shell_volume, so that is needed for volume normalization. For |
---|
358 | # solid shapes we can use shell_volume as well since it is equal |
---|
359 | # to form volume. |
---|
360 | combined_scale = scale*volfrac/shell_volume |
---|
361 | |
---|
362 | # Combine form factor and structure factor |
---|
363 | #print("beta", beta_mode, F1, F2, S) |
---|
364 | PS = F2 + F1**2*(S-1) if beta_mode else F2*S |
---|
365 | final_result = combined_scale*PS + background |
---|
366 | |
---|
367 | # Capture intermediate values so user can see them. These are |
---|
368 | # returned as a lazy evaluator since they are only needed in the |
---|
369 | # GUI, and not for each evaluation during a fit. |
---|
370 | # TODO: return the results structure with the final results |
---|
371 | # That way the model calcs are idempotent. Further, we can |
---|
372 | # generalize intermediates to various other model types if we put it |
---|
373 | # kernel calling interface. Could do this as an "optional" |
---|
374 | # return value in the caller, though in that case we could return |
---|
375 | # the results directly rather than through a lazy evaluator. |
---|
376 | self.results = lambda: _intermediates( |
---|
377 | F1, F2, S, combined_scale, radius_effective, beta_mode) |
---|
378 | |
---|
379 | return final_result |
---|
380 | |
---|
381 | Iq.__doc__ = Kernel.Iq.__doc__ |
---|
382 | __call__ = Iq |
---|
383 | |
---|
384 | def release(self): |
---|
385 | # type: () -> None |
---|
386 | """Free resources associated with the kernel.""" |
---|
387 | self.p_kernel.release() |
---|
388 | self.s_kernel.release() |
---|