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 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, dispersion_mesh |
---|
21 | |
---|
22 | try: |
---|
23 | from typing import Tuple |
---|
24 | except ImportError: |
---|
25 | pass |
---|
26 | |
---|
27 | # TODO: make estimates available to constraints |
---|
28 | #ESTIMATED_PARAMETERS = [ |
---|
29 | # ["est_radius_effective", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], |
---|
30 | # ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], |
---|
31 | #] |
---|
32 | |
---|
33 | ER_ID = "radius_effective" |
---|
34 | VF_ID = "volfraction" |
---|
35 | |
---|
36 | # TODO: core_shell_sphere model has suppressed the volume ratio calculation |
---|
37 | # revert it after making VR and ER available at run time as constraints. |
---|
38 | def make_product_info(p_info, s_info): |
---|
39 | # type: (ModelInfo, ModelInfo) -> ModelInfo |
---|
40 | """ |
---|
41 | Create info block for product model. |
---|
42 | """ |
---|
43 | # Make sure effective radius is the first parameter and |
---|
44 | # make sure volume fraction is the second parameter of the |
---|
45 | # structure factor calculator. Structure factors should not |
---|
46 | # have any magnetic parameters |
---|
47 | assert(s_info.parameters.kernel_parameters[0].id == ER_ID) |
---|
48 | assert(s_info.parameters.kernel_parameters[1].id == VF_ID) |
---|
49 | assert(s_info.parameters.magnetism_index == []) |
---|
50 | p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters |
---|
51 | s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters |
---|
52 | p_set = set(p.id for p in p_pars.call_parameters) |
---|
53 | s_set = set(p.id for p in s_pars.call_parameters) |
---|
54 | |
---|
55 | if p_set & s_set: |
---|
56 | # there is some overlap between the parameter names; tag the |
---|
57 | # overlapping S parameters with name_S. |
---|
58 | # Skip the first parameter of s, which is effective radius |
---|
59 | s_list = [(suffix_parameter(par) if par.id in p_set else par) |
---|
60 | for par in s_pars.kernel_parameters[1:]] |
---|
61 | else: |
---|
62 | # Skip the first parameter of s, which is effective radius |
---|
63 | s_list = s_pars.kernel_parameters[1:] |
---|
64 | translate_name = dict((old.id, new.id) for old, new |
---|
65 | in zip(s_pars.kernel_parameters[1:], s_list)) |
---|
66 | demo = {} |
---|
67 | demo.update((k, v) for k, v in p_info.demo.items() |
---|
68 | if k not in ("background", "scale")) |
---|
69 | demo.update((translate_name[k], v) for k, v in s_info.demo.items() |
---|
70 | if k not in ("background", "scale") and not k.startswith(ER_ID)) |
---|
71 | combined_pars = p_pars.kernel_parameters + s_list |
---|
72 | parameters = ParameterTable(combined_pars) |
---|
73 | parameters.max_pd = p_pars.max_pd + s_pars.max_pd |
---|
74 | |
---|
75 | model_info = ModelInfo() |
---|
76 | model_info.id = '*'.join((p_id, s_id)) |
---|
77 | model_info.name = '*'.join((p_name, s_name)) |
---|
78 | model_info.filename = None |
---|
79 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
---|
80 | model_info.description = model_info.title |
---|
81 | model_info.docs = model_info.title |
---|
82 | model_info.category = "custom" |
---|
83 | model_info.parameters = parameters |
---|
84 | #model_info.single = p_info.single and s_info.single |
---|
85 | model_info.structure_factor = False |
---|
86 | model_info.variant_info = None |
---|
87 | #model_info.tests = [] |
---|
88 | #model_info.source = [] |
---|
89 | # Iq, Iqxy, form_volume, ER, VR and sesans |
---|
90 | # Remember the component info blocks so we can build the model |
---|
91 | model_info.composition = ('product', [p_info, s_info]) |
---|
92 | model_info.demo = {} |
---|
93 | return model_info |
---|
94 | |
---|
95 | def suffix_parameter(par, suffix): |
---|
96 | par = copy(par) |
---|
97 | par.name = par.name + " S" |
---|
98 | par.id = par.id + "_S" |
---|
99 | return par |
---|
100 | |
---|
101 | class ProductModel(KernelModel): |
---|
102 | def __init__(self, model_info, P, S): |
---|
103 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
---|
104 | self.info = model_info |
---|
105 | self.P = P |
---|
106 | self.S = S |
---|
107 | |
---|
108 | def make_kernel(self, q_vectors): |
---|
109 | # type: (List[np.ndarray]) -> Kernel |
---|
110 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
111 | # are only needed once. It would mess up modularity quite a bit to |
---|
112 | # handle this optimally, especially since there are many cases where |
---|
113 | # separate q vectors are needed (e.g., form in python and structure |
---|
114 | # in opencl; or both in opencl, but one in single precision and the |
---|
115 | # other in double precision). |
---|
116 | p_kernel = self.P.make_kernel(q_vectors) |
---|
117 | s_kernel = self.S.make_kernel(q_vectors) |
---|
118 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
119 | |
---|
120 | def release(self): |
---|
121 | # type: (None) -> None |
---|
122 | """ |
---|
123 | Free resources associated with the model. |
---|
124 | """ |
---|
125 | self.P.release() |
---|
126 | self.S.release() |
---|
127 | |
---|
128 | |
---|
129 | class ProductKernel(Kernel): |
---|
130 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
131 | # type: (ModelInfo, Kernel, Kernel) -> None |
---|
132 | self.info = model_info |
---|
133 | self.p_kernel = p_kernel |
---|
134 | self.s_kernel = s_kernel |
---|
135 | self.dtype = p_kernel.dtype |
---|
136 | self.results = [] # type: List[np.ndarray] |
---|
137 | |
---|
138 | def __call__(self, call_details, values, cutoff, magnetic): |
---|
139 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
---|
140 | p_info, s_info = self.info.composition[1] |
---|
141 | |
---|
142 | # if there are magnetic parameters, they will only be on the |
---|
143 | # form factor P, not the structure factor S. |
---|
144 | nmagnetic = len(self.info.parameters.magnetism_index) |
---|
145 | if nmagnetic: |
---|
146 | spin_index = self.info.parameters.npars + 2 |
---|
147 | magnetism = values[spin_index: spin_index+3+3*nmagnetic] |
---|
148 | else: |
---|
149 | magnetism = [] |
---|
150 | nvalues = self.info.parameters.nvalues |
---|
151 | nweights = call_details.num_weights |
---|
152 | weights = values[nvalues:nvalues + 2*nweights] |
---|
153 | |
---|
154 | # Construct the calling parameters for P. |
---|
155 | p_npars = p_info.parameters.npars |
---|
156 | p_length = call_details.length[:p_npars] |
---|
157 | p_offset = call_details.offset[:p_npars] |
---|
158 | p_details = make_details(p_info, p_length, p_offset, nweights) |
---|
159 | # Set p scale to the volume fraction in s, which is the first of the |
---|
160 | # 'S' parameters in the parameter list, or 2+np in 0-origin. |
---|
161 | volfrac = values[2+p_npars] |
---|
162 | p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights] |
---|
163 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
---|
164 | p_values.append([0.]*spacer) |
---|
165 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
---|
166 | |
---|
167 | # Call ER and VR for P since these are needed for S. |
---|
168 | p_er, p_vr = calc_er_vr(p_info, p_details, p_values) |
---|
169 | s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) |
---|
170 | #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) |
---|
171 | |
---|
172 | # Construct the calling parameters for S. |
---|
173 | # The effective radius is not in the combined parameter list, so |
---|
174 | # the number of 'S' parameters is one less than expected. The |
---|
175 | # computed effective radius needs to be added into the weights |
---|
176 | # vector, especially since it is a polydisperse parameter in the |
---|
177 | # stand-alone structure factor models. We will added it at the |
---|
178 | # end so the remaining offsets don't need to change. |
---|
179 | s_npars = s_info.parameters.npars-1 |
---|
180 | s_length = call_details.length[p_npars:p_npars+s_npars] |
---|
181 | s_offset = call_details.offset[p_npars:p_npars+s_npars] |
---|
182 | s_length = np.hstack((1, s_length)) |
---|
183 | s_offset = np.hstack((nweights, s_offset)) |
---|
184 | s_details = make_details(s_info, s_length, s_offset, nweights+1) |
---|
185 | v, w = weights[:nweights], weights[nweights:] |
---|
186 | s_values = [[1., 0., p_er, s_vr], |
---|
187 | # er and vf already included, so start one past the |
---|
188 | # volfrac parameter, and go two less than the number |
---|
189 | # of S parameters. |
---|
190 | values[2+p_npars+2:2+p_npars+s_npars-1], |
---|
191 | # no magnetism parameters to include for S |
---|
192 | # add er into the (value, weights) pairs |
---|
193 | v, [p_er], w, [1.0]] |
---|
194 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
---|
195 | s_values.append([0.]*spacer) |
---|
196 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
---|
197 | |
---|
198 | # Call the kernels |
---|
199 | p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) |
---|
200 | s_result = self.s_kernel(s_details, s_values, cutoff, False) |
---|
201 | |
---|
202 | #call_details.show(values) |
---|
203 | #print("values", values) |
---|
204 | #p_details.show(p_values) |
---|
205 | #print("=>", p_result) |
---|
206 | #print("p val", s_values) |
---|
207 | #s_details.show(s_values) |
---|
208 | #print("=>", s_result) |
---|
209 | |
---|
210 | # remember the parts for plotting later |
---|
211 | self.results = [p_result, s_result] |
---|
212 | |
---|
213 | #import pylab as plt |
---|
214 | #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') |
---|
215 | #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') |
---|
216 | #plt.figure() |
---|
217 | |
---|
218 | return values[0]*(p_result*s_result) + values[1] |
---|
219 | |
---|
220 | def release(self): |
---|
221 | # type: () -> None |
---|
222 | self.p_kernel.release() |
---|
223 | self.s_kernel.release() |
---|
224 | |
---|
225 | |
---|
226 | def calc_er_vr(model_info, call_details, values): |
---|
227 | # type: (ModelInfo, ParameterSet) -> Tuple[float, float] |
---|
228 | |
---|
229 | if model_info.ER is None and model_info.VR is None: |
---|
230 | return 1.0, 1.0 |
---|
231 | |
---|
232 | nvalues = model_info.parameters.nvalues |
---|
233 | value = values[nvalues:nvalues + call_details.num_weights] |
---|
234 | weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights] |
---|
235 | npars = model_info.parameters.npars |
---|
236 | pairs = [(value[offset:offset+length], weight[offset:offset+length]) |
---|
237 | for p, offset, length |
---|
238 | in zip(model_info.parameters.call_parameters[2:2+npars], |
---|
239 | call_details.offset, |
---|
240 | call_details.length) |
---|
241 | if p.type == 'volume'] |
---|
242 | value, weight = dispersion_mesh(model_info, pairs) |
---|
243 | |
---|
244 | if model_info.ER is not None: |
---|
245 | individual_radii = model_info.ER(*value) |
---|
246 | radius_effective = np.sum(weight*individual_radii) / np.sum(weight) |
---|
247 | else: |
---|
248 | radius_effective = 1.0 |
---|
249 | |
---|
250 | if model_info.VR is not None: |
---|
251 | whole, part = model_info.VR(*value) |
---|
252 | volume_ratio = np.sum(weight*part)/np.sum(weight*whole) |
---|
253 | else: |
---|
254 | volume_ratio = 1.0 |
---|
255 | |
---|
256 | return radius_effective, volume_ratio |
---|