| 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 | if not s_info.parameters.kernel_parameters[0].id == ER_ID: |
|---|
| 48 | raise TypeError("S needs %s as first parameter"%ER_ID) |
|---|
| 49 | if not s_info.parameters.kernel_parameters[1].id == VF_ID: |
|---|
| 50 | raise TypeError("S needs %s as second parameter"%VF_ID) |
|---|
| 51 | if not s_info.parameters.magnetism_index == []: |
|---|
| 52 | raise TypeError("S should not have SLD parameters") |
|---|
| 53 | p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters |
|---|
| 54 | s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters |
|---|
| 55 | |
|---|
| 56 | # Create list of parameters for the combined model. Skip the first |
|---|
| 57 | # parameter of S, which we verified above is effective radius. If there |
|---|
| 58 | # are any names in P that overlap with those in S, modify the name in S |
|---|
| 59 | # to distinguish it. |
|---|
| 60 | p_set = set(p.id for p in p_pars.kernel_parameters) |
|---|
| 61 | s_list = [(_tag_parameter(par) if par.id in p_set else par) |
|---|
| 62 | for par in s_pars.kernel_parameters[1:]] |
|---|
| 63 | # Check if still a collision after renaming. This could happen if for |
|---|
| 64 | # example S has volfrac and P has both volfrac and volfrac_S. |
|---|
| 65 | if any(p.id in p_set for p in s_list): |
|---|
| 66 | raise TypeError("name collision: P has P.name and P.name_S while S has S.name") |
|---|
| 67 | |
|---|
| 68 | translate_name = dict((old.id, new.id) for old, new |
|---|
| 69 | in zip(s_pars.kernel_parameters[1:], s_list)) |
|---|
| 70 | demo = {} |
|---|
| 71 | demo.update(p_info.demo.items()) |
|---|
| 72 | demo.update((translate_name[k], v) for k, v in s_info.demo.items() |
|---|
| 73 | if k not in ("background", "scale") and not k.startswith(ER_ID)) |
|---|
| 74 | combined_pars = p_pars.kernel_parameters + s_list |
|---|
| 75 | parameters = ParameterTable(combined_pars) |
|---|
| 76 | parameters.max_pd = p_pars.max_pd + s_pars.max_pd |
|---|
| 77 | |
|---|
| 78 | model_info = ModelInfo() |
|---|
| 79 | model_info.id = '*'.join((p_id, s_id)) |
|---|
| 80 | model_info.name = '*'.join((p_name, s_name)) |
|---|
| 81 | model_info.filename = None |
|---|
| 82 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
|---|
| 83 | model_info.description = model_info.title |
|---|
| 84 | model_info.docs = model_info.title |
|---|
| 85 | model_info.category = "custom" |
|---|
| 86 | model_info.parameters = parameters |
|---|
| 87 | #model_info.single = p_info.single and s_info.single |
|---|
| 88 | model_info.structure_factor = False |
|---|
| 89 | model_info.variant_info = None |
|---|
| 90 | #model_info.tests = [] |
|---|
| 91 | #model_info.source = [] |
|---|
| 92 | # Iq, Iqxy, form_volume, ER, VR and sesans |
|---|
| 93 | # Remember the component info blocks so we can build the model |
|---|
| 94 | model_info.composition = ('product', [p_info, s_info]) |
|---|
| 95 | # TODO: delegate random to p_info, s_info |
|---|
| 96 | #model_info.random = lambda: {} |
|---|
| 97 | model_info.demo = demo |
|---|
| 98 | |
|---|
| 99 | ## Show the parameter table with the demo values |
|---|
| 100 | #from .compare import get_pars, parlist |
|---|
| 101 | #print("==== %s ====="%model_info.name) |
|---|
| 102 | #values = get_pars(model_info, use_demo=True) |
|---|
| 103 | #print(parlist(model_info, values, is2d=True)) |
|---|
| 104 | return model_info |
|---|
| 105 | |
|---|
| 106 | def _tag_parameter(par): |
|---|
| 107 | """ |
|---|
| 108 | Tag the parameter name with _S to indicate that the parameter comes from |
|---|
| 109 | the structure factor parameter set. This is only necessary if the |
|---|
| 110 | form factor model includes a parameter of the same name as a parameter |
|---|
| 111 | in the structure factor. |
|---|
| 112 | """ |
|---|
| 113 | par = copy(par) |
|---|
| 114 | # Protect against a vector parameter in S by appending the vector length |
|---|
| 115 | # to the renamed parameter. Note: haven't tested this since no existing |
|---|
| 116 | # structure factor models contain vector parameters. |
|---|
| 117 | vector_length = par.name[len(par.id):] |
|---|
| 118 | par.id = par.id + "_S" |
|---|
| 119 | par.name = par.id + vector_length |
|---|
| 120 | return par |
|---|
| 121 | |
|---|
| 122 | class ProductModel(KernelModel): |
|---|
| 123 | def __init__(self, model_info, P, S): |
|---|
| 124 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
|---|
| 125 | self.info = model_info |
|---|
| 126 | self.P = P |
|---|
| 127 | self.S = S |
|---|
| 128 | |
|---|
| 129 | def make_kernel(self, q_vectors): |
|---|
| 130 | # type: (List[np.ndarray]) -> Kernel |
|---|
| 131 | # Note: may be sending the q_vectors to the GPU twice even though they |
|---|
| 132 | # are only needed once. It would mess up modularity quite a bit to |
|---|
| 133 | # handle this optimally, especially since there are many cases where |
|---|
| 134 | # separate q vectors are needed (e.g., form in python and structure |
|---|
| 135 | # in opencl; or both in opencl, but one in single precision and the |
|---|
| 136 | # other in double precision). |
|---|
| 137 | p_kernel = self.P.make_kernel(q_vectors) |
|---|
| 138 | s_kernel = self.S.make_kernel(q_vectors) |
|---|
| 139 | return ProductKernel(self.info, p_kernel, s_kernel) |
|---|
| 140 | |
|---|
| 141 | def release(self): |
|---|
| 142 | # type: (None) -> None |
|---|
| 143 | """ |
|---|
| 144 | Free resources associated with the model. |
|---|
| 145 | """ |
|---|
| 146 | self.P.release() |
|---|
| 147 | self.S.release() |
|---|
| 148 | |
|---|
| 149 | |
|---|
| 150 | class ProductKernel(Kernel): |
|---|
| 151 | def __init__(self, model_info, p_kernel, s_kernel): |
|---|
| 152 | # type: (ModelInfo, Kernel, Kernel) -> None |
|---|
| 153 | self.info = model_info |
|---|
| 154 | self.p_kernel = p_kernel |
|---|
| 155 | self.s_kernel = s_kernel |
|---|
| 156 | self.dtype = p_kernel.dtype |
|---|
| 157 | self.results = [] # type: List[np.ndarray] |
|---|
| 158 | |
|---|
| 159 | def __call__(self, call_details, values, cutoff, magnetic): |
|---|
| 160 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
|---|
| 161 | p_info, s_info = self.info.composition[1] |
|---|
| 162 | |
|---|
| 163 | # if there are magnetic parameters, they will only be on the |
|---|
| 164 | # form factor P, not the structure factor S. |
|---|
| 165 | nmagnetic = len(self.info.parameters.magnetism_index) |
|---|
| 166 | if nmagnetic: |
|---|
| 167 | spin_index = self.info.parameters.npars + 2 |
|---|
| 168 | magnetism = values[spin_index: spin_index+3+3*nmagnetic] |
|---|
| 169 | else: |
|---|
| 170 | magnetism = [] |
|---|
| 171 | nvalues = self.info.parameters.nvalues |
|---|
| 172 | nweights = call_details.num_weights |
|---|
| 173 | weights = values[nvalues:nvalues + 2*nweights] |
|---|
| 174 | |
|---|
| 175 | # Construct the calling parameters for P. |
|---|
| 176 | p_npars = p_info.parameters.npars |
|---|
| 177 | p_length = call_details.length[:p_npars] |
|---|
| 178 | p_offset = call_details.offset[:p_npars] |
|---|
| 179 | p_details = make_details(p_info, p_length, p_offset, nweights) |
|---|
| 180 | # Set p scale to the volume fraction in s, which is the first of the |
|---|
| 181 | # 'S' parameters in the parameter list, or 2+np in 0-origin. |
|---|
| 182 | volfrac = values[2+p_npars] |
|---|
| 183 | p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights] |
|---|
| 184 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
|---|
| 185 | p_values.append([0.]*spacer) |
|---|
| 186 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
|---|
| 187 | |
|---|
| 188 | # Call ER and VR for P since these are needed for S. |
|---|
| 189 | p_er, p_vr = calc_er_vr(p_info, p_details, p_values) |
|---|
| 190 | s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) |
|---|
| 191 | #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) |
|---|
| 192 | |
|---|
| 193 | # Construct the calling parameters for S. |
|---|
| 194 | # The effective radius is not in the combined parameter list, so |
|---|
| 195 | # the number of 'S' parameters is one less than expected. The |
|---|
| 196 | # computed effective radius needs to be added into the weights |
|---|
| 197 | # vector, especially since it is a polydisperse parameter in the |
|---|
| 198 | # stand-alone structure factor models. We will added it at the |
|---|
| 199 | # end so the remaining offsets don't need to change. |
|---|
| 200 | s_npars = s_info.parameters.npars-1 |
|---|
| 201 | s_length = call_details.length[p_npars:p_npars+s_npars] |
|---|
| 202 | s_offset = call_details.offset[p_npars:p_npars+s_npars] |
|---|
| 203 | s_length = np.hstack((1, s_length)) |
|---|
| 204 | s_offset = np.hstack((nweights, s_offset)) |
|---|
| 205 | s_details = make_details(s_info, s_length, s_offset, nweights+1) |
|---|
| 206 | v, w = weights[:nweights], weights[nweights:] |
|---|
| 207 | s_values = [ |
|---|
| 208 | # scale=1, background=0, radius_effective=p_er, volfraction=s_vr |
|---|
| 209 | [1., 0., p_er, s_vr], |
|---|
| 210 | # structure factor parameters start after scale, background and |
|---|
| 211 | # all the form factor parameters. Skip the volfraction parameter |
|---|
| 212 | # as well, since it is computed elsewhere, and go to the end of the |
|---|
| 213 | # parameter list. |
|---|
| 214 | values[2+p_npars+1:2+p_npars+s_npars], |
|---|
| 215 | # no magnetism parameters to include for S |
|---|
| 216 | # add er into the (value, weights) pairs |
|---|
| 217 | v, [p_er], w, [1.0] |
|---|
| 218 | ] |
|---|
| 219 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
|---|
| 220 | s_values.append([0.]*spacer) |
|---|
| 221 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
|---|
| 222 | |
|---|
| 223 | # Call the kernels |
|---|
| 224 | p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) |
|---|
| 225 | s_result = self.s_kernel(s_details, s_values, cutoff, False) |
|---|
| 226 | |
|---|
| 227 | #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) |
|---|
| 228 | #call_details.show(values) |
|---|
| 229 | #print("values", values) |
|---|
| 230 | #p_details.show(p_values) |
|---|
| 231 | #print("=>", p_result) |
|---|
| 232 | #s_details.show(s_values) |
|---|
| 233 | #print("=>", s_result) |
|---|
| 234 | |
|---|
| 235 | # remember the parts for plotting later |
|---|
| 236 | self.results = [p_result, s_result] |
|---|
| 237 | |
|---|
| 238 | #import pylab as plt |
|---|
| 239 | #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') |
|---|
| 240 | #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') |
|---|
| 241 | #plt.figure() |
|---|
| 242 | |
|---|
| 243 | return values[0]*(p_result*s_result) + values[1] |
|---|
| 244 | |
|---|
| 245 | def release(self): |
|---|
| 246 | # type: () -> None |
|---|
| 247 | self.p_kernel.release() |
|---|
| 248 | self.s_kernel.release() |
|---|
| 249 | |
|---|
| 250 | |
|---|
| 251 | def calc_er_vr(model_info, call_details, values): |
|---|
| 252 | # type: (ModelInfo, ParameterSet) -> Tuple[float, float] |
|---|
| 253 | |
|---|
| 254 | if model_info.ER is None and model_info.VR is None: |
|---|
| 255 | return 1.0, 1.0 |
|---|
| 256 | |
|---|
| 257 | nvalues = model_info.parameters.nvalues |
|---|
| 258 | value = values[nvalues:nvalues + call_details.num_weights] |
|---|
| 259 | weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights] |
|---|
| 260 | npars = model_info.parameters.npars |
|---|
| 261 | pairs = [(value[offset:offset+length], weight[offset:offset+length]) |
|---|
| 262 | for p, offset, length |
|---|
| 263 | in zip(model_info.parameters.call_parameters[2:2+npars], |
|---|
| 264 | call_details.offset, |
|---|
| 265 | call_details.length) |
|---|
| 266 | if p.type == 'volume'] |
|---|
| 267 | value, weight = dispersion_mesh(model_info, pairs) |
|---|
| 268 | |
|---|
| 269 | if model_info.ER is not None: |
|---|
| 270 | individual_radii = model_info.ER(*value) |
|---|
| 271 | radius_effective = np.sum(weight*individual_radii) / np.sum(weight) |
|---|
| 272 | else: |
|---|
| 273 | radius_effective = 1.0 |
|---|
| 274 | |
|---|
| 275 | if model_info.VR is not None: |
|---|
| 276 | whole, part = model_info.VR(*value) |
|---|
| 277 | volume_ratio = np.sum(weight*part)/np.sum(weight*whole) |
|---|
| 278 | else: |
|---|
| 279 | volume_ratio = 1.0 |
|---|
| 280 | |
|---|
| 281 | return radius_effective, volume_ratio |
|---|
