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 | The P@S models is somewhat complicated because there are many special |
---|
14 | parameters that need to be handled in particular ways. Much of the |
---|
15 | code is used to figure out what special parameters we have, where to |
---|
16 | find them in the P@S model inputs and how to distribute them to the underlying |
---|
17 | P and S model calculators. |
---|
18 | |
---|
19 | The parameter packet received by the P@S is a :class:`details.CallDetails` |
---|
20 | structure along with a data vector. The CallDetails structure indicates which |
---|
21 | parameters are polydisperse, the length of the distribution, and where to |
---|
22 | find it in the data vector. The distributions are ordered from longest to |
---|
23 | shortest, with length 1 distributions filling out the distribution set. That |
---|
24 | way the kernel calculator doesn't have to check if it needs another nesting |
---|
25 | level since it is always there. The data vector consists of a list of target |
---|
26 | values for the parameters, followed by a concatenation of the distribution |
---|
27 | values, and then followed by a concatenation of the distribution weights. |
---|
28 | Given the combined details and data for P@S, we must decompose them in to |
---|
29 | details for P and details for S separately, which unfortunately requires |
---|
30 | intimate knowledge of the data structures and tricky code. |
---|
31 | |
---|
32 | The special parameters are: |
---|
33 | |
---|
34 | * *scale* and *background*: |
---|
35 | First two parameters of the value list in each of P, S and P@S. |
---|
36 | When decomposing P@S parameters, ignore *scale* and *background*, |
---|
37 | instead using 1 and 0 for the first two slots of both P and S. |
---|
38 | After calling P and S individually, the results are combined as |
---|
39 | :code:`volfraction*scale*P*S + background`. The *scale* and |
---|
40 | *background* do not show up in the polydispersity structure so |
---|
41 | they are easy to handle. |
---|
42 | |
---|
43 | * *volfraction*: |
---|
44 | Always the first parameter of S, but it may also be in P. If it is in P, |
---|
45 | then *P.volfraction* is used in the combined P@S list, and |
---|
46 | *S.volfraction* is elided, otherwise *S.volfraction* is used. If we |
---|
47 | are using *volfraction* from P we can treat it like all the other P |
---|
48 | parameters when calling P, but when calling S we need to insert the |
---|
49 | *P.volfraction* into data vector for S and assign a slot of length 1 |
---|
50 | in the distribution. Because we are using the original layout of the |
---|
51 | distribution vectors from P@S, but copying it into private data |
---|
52 | vectors for S and P, we are free to "borrow" a P slots to store the |
---|
53 | missing *S.volfraction* distribution. We use the *P.volfraction* |
---|
54 | slot itself but any slot will work. |
---|
55 | |
---|
56 | For hollow shapes, *volfraction* represents the volume fraction of |
---|
57 | material but S needs the volume fraction enclosed by the shape. The |
---|
58 | answer is to scale the user specified volume fraction by the form:shell |
---|
59 | ratio computed from the average form volume and average shell volume |
---|
60 | returned from P. Use the original *volfraction* divided by *shell_volume* |
---|
61 | to compute the number density, and scale P@S by that to get absolute |
---|
62 | scaling on the final *I(q)*. The *scale* for P@S should therefore usually |
---|
63 | be one. |
---|
64 | |
---|
65 | * *radius_effective*: |
---|
66 | Always the second parameter of S and always part of P@S, but never in P. |
---|
67 | The value may be calculated using *P.radius_effective()* or it |
---|
68 | may be set to the *radius_effective* value in P@S, depending on |
---|
69 | *radius_effective_mode*. If part of S, the value may be polydisperse. |
---|
70 | If calculated by P, then it will be the weighted average of effective |
---|
71 | radii computed for the polydisperse shape parameters. |
---|
72 | |
---|
73 | * *structure_factor_mode* |
---|
74 | If P@S supports beta approximation (i.e., if it has the *Fq* function that |
---|
75 | returns <FF*> and <F><F*>), then *structure_factor_mode* will be added |
---|
76 | to the P@S parameters right after the S parameters. This mode may be 0 |
---|
77 | for the monodisperse approximation or 1 for the beta approximation. We |
---|
78 | will add more values here as we implemented more complicated operations, |
---|
79 | but for now P and S must be computed separately. If beta, then we |
---|
80 | return *I = scale volfrac/volume ( <FF> + <F>^2 (S-1)) + background*. |
---|
81 | If not beta then return *I = scale/volume P S + background* . In both |
---|
82 | cases, return the appropriate immediate values. |
---|
83 | |
---|
84 | * *radius_effective_mode* |
---|
85 | If P defines the *radius_effective* function (and therefore |
---|
86 | *P.info.radius_effective_modes* is a list of effective radius modes), |
---|
87 | then *radius_effective_mode* will be the final parameter in P@S. Mode |
---|
88 | will be zero if *radius_effective* is defined by the user using the S |
---|
89 | parameter; any other value and the *radius_effective* parameter will be |
---|
90 | filled in from the value computed in P. In the latter case, the |
---|
91 | polydispersity information for *S.radius_effective* will need to be |
---|
92 | suppressed, with pd length set to 1, the first value set to the |
---|
93 | effective radius and the first weight set to 1. Do this after composing |
---|
94 | the S data vector so the inputs are left untouched. |
---|
95 | |
---|
96 | * *regular parameters* |
---|
97 | The regular P parameters form a block of length *P.info.npars* at the |
---|
98 | start of the data vector (after scale and background). These will be |
---|
99 | followed by *S.effective_radius*, and *S.volfraction* (if *P.volfraction* |
---|
100 | is absent), and then the regular S parameters. The P and S blocks can |
---|
101 | be copied as a group into the respective P and S data vectors. |
---|
102 | We can copy the distribution value and weight vectors untouched to both |
---|
103 | the P and S data vectors since they are referenced by offset and length. |
---|
104 | We can update the radius_effective slots in the P data vector with |
---|
105 | *P.radius_effective()* if needed. |
---|
106 | |
---|
107 | * *magnetic parameters* |
---|
108 | For each P parameter that is an SLD there will be a set of three magnetic |
---|
109 | parameters tacked on to P@S after the regular P and S and after the |
---|
110 | special *structure_factor_mode* and *radius_effective_mode*. These |
---|
111 | can be copied as a group after the regular P parameters. There won't |
---|
112 | be any magnetic S parameters. |
---|
113 | |
---|
114 | """ |
---|
115 | from __future__ import print_function, division |
---|
116 | |
---|
117 | from collections import OrderedDict |
---|
118 | |
---|
119 | from copy import copy |
---|
120 | import numpy as np # type: ignore |
---|
121 | |
---|
122 | from .modelinfo import ParameterTable, ModelInfo, parse_parameter |
---|
123 | from .kernel import KernelModel, Kernel |
---|
124 | from .details import make_details |
---|
125 | |
---|
126 | # pylint: disable=unused-import |
---|
127 | try: |
---|
128 | from typing import Tuple, Callable, Union |
---|
129 | except ImportError: |
---|
130 | pass |
---|
131 | else: |
---|
132 | from .modelinfo import ParameterSet, Parameter |
---|
133 | # pylint: enable=unused-import |
---|
134 | |
---|
135 | # TODO: make shape averages available to constraints |
---|
136 | #ESTIMATED_PARAMETERS = [ |
---|
137 | # ["mean_radius_effective", "A", 0.0, [0, np.inf], "", "mean effective radius"], |
---|
138 | # ["mean_volume", "A", 0.0, [0, np.inf], "", "mean volume"], |
---|
139 | # ["mean_volume_ratio", "", 1.0, [0, np.inf], "", "mean form: mean shell volume ratio"], |
---|
140 | #] |
---|
141 | STRUCTURE_MODE_ID = "structure_factor_mode" |
---|
142 | RADIUS_MODE_ID = "radius_effective_mode" |
---|
143 | RADIUS_ID = "radius_effective" |
---|
144 | VOLFRAC_ID = "volfraction" |
---|
145 | def make_extra_pars(p_info): |
---|
146 | # type: (ModelInfo) -> List[Parameter] |
---|
147 | """ |
---|
148 | Create parameters for structure factor and effective radius modes. |
---|
149 | """ |
---|
150 | pars = [] |
---|
151 | if p_info.have_Fq: |
---|
152 | par = parse_parameter( |
---|
153 | STRUCTURE_MODE_ID, |
---|
154 | "", |
---|
155 | 0, |
---|
156 | [["P*S", "P*(1+beta*(S-1))"]], |
---|
157 | "", |
---|
158 | "Structure factor calculation") |
---|
159 | pars.append(par) |
---|
160 | if p_info.radius_effective_modes is not None: |
---|
161 | par = parse_parameter( |
---|
162 | RADIUS_MODE_ID, |
---|
163 | "", |
---|
164 | 1, |
---|
165 | [["unconstrained"] + p_info.radius_effective_modes], |
---|
166 | "", |
---|
167 | "Effective radius calculation") |
---|
168 | pars.append(par) |
---|
169 | return pars |
---|
170 | |
---|
171 | def make_product_info(p_info, s_info): |
---|
172 | # type: (ModelInfo, ModelInfo) -> ModelInfo |
---|
173 | """ |
---|
174 | Create info block for product model. |
---|
175 | """ |
---|
176 | # Make sure effective radius is the first parameter and |
---|
177 | # make sure volume fraction is the second parameter of the |
---|
178 | # structure factor calculator. Structure factors should not |
---|
179 | # have any magnetic parameters |
---|
180 | if not len(s_info.parameters.kernel_parameters) >= 2: |
---|
181 | raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID)) |
---|
182 | if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID: |
---|
183 | raise TypeError("S needs {} as first parameter".format(RADIUS_ID)) |
---|
184 | if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID: |
---|
185 | raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID)) |
---|
186 | if not s_info.parameters.magnetism_index == []: |
---|
187 | raise TypeError("S should not have SLD parameters") |
---|
188 | if RADIUS_ID in p_info.parameters: |
---|
189 | raise TypeError("P should not have {}".format(RADIUS_ID)) |
---|
190 | p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters |
---|
191 | s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters |
---|
192 | p_has_volfrac = VOLFRAC_ID in p_info.parameters |
---|
193 | |
---|
194 | # Create list of parameters for the combined model. If a name in |
---|
195 | # S overlaps a name in P, tag the S parameter name to distinguish it. |
---|
196 | # If the tagged name also collides it will be caught by the parameter |
---|
197 | # table builder. Similarly if any special names are abused. Need the |
---|
198 | # pairs to create the translation table for random model generation. |
---|
199 | p_set = set(p.id for p in p_pars.kernel_parameters) |
---|
200 | s_pairs = [(par, (_tag_parameter(par) if par.id in p_set else par)) |
---|
201 | for par in s_pars.kernel_parameters |
---|
202 | # Note: exclude volfraction from s_list if volfraction in p |
---|
203 | if par.id != VOLFRAC_ID or not p_has_volfrac] |
---|
204 | s_list = [pair[0] for pair in s_pairs] |
---|
205 | |
---|
206 | # Build combined parameter table |
---|
207 | combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) |
---|
208 | parameters = ParameterTable(combined_pars) |
---|
209 | # Allow for the scenario in which each component has all its PD parameters |
---|
210 | # active simultaneously. details.make_details() will throw an error if |
---|
211 | # too many are used from any one component. |
---|
212 | parameters.Pumax_pd = p_pars.max_pd + s_pars.max_pd |
---|
213 | |
---|
214 | # TODO: does user-defined polydisperse S.radius_effective make sense? |
---|
215 | # make sure effective radius is not a polydisperse parameter in product |
---|
216 | #s_list[0] = copy(s_list[0]) |
---|
217 | #s_list[0].polydisperse = False |
---|
218 | |
---|
219 | s_translate = {old.id: new.id for old, new in s_pairs} |
---|
220 | def random(): |
---|
221 | """Random set of model parameters for product model""" |
---|
222 | combined_pars = p_info.random() |
---|
223 | combined_pars.update((s_translate[k], v) |
---|
224 | for k, v in s_info.random().items() |
---|
225 | if k in s_translate) |
---|
226 | return combined_pars |
---|
227 | |
---|
228 | model_info = ModelInfo() |
---|
229 | model_info.id = '@'.join((p_id, s_id)) |
---|
230 | model_info.name = '@'.join((p_name, s_name)) |
---|
231 | model_info.filename = None |
---|
232 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
---|
233 | model_info.description = model_info.title |
---|
234 | model_info.docs = model_info.title |
---|
235 | model_info.category = "custom" |
---|
236 | model_info.parameters = parameters |
---|
237 | model_info.random = random |
---|
238 | #model_info.single = p_info.single and s_info.single |
---|
239 | model_info.structure_factor = False |
---|
240 | model_info.variant_info = None |
---|
241 | #model_info.tests = [] |
---|
242 | #model_info.source = [] |
---|
243 | # Remember the component info blocks so we can build the model |
---|
244 | model_info.composition = ('product', [p_info, s_info]) |
---|
245 | model_info.hidden = p_info.hidden |
---|
246 | if getattr(p_info, 'profile', None) is not None: |
---|
247 | profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) |
---|
248 | def profile(**kwargs): |
---|
249 | """Return SLD profile of the form factor as a function of radius.""" |
---|
250 | # extract the profile args |
---|
251 | kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars) |
---|
252 | return p_info.profile(**kwargs) |
---|
253 | else: |
---|
254 | profile = None |
---|
255 | model_info.profile = profile |
---|
256 | model_info.profile_axes = p_info.profile_axes |
---|
257 | |
---|
258 | # TODO: delegate random to p_info, s_info |
---|
259 | #model_info.random = lambda: {} |
---|
260 | |
---|
261 | ## Show the parameter table |
---|
262 | #from .compare import get_pars, parlist |
---|
263 | #print("==== %s ====="%model_info.name) |
---|
264 | #values = get_pars(model_info) |
---|
265 | #print(parlist(model_info, values, is2d=True)) |
---|
266 | return model_info |
---|
267 | |
---|
268 | def _tag_parameter(par): |
---|
269 | """ |
---|
270 | Tag the parameter name with _S to indicate that the parameter comes from |
---|
271 | the structure factor parameter set. This is only necessary if the |
---|
272 | form factor model includes a parameter of the same name as a parameter |
---|
273 | in the structure factor. |
---|
274 | """ |
---|
275 | par = copy(par) |
---|
276 | # Protect against a vector parameter in S by appending the vector length |
---|
277 | # to the renamed parameter. Note: haven't tested this since no existing |
---|
278 | # structure factor models contain vector parameters. |
---|
279 | vector_length = par.name[len(par.id):] |
---|
280 | par.id = par.id + "_S" |
---|
281 | par.name = par.id + vector_length |
---|
282 | return par |
---|
283 | |
---|
284 | def _intermediates( |
---|
285 | F, # type: np.ndarray |
---|
286 | Fsq, # type: np.ndarray |
---|
287 | S, # type: np.ndarray |
---|
288 | scale, # type: float |
---|
289 | radius_effective, # type: float |
---|
290 | beta_mode, # type: bool |
---|
291 | ): |
---|
292 | # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] |
---|
293 | """ |
---|
294 | Returns intermediate results for beta approximation-enabled product. |
---|
295 | The result may be an array or a float. |
---|
296 | """ |
---|
297 | if beta_mode: |
---|
298 | # TODO: 1. include calculated Q vector |
---|
299 | # TODO: 2. consider implications if there are intermediate results in P(Q) |
---|
300 | parts = OrderedDict(( |
---|
301 | ("P(Q)", scale*Fsq), |
---|
302 | ("S(Q)", S), |
---|
303 | ("beta(Q)", F**2 / Fsq), |
---|
304 | ("S_eff(Q)", 1 + (F**2 / Fsq)*(S-1)), |
---|
305 | ("effective_radius", radius_effective), |
---|
306 | # ("I(Q)", scale*(Fsq + (F**2)*(S-1)) + bg), |
---|
307 | )) |
---|
308 | else: |
---|
309 | parts = OrderedDict(( |
---|
310 | ("P(Q)", scale*Fsq), |
---|
311 | ("S(Q)", S), |
---|
312 | ("effective_radius", radius_effective), |
---|
313 | )) |
---|
314 | return parts |
---|
315 | |
---|
316 | class ProductModel(KernelModel): |
---|
317 | """ |
---|
318 | Model definition for product model. |
---|
319 | """ |
---|
320 | def __init__(self, model_info, P, S): |
---|
321 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
---|
322 | #: Combined info plock for the product model |
---|
323 | self.info = model_info |
---|
324 | #: Form factor modelling individual particles. |
---|
325 | self.P = P |
---|
326 | #: Structure factor modelling interaction between particles. |
---|
327 | self.S = S |
---|
328 | |
---|
329 | #: Model precision. This is not really relevant, since it is the |
---|
330 | #: individual P and S models that control the effective dtype, |
---|
331 | #: converting the q-vectors to the correct type when the kernels |
---|
332 | #: for each are created. Ideally this should be set to the more |
---|
333 | #: precise type to avoid loss of precision, but precision in q is |
---|
334 | #: not critical (single is good enough for our purposes), so it just |
---|
335 | #: uses the precision of the form factor. |
---|
336 | self.dtype = P.dtype # type: np.dtype |
---|
337 | |
---|
338 | def make_kernel(self, q_vectors): |
---|
339 | # type: (List[np.ndarray]) -> Kernel |
---|
340 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
341 | # are only needed once. It would mess up modularity quite a bit to |
---|
342 | # handle this optimally, especially since there are many cases where |
---|
343 | # separate q vectors are needed (e.g., form in python and structure |
---|
344 | # in opencl; or both in opencl, but one in single precision and the |
---|
345 | # other in double precision). |
---|
346 | |
---|
347 | p_kernel = self.P.make_kernel(q_vectors) |
---|
348 | s_kernel = self.S.make_kernel(q_vectors) |
---|
349 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
350 | make_kernel.__doc__ = KernelModel.make_kernel.__doc__ |
---|
351 | |
---|
352 | def release(self): |
---|
353 | # type: (None) -> None |
---|
354 | """ |
---|
355 | Free resources associated with the model. |
---|
356 | """ |
---|
357 | self.P.release() |
---|
358 | self.S.release() |
---|
359 | |
---|
360 | |
---|
361 | class ProductKernel(Kernel): |
---|
362 | """ |
---|
363 | Instantiated kernel for product model. |
---|
364 | """ |
---|
365 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
366 | # type: (ModelInfo, Kernel, Kernel) -> None |
---|
367 | self.info = model_info |
---|
368 | self.p_kernel = p_kernel |
---|
369 | self.s_kernel = s_kernel |
---|
370 | self.dtype = p_kernel.dtype |
---|
371 | self.results = [] # type: List[np.ndarray] |
---|
372 | |
---|
373 | # Find index of volfraction parameter in parameter list |
---|
374 | for k, p in enumerate(model_info.parameters.call_parameters): |
---|
375 | if p.id == VOLFRAC_ID: |
---|
376 | self._volfrac_index = k |
---|
377 | break |
---|
378 | else: |
---|
379 | raise RuntimeError("no %s parameter in %s"%(VOLFRAC_ID, self)) |
---|
380 | |
---|
381 | p_info, s_info = self.info.composition[1] |
---|
382 | p_npars = p_info.parameters.npars |
---|
383 | s_npars = s_info.parameters.npars |
---|
384 | |
---|
385 | have_beta_mode = p_info.have_Fq |
---|
386 | have_er_mode = p_info.radius_effective_modes is not None |
---|
387 | volfrac_in_p = self._volfrac_index < p_npars + 2 # scale & background |
---|
388 | |
---|
389 | # Slices into the details length/offset structure for P@S. |
---|
390 | # Made complicated by the possibly missing volfraction in S. |
---|
391 | self._p_detail_slice = slice(0, p_npars) |
---|
392 | self._s_detail_slice = slice(p_npars, p_npars+s_npars-volfrac_in_p) |
---|
393 | self._volfrac_in_p = volfrac_in_p |
---|
394 | |
---|
395 | # P block from data vector, without scale and background |
---|
396 | first_p = 2 |
---|
397 | last_p = p_npars + 2 |
---|
398 | self._p_value_slice = slice(first_p, last_p) |
---|
399 | |
---|
400 | # radius_effective is the first parameter in S from the data vector. |
---|
401 | self._er_index = last_p |
---|
402 | |
---|
403 | # S block from data vector, without scale, background, volfrac or er. |
---|
404 | first_s = last_p + 2 - volfrac_in_p |
---|
405 | last_s = first_s + s_npars - 2 |
---|
406 | self._s_value_slice = slice(first_s, last_s) |
---|
407 | |
---|
408 | # S distribution block in S data vector starts after all S values |
---|
409 | self._s_dist_slice = slice(2 + s_npars, None) |
---|
410 | |
---|
411 | # structure_factor_mode is the first parameter after P and S. Skip |
---|
412 | # 2 for scale and background, and subtract 1 in case there is no |
---|
413 | # volfraction in S. |
---|
414 | self._beta_mode_index = last_s if have_beta_mode else 0 |
---|
415 | |
---|
416 | # radius_effective_mode is the second parameter after P and S |
---|
417 | # unless structure_factor_mode isn't available, in which case it |
---|
418 | # is first. |
---|
419 | self._er_mode_index = last_s + have_beta_mode if have_er_mode else 0 |
---|
420 | |
---|
421 | # Magnetic parameters are after everything else. If they exist, |
---|
422 | # they will only be for form factor P, not structure factor S. |
---|
423 | first_mag = last_s + have_beta_mode + have_er_mode |
---|
424 | mag_pars = 3*p_info.parameters.nmagnetic |
---|
425 | last_mag = first_mag + (mag_pars + 3 if mag_pars else 0) |
---|
426 | self._magentic_slice = slice(first_mag, last_mag) |
---|
427 | |
---|
428 | def Iq(self, call_details, values, cutoff, magnetic): |
---|
429 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
---|
430 | p_info, s_info = self.info.composition[1] |
---|
431 | |
---|
432 | # Retrieve values from the data vector |
---|
433 | scale, background = values[0], values[1] |
---|
434 | volfrac = values[self._volfrac_index] |
---|
435 | er_mode = (int(values[self._er_mode_index]) |
---|
436 | if self._er_mode_index > 0 else 0) |
---|
437 | beta_mode = (values[self._beta_mode_index] > 0 |
---|
438 | if self._beta_mode_index > 0 else False) |
---|
439 | |
---|
440 | nvalues = self.info.parameters.nvalues |
---|
441 | nweights = call_details.num_weights |
---|
442 | weights = values[nvalues:nvalues + 2*nweights] |
---|
443 | |
---|
444 | # Can't do 2d and beta_mode just yet |
---|
445 | if beta_mode and self.p_kernel.dim == '2d': |
---|
446 | raise NotImplementedError("beta not yet supported for 2D") |
---|
447 | |
---|
448 | # Construct the calling parameters for P. |
---|
449 | p_length = call_details.length[self._p_detail_slice] |
---|
450 | p_offset = call_details.offset[self._p_detail_slice] |
---|
451 | p_details = make_details(p_info, p_length, p_offset, nweights) |
---|
452 | p_values = [ |
---|
453 | [1., 0.], # scale=1, background=0, |
---|
454 | values[self._p_value_slice], |
---|
455 | values[self._magentic_slice], |
---|
456 | weights] |
---|
457 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
---|
458 | p_values.append([0.]*spacer) |
---|
459 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
---|
460 | |
---|
461 | # Call the form factor kernel to compute <F> and <F^2>. |
---|
462 | # If the model doesn't support Fq the returned <F> will be None. |
---|
463 | F, Fsq, radius_effective, shell_volume, volume_ratio \ |
---|
464 | = self.p_kernel.Fq(p_details, p_values, cutoff, magnetic, er_mode) |
---|
465 | |
---|
466 | # TODO: async call to the GPU |
---|
467 | |
---|
468 | # Construct the calling parameters for S. |
---|
469 | s_length = call_details.length[self._s_detail_slice] |
---|
470 | s_offset = call_details.offset[self._s_detail_slice] |
---|
471 | if self._volfrac_in_p: |
---|
472 | # Volfrac is in P and missing from S so insert a slot for it. Say |
---|
473 | # the distribution is length 1 and use the slot for volfraction |
---|
474 | # from the P distribution. |
---|
475 | s_length = np.insert(s_length, 1, 1) |
---|
476 | s_offset = np.insert(s_offset, 1, p_offset[self._volfrac_index - 2]) |
---|
477 | if er_mode > 0: |
---|
478 | # If effective_radius comes from P, make sure it is monodisperse. |
---|
479 | # Weight is set to 1 later, after the value array is created |
---|
480 | s_length[0] = 1 |
---|
481 | s_details = make_details(s_info, s_length, s_offset, nweights) |
---|
482 | s_values = [ |
---|
483 | [1., # scale=1 |
---|
484 | 0., # background=0, |
---|
485 | values[self._er_index], # S.radius_effective; may be replaced by P |
---|
486 | 0.], # volfraction; will be replaced by volfrac * volume_ratio |
---|
487 | # followed by S parameters after effective_radius and volfraction |
---|
488 | values[self._s_value_slice], |
---|
489 | weights, |
---|
490 | ] |
---|
491 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
---|
492 | s_values.append([0.]*spacer) |
---|
493 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
---|
494 | |
---|
495 | # Plug R_eff from the form factor into structure factor parameters |
---|
496 | # and scale volume fraction by form:shell volume ratio. These changes |
---|
497 | # needs to be both in the initial value slot as well as the |
---|
498 | # polydispersity distribution slot in the values array due to |
---|
499 | # implementation details in kernel_iq.c. |
---|
500 | #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g" |
---|
501 | # % (radius_type, radius_effective, volfrac, volume_ratio)) |
---|
502 | s_dist = s_values[self._s_dist_slice] |
---|
503 | if er_mode > 0: |
---|
504 | # set the value to the model R_eff and set the weight to 1 |
---|
505 | s_values[2] = s_dist[s_offset[0]] = radius_effective |
---|
506 | s_dist[s_offset[0]+nweights] = 1.0 |
---|
507 | s_values[3] = s_dist[s_offset[1]] = volfrac*volume_ratio |
---|
508 | s_dist[s_offset[1]+nweights] = 1.0 |
---|
509 | |
---|
510 | # Call the structure factor kernel to compute S. |
---|
511 | S = self.s_kernel.Iq(s_details, s_values, cutoff, False) |
---|
512 | #print("P", Fsq[:10]) |
---|
513 | #print("S", S[:10]) |
---|
514 | #print(radius_effective, volfrac*volume_ratio) |
---|
515 | |
---|
516 | # Combine form factor and structure factor |
---|
517 | #print("beta", beta_mode, F, Fsq, S) |
---|
518 | PS = Fsq + F**2*(S-1) if beta_mode else Fsq*S |
---|
519 | |
---|
520 | # Determine overall scale factor. Hollow shapes are weighted by |
---|
521 | # shell_volume, so that is needed for number density estimation. |
---|
522 | # For solid shapes we can use shell_volume as well since it is |
---|
523 | # equal to form volume. If P already has a volfraction parameter, |
---|
524 | # then assume that it is already on absolute scale, and don't |
---|
525 | # include volfrac in the combined_scale. |
---|
526 | combined_scale = scale*(volfrac if not self._volfrac_in_p else 1.0) |
---|
527 | final_result = combined_scale/shell_volume*PS + background |
---|
528 | |
---|
529 | # Capture intermediate values so user can see them. These are |
---|
530 | # returned as a lazy evaluator since they are only needed in the |
---|
531 | # GUI, and not for each evaluation during a fit. |
---|
532 | # TODO: return the results structure with the final results |
---|
533 | # That way the model calcs are idempotent. Further, we can |
---|
534 | # generalize intermediates to various other model types if we put it |
---|
535 | # kernel calling interface. Could do this as an "optional" |
---|
536 | # return value in the caller, though in that case we could return |
---|
537 | # the results directly rather than through a lazy evaluator. |
---|
538 | self.results = lambda: _intermediates( |
---|
539 | F, Fsq, S, combined_scale, radius_effective, beta_mode) |
---|
540 | |
---|
541 | return final_result |
---|
542 | |
---|
543 | Iq.__doc__ = Kernel.Iq.__doc__ |
---|
544 | __call__ = Iq |
---|
545 | |
---|
546 | def release(self): |
---|
547 | # type: () -> None |
---|
548 | """Free resources associated with the kernel.""" |
---|
549 | self.p_kernel.release() |
---|
550 | self.s_kernel.release() |
---|