[17bbadd] | 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 |
---|
[6dc78e4] | 11 | *make_product_info(P, S)*. |
---|
[17bbadd] | 12 | """ |
---|
[6dc78e4] | 13 | from __future__ import print_function, division |
---|
| 14 | |
---|
| 15 | from copy import copy |
---|
[7ae2b7f] | 16 | import numpy as np # type: ignore |
---|
[17bbadd] | 17 | |
---|
[01c8d9e] | 18 | from .modelinfo import ParameterTable, ModelInfo, Parameter |
---|
[9eb3632] | 19 | from .kernel import KernelModel, Kernel |
---|
[6dc78e4] | 20 | from .details import make_details, dispersion_mesh |
---|
[f619de7] | 21 | |
---|
[2d81cfe] | 22 | # pylint: disable=unused-import |
---|
[f619de7] | 23 | try: |
---|
| 24 | from typing import Tuple |
---|
| 25 | except ImportError: |
---|
| 26 | pass |
---|
[2d81cfe] | 27 | else: |
---|
| 28 | from .modelinfo import ParameterSet |
---|
| 29 | # pylint: enable=unused-import |
---|
[17bbadd] | 30 | |
---|
[6d6508e] | 31 | # TODO: make estimates available to constraints |
---|
| 32 | #ESTIMATED_PARAMETERS = [ |
---|
[6dc78e4] | 33 | # ["est_radius_effective", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], |
---|
[6d6508e] | 34 | # ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], |
---|
| 35 | #] |
---|
[17bbadd] | 36 | |
---|
[6dc78e4] | 37 | ER_ID = "radius_effective" |
---|
| 38 | VF_ID = "volfraction" |
---|
[c036ddb] | 39 | BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "", |
---|
| 40 | "Structure factor dispersion calculation mode") |
---|
[6dc78e4] | 41 | |
---|
[3c6d5bc] | 42 | # TODO: core_shell_sphere model has suppressed the volume ratio calculation |
---|
| 43 | # revert it after making VR and ER available at run time as constraints. |
---|
[17bbadd] | 44 | def make_product_info(p_info, s_info): |
---|
[f619de7] | 45 | # type: (ModelInfo, ModelInfo) -> ModelInfo |
---|
[17bbadd] | 46 | """ |
---|
| 47 | Create info block for product model. |
---|
| 48 | """ |
---|
[6dc78e4] | 49 | # Make sure effective radius is the first parameter and |
---|
| 50 | # make sure volume fraction is the second parameter of the |
---|
| 51 | # structure factor calculator. Structure factors should not |
---|
| 52 | # have any magnetic parameters |
---|
[058460c] | 53 | if not len(s_info.parameters.kernel_parameters) >= 2: |
---|
| 54 | raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID)) |
---|
[f88e248] | 55 | if not s_info.parameters.kernel_parameters[0].id == ER_ID: |
---|
[058460c] | 56 | raise TypeError("S needs {} as first parameter".format(ER_ID)) |
---|
[f88e248] | 57 | if not s_info.parameters.kernel_parameters[1].id == VF_ID: |
---|
[058460c] | 58 | raise TypeError("S needs {} as second parameter".format(VF_ID)) |
---|
[f88e248] | 59 | if not s_info.parameters.magnetism_index == []: |
---|
| 60 | raise TypeError("S should not have SLD parameters") |
---|
[f619de7] | 61 | p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters |
---|
| 62 | s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters |
---|
[6d6508e] | 63 | |
---|
[f88e248] | 64 | # Create list of parameters for the combined model. Skip the first |
---|
| 65 | # parameter of S, which we verified above is effective radius. If there |
---|
| 66 | # are any names in P that overlap with those in S, modify the name in S |
---|
| 67 | # to distinguish it. |
---|
| 68 | p_set = set(p.id for p in p_pars.kernel_parameters) |
---|
| 69 | s_list = [(_tag_parameter(par) if par.id in p_set else par) |
---|
| 70 | for par in s_pars.kernel_parameters[1:]] |
---|
| 71 | # Check if still a collision after renaming. This could happen if for |
---|
| 72 | # example S has volfrac and P has both volfrac and volfrac_S. |
---|
| 73 | if any(p.id in p_set for p in s_list): |
---|
| 74 | raise TypeError("name collision: P has P.name and P.name_S while S has S.name") |
---|
| 75 | |
---|
[6dc78e4] | 76 | translate_name = dict((old.id, new.id) for old, new |
---|
| 77 | in zip(s_pars.kernel_parameters[1:], s_list)) |
---|
[c036ddb] | 78 | beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else [] |
---|
| 79 | combined_pars = p_pars.kernel_parameters + s_list + beta |
---|
[f619de7] | 80 | parameters = ParameterTable(combined_pars) |
---|
[6dc78e4] | 81 | parameters.max_pd = p_pars.max_pd + s_pars.max_pd |
---|
[765eb0e] | 82 | def random(): |
---|
| 83 | combined_pars = p_info.random() |
---|
| 84 | s_names = set(par.id for par in s_pars.kernel_parameters[1:]) |
---|
| 85 | combined_pars.update((translate_name[k], v) |
---|
[2d81cfe] | 86 | for k, v in s_info.random().items() |
---|
| 87 | if k in s_names) |
---|
[765eb0e] | 88 | return combined_pars |
---|
[6d6508e] | 89 | |
---|
| 90 | model_info = ModelInfo() |
---|
[6a5ccfb] | 91 | model_info.id = '@'.join((p_id, s_id)) |
---|
| 92 | model_info.name = '@'.join((p_name, s_name)) |
---|
[6d6508e] | 93 | model_info.filename = None |
---|
| 94 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
---|
| 95 | model_info.description = model_info.title |
---|
| 96 | model_info.docs = model_info.title |
---|
| 97 | model_info.category = "custom" |
---|
[f619de7] | 98 | model_info.parameters = parameters |
---|
[765eb0e] | 99 | model_info.random = random |
---|
[6d6508e] | 100 | #model_info.single = p_info.single and s_info.single |
---|
| 101 | model_info.structure_factor = False |
---|
| 102 | model_info.variant_info = None |
---|
| 103 | #model_info.tests = [] |
---|
| 104 | #model_info.source = [] |
---|
[fcd7bbd] | 105 | # Iq, Iqxy, form_volume, ER, VR and sesans |
---|
[6dc78e4] | 106 | # Remember the component info blocks so we can build the model |
---|
[6d6508e] | 107 | model_info.composition = ('product', [p_info, s_info]) |
---|
[1f35235] | 108 | model_info.control = p_info.control |
---|
[439ffcd] | 109 | model_info.hidden = p_info.hidden |
---|
[ee95012] | 110 | if getattr(p_info, 'profile', None) is not None: |
---|
[edb0f85] | 111 | profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) |
---|
[ee95012] | 112 | def profile(**kwargs): |
---|
[edb0f85] | 113 | # extract the profile args |
---|
| 114 | kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars) |
---|
| 115 | return p_info.profile(**kwargs) |
---|
[ee95012] | 116 | else: |
---|
| 117 | profile = None |
---|
| 118 | model_info.profile = profile |
---|
[439ffcd] | 119 | model_info.profile_axes = p_info.profile_axes |
---|
[edb0f85] | 120 | |
---|
[8f04da4] | 121 | # TODO: delegate random to p_info, s_info |
---|
| 122 | #model_info.random = lambda: {} |
---|
[f88e248] | 123 | |
---|
[765eb0e] | 124 | ## Show the parameter table |
---|
[f88e248] | 125 | #from .compare import get_pars, parlist |
---|
| 126 | #print("==== %s ====="%model_info.name) |
---|
[765eb0e] | 127 | #values = get_pars(model_info) |
---|
[f88e248] | 128 | #print(parlist(model_info, values, is2d=True)) |
---|
[17bbadd] | 129 | return model_info |
---|
| 130 | |
---|
[f88e248] | 131 | def _tag_parameter(par): |
---|
| 132 | """ |
---|
| 133 | Tag the parameter name with _S to indicate that the parameter comes from |
---|
| 134 | the structure factor parameter set. This is only necessary if the |
---|
| 135 | form factor model includes a parameter of the same name as a parameter |
---|
| 136 | in the structure factor. |
---|
| 137 | """ |
---|
[6dc78e4] | 138 | par = copy(par) |
---|
[f88e248] | 139 | # Protect against a vector parameter in S by appending the vector length |
---|
| 140 | # to the renamed parameter. Note: haven't tested this since no existing |
---|
| 141 | # structure factor models contain vector parameters. |
---|
| 142 | vector_length = par.name[len(par.id):] |
---|
[6dc78e4] | 143 | par.id = par.id + "_S" |
---|
[f88e248] | 144 | par.name = par.id + vector_length |
---|
[6dc78e4] | 145 | return par |
---|
| 146 | |
---|
[f619de7] | 147 | class ProductModel(KernelModel): |
---|
[72a081d] | 148 | def __init__(self, model_info, P, S): |
---|
[f619de7] | 149 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
---|
[146793b] | 150 | #: Combined info plock for the product model |
---|
[72a081d] | 151 | self.info = model_info |
---|
[146793b] | 152 | #: Form factor modelling individual particles. |
---|
[17bbadd] | 153 | self.P = P |
---|
[146793b] | 154 | #: Structure factor modelling interaction between particles. |
---|
[17bbadd] | 155 | self.S = S |
---|
[c036ddb] | 156 | |
---|
[146793b] | 157 | #: Model precision. This is not really relevant, since it is the |
---|
| 158 | #: individual P and S models that control the effective dtype, |
---|
| 159 | #: converting the q-vectors to the correct type when the kernels |
---|
| 160 | #: for each are created. Ideally this should be set to the more |
---|
| 161 | #: precise type to avoid loss of precision, but precision in q is |
---|
| 162 | #: not critical (single is good enough for our purposes), so it just |
---|
| 163 | #: uses the precision of the form factor. |
---|
| 164 | self.dtype = P.dtype # type: np.dtype |
---|
[17bbadd] | 165 | |
---|
[6dc78e4] | 166 | def make_kernel(self, q_vectors): |
---|
[f619de7] | 167 | # type: (List[np.ndarray]) -> Kernel |
---|
[17bbadd] | 168 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
| 169 | # are only needed once. It would mess up modularity quite a bit to |
---|
| 170 | # handle this optimally, especially since there are many cases where |
---|
| 171 | # separate q vectors are needed (e.g., form in python and structure |
---|
| 172 | # in opencl; or both in opencl, but one in single precision and the |
---|
| 173 | # other in double precision). |
---|
[c036ddb] | 174 | |
---|
[f619de7] | 175 | p_kernel = self.P.make_kernel(q_vectors) |
---|
| 176 | s_kernel = self.S.make_kernel(q_vectors) |
---|
[17bbadd] | 177 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
| 178 | |
---|
| 179 | def release(self): |
---|
[f619de7] | 180 | # type: (None) -> None |
---|
[17bbadd] | 181 | """ |
---|
| 182 | Free resources associated with the model. |
---|
| 183 | """ |
---|
| 184 | self.P.release() |
---|
| 185 | self.S.release() |
---|
| 186 | |
---|
| 187 | |
---|
[f619de7] | 188 | class ProductKernel(Kernel): |
---|
[17bbadd] | 189 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
[f619de7] | 190 | # type: (ModelInfo, Kernel, Kernel) -> None |
---|
[17bbadd] | 191 | self.info = model_info |
---|
| 192 | self.p_kernel = p_kernel |
---|
| 193 | self.s_kernel = s_kernel |
---|
[6dc78e4] | 194 | self.dtype = p_kernel.dtype |
---|
| 195 | self.results = [] # type: List[np.ndarray] |
---|
| 196 | |
---|
| 197 | def __call__(self, call_details, values, cutoff, magnetic): |
---|
| 198 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
---|
| 199 | p_info, s_info = self.info.composition[1] |
---|
[c036ddb] | 200 | |
---|
[6dc78e4] | 201 | # if there are magnetic parameters, they will only be on the |
---|
| 202 | # form factor P, not the structure factor S. |
---|
| 203 | nmagnetic = len(self.info.parameters.magnetism_index) |
---|
| 204 | if nmagnetic: |
---|
| 205 | spin_index = self.info.parameters.npars + 2 |
---|
| 206 | magnetism = values[spin_index: spin_index+3+3*nmagnetic] |
---|
| 207 | else: |
---|
| 208 | magnetism = [] |
---|
| 209 | nvalues = self.info.parameters.nvalues |
---|
| 210 | nweights = call_details.num_weights |
---|
| 211 | weights = values[nvalues:nvalues + 2*nweights] |
---|
[c036ddb] | 212 | |
---|
[6dc78e4] | 213 | # Construct the calling parameters for P. |
---|
| 214 | p_npars = p_info.parameters.npars |
---|
| 215 | p_length = call_details.length[:p_npars] |
---|
| 216 | p_offset = call_details.offset[:p_npars] |
---|
| 217 | p_details = make_details(p_info, p_length, p_offset, nweights) |
---|
[c036ddb] | 218 | |
---|
[6dc78e4] | 219 | # Set p scale to the volume fraction in s, which is the first of the |
---|
| 220 | # 'S' parameters in the parameter list, or 2+np in 0-origin. |
---|
| 221 | volfrac = values[2+p_npars] |
---|
| 222 | p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights] |
---|
| 223 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
---|
| 224 | p_values.append([0.]*spacer) |
---|
| 225 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
---|
[c036ddb] | 226 | |
---|
[6dc78e4] | 227 | # Call ER and VR for P since these are needed for S. |
---|
| 228 | p_er, p_vr = calc_er_vr(p_info, p_details, p_values) |
---|
| 229 | s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) |
---|
| 230 | #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) |
---|
[c036ddb] | 231 | |
---|
[6dc78e4] | 232 | # Construct the calling parameters for S. |
---|
| 233 | # The effective radius is not in the combined parameter list, so |
---|
| 234 | # the number of 'S' parameters is one less than expected. The |
---|
| 235 | # computed effective radius needs to be added into the weights |
---|
| 236 | # vector, especially since it is a polydisperse parameter in the |
---|
| 237 | # stand-alone structure factor models. We will added it at the |
---|
| 238 | # end so the remaining offsets don't need to change. |
---|
| 239 | s_npars = s_info.parameters.npars-1 |
---|
| 240 | s_length = call_details.length[p_npars:p_npars+s_npars] |
---|
| 241 | s_offset = call_details.offset[p_npars:p_npars+s_npars] |
---|
| 242 | s_length = np.hstack((1, s_length)) |
---|
| 243 | s_offset = np.hstack((nweights, s_offset)) |
---|
| 244 | s_details = make_details(s_info, s_length, s_offset, nweights+1) |
---|
| 245 | v, w = weights[:nweights], weights[nweights:] |
---|
[9951a86] | 246 | s_values = [ |
---|
| 247 | # scale=1, background=0, radius_effective=p_er, volfraction=s_vr |
---|
| 248 | [1., 0., p_er, s_vr], |
---|
| 249 | # structure factor parameters start after scale, background and |
---|
| 250 | # all the form factor parameters. Skip the volfraction parameter |
---|
| 251 | # as well, since it is computed elsewhere, and go to the end of the |
---|
| 252 | # parameter list. |
---|
| 253 | values[2+p_npars+1:2+p_npars+s_npars], |
---|
| 254 | # no magnetism parameters to include for S |
---|
| 255 | # add er into the (value, weights) pairs |
---|
| 256 | v, [p_er], w, [1.0] |
---|
| 257 | ] |
---|
[6dc78e4] | 258 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
---|
| 259 | s_values.append([0.]*spacer) |
---|
| 260 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
---|
[c036ddb] | 261 | |
---|
[01c8d9e] | 262 | # beta mode is the first parameter after the structure factor pars |
---|
| 263 | beta_index = 2+p_npars+s_npars |
---|
| 264 | beta_mode = values[beta_index] |
---|
[c036ddb] | 265 | |
---|
[6dc78e4] | 266 | # Call the kernels |
---|
[c036ddb] | 267 | s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) |
---|
| 268 | scale, background = values[0], values[1] |
---|
| 269 | if beta_mode: |
---|
[01c8d9e] | 270 | F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) |
---|
[c036ddb] | 271 | combined_scale = scale*volfrac/volume_avg |
---|
| 272 | # Define lazy results based on intermediate values. |
---|
| 273 | # The return value for the calculation should be an ordered |
---|
| 274 | # dictionary containing any result the user might want to see |
---|
| 275 | # at the end of the calculation, including scalars, strings, |
---|
| 276 | # and plottable data. Don't want to build this structure during |
---|
| 277 | # fits, only when displaying the final result (or a one-off |
---|
| 278 | # computation which asks for it). |
---|
| 279 | # Do not use the current hack of storing the intermediate values |
---|
| 280 | # in self.results since that leads to awkward threading issues. |
---|
| 281 | # Instead return the function along with the bundle of inputs. |
---|
| 282 | # P and Q may themselves have intermediate results they want to |
---|
| 283 | # include, such as A and B if P = A + B. Might use this mechanism |
---|
| 284 | # to return the computed effective radius as well. |
---|
| 285 | #def lazy_results(Q, S, F1, F2, scale): |
---|
| 286 | # """ |
---|
| 287 | # beta = F1**2 / F2 # what about divide by zero errors? |
---|
| 288 | # return { |
---|
| 289 | # 'P' : Data1D(Q, scale*F2), |
---|
| 290 | # 'beta': Data1D(Q, beta), |
---|
| 291 | # 'S' : Data1D(Q, S), |
---|
| 292 | # 'Seff': Data1D(Q, 1 + beta*(S-1)), |
---|
| 293 | # 'I' : Data1D(Q, scale*(F2 + (F1**2)*(S-1)) + background), |
---|
| 294 | # } |
---|
| 295 | #lazy_pars = s_result, F1, F2, combined_scale |
---|
| 296 | self.results = [F2, s_result] |
---|
| 297 | final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background |
---|
[01c8d9e] | 298 | else: |
---|
| 299 | p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) |
---|
[c036ddb] | 300 | # remember the parts for plotting later |
---|
| 301 | self.results = [p_result, s_result] |
---|
| 302 | final_result = scale*(p_result*s_result) + background |
---|
| 303 | |
---|
[6dc78e4] | 304 | #call_details.show(values) |
---|
| 305 | #print("values", values) |
---|
| 306 | #p_details.show(p_values) |
---|
| 307 | #print("=>", p_result) |
---|
| 308 | #s_details.show(s_values) |
---|
| 309 | #print("=>", s_result) |
---|
| 310 | #import pylab as plt |
---|
| 311 | #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') |
---|
| 312 | #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') |
---|
| 313 | #plt.figure() |
---|
[01c8d9e] | 314 | return final_result |
---|
[17bbadd] | 315 | |
---|
| 316 | def release(self): |
---|
[f619de7] | 317 | # type: () -> None |
---|
[17bbadd] | 318 | self.p_kernel.release() |
---|
[f619de7] | 319 | self.s_kernel.release() |
---|
[17bbadd] | 320 | |
---|
[6dc78e4] | 321 | |
---|
| 322 | def calc_er_vr(model_info, call_details, values): |
---|
| 323 | # type: (ModelInfo, ParameterSet) -> Tuple[float, float] |
---|
| 324 | |
---|
[f619de7] | 325 | if model_info.ER is None and model_info.VR is None: |
---|
| 326 | return 1.0, 1.0 |
---|
| 327 | |
---|
[6dc78e4] | 328 | nvalues = model_info.parameters.nvalues |
---|
| 329 | value = values[nvalues:nvalues + call_details.num_weights] |
---|
| 330 | weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights] |
---|
| 331 | npars = model_info.parameters.npars |
---|
[ce99754] | 332 | # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the |
---|
| 333 | # dispersion mesh code doesn't actually care about the nominal parameter |
---|
| 334 | # value, p, so set it to None. |
---|
| 335 | pairs = [(None, value[offset:offset+length], weight[offset:offset+length]) |
---|
[6dc78e4] | 336 | for p, offset, length |
---|
| 337 | in zip(model_info.parameters.call_parameters[2:2+npars], |
---|
| 338 | call_details.offset, |
---|
| 339 | call_details.length) |
---|
| 340 | if p.type == 'volume'] |
---|
| 341 | value, weight = dispersion_mesh(model_info, pairs) |
---|
[6d6508e] | 342 | |
---|
[f619de7] | 343 | if model_info.ER is not None: |
---|
| 344 | individual_radii = model_info.ER(*value) |
---|
[6dc78e4] | 345 | radius_effective = np.sum(weight*individual_radii) / np.sum(weight) |
---|
[f619de7] | 346 | else: |
---|
[6dc78e4] | 347 | radius_effective = 1.0 |
---|
[f619de7] | 348 | |
---|
| 349 | if model_info.VR is not None: |
---|
| 350 | whole, part = model_info.VR(*value) |
---|
| 351 | volume_ratio = np.sum(weight*part)/np.sum(weight*whole) |
---|
| 352 | else: |
---|
| 353 | volume_ratio = 1.0 |
---|
[6d6508e] | 354 | |
---|
[6dc78e4] | 355 | return radius_effective, volume_ratio |
---|