[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)*. |
---|
[b2f0e5d] | 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 | |
---|
[17bbadd] | 114 | """ |
---|
[6dc78e4] | 115 | from __future__ import print_function, division |
---|
| 116 | |
---|
[d32de68] | 117 | from collections import OrderedDict |
---|
| 118 | |
---|
[6dc78e4] | 119 | from copy import copy |
---|
[7ae2b7f] | 120 | import numpy as np # type: ignore |
---|
[17bbadd] | 121 | |
---|
[b297ba9] | 122 | from .modelinfo import ParameterTable, ModelInfo, parse_parameter |
---|
[9eb3632] | 123 | from .kernel import KernelModel, Kernel |
---|
[b297ba9] | 124 | from .details import make_details |
---|
[f619de7] | 125 | |
---|
[2d81cfe] | 126 | # pylint: disable=unused-import |
---|
[f619de7] | 127 | try: |
---|
[aa44a6a] | 128 | from typing import Tuple, Callable, Union |
---|
[f619de7] | 129 | except ImportError: |
---|
| 130 | pass |
---|
[2d81cfe] | 131 | else: |
---|
[b297ba9] | 132 | from .modelinfo import ParameterSet, Parameter |
---|
[2d81cfe] | 133 | # pylint: enable=unused-import |
---|
[17bbadd] | 134 | |
---|
[a34b811] | 135 | # TODO: make shape averages available to constraints |
---|
[6d6508e] | 136 | #ESTIMATED_PARAMETERS = [ |
---|
[a34b811] | 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"], |
---|
[6d6508e] | 140 | #] |
---|
[0b8a1fc] | 141 | STRUCTURE_MODE_ID = "structure_factor_mode" |
---|
| 142 | RADIUS_MODE_ID = "radius_effective_mode" |
---|
[6e7ba14] | 143 | RADIUS_ID = "radius_effective" |
---|
| 144 | VOLFRAC_ID = "volfraction" |
---|
| 145 | def make_extra_pars(p_info): |
---|
[b297ba9] | 146 | # type: (ModelInfo) -> List[Parameter] |
---|
| 147 | """ |
---|
[a34b811] | 148 | Create parameters for structure factor and effective radius modes. |
---|
[b297ba9] | 149 | """ |
---|
[6e7ba14] | 150 | pars = [] |
---|
| 151 | if p_info.have_Fq: |
---|
[c0131d44] | 152 | par = parse_parameter( |
---|
[b297ba9] | 153 | STRUCTURE_MODE_ID, |
---|
| 154 | "", |
---|
| 155 | 0, |
---|
| 156 | [["P*S", "P*(1+beta*(S-1))"]], |
---|
| 157 | "", |
---|
| 158 | "Structure factor calculation") |
---|
[6e7ba14] | 159 | pars.append(par) |
---|
[a34b811] | 160 | if p_info.radius_effective_modes is not None: |
---|
[2773c66] | 161 | par = parse_parameter( |
---|
[b297ba9] | 162 | RADIUS_MODE_ID, |
---|
| 163 | "", |
---|
| 164 | 1, |
---|
[a34b811] | 165 | [["unconstrained"] + p_info.radius_effective_modes], |
---|
[b297ba9] | 166 | "", |
---|
| 167 | "Effective radius calculation") |
---|
[6e7ba14] | 168 | pars.append(par) |
---|
| 169 | return pars |
---|
[6dc78e4] | 170 | |
---|
[17bbadd] | 171 | def make_product_info(p_info, s_info): |
---|
[f619de7] | 172 | # type: (ModelInfo, ModelInfo) -> ModelInfo |
---|
[17bbadd] | 173 | """ |
---|
| 174 | Create info block for product model. |
---|
| 175 | """ |
---|
[6dc78e4] | 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 |
---|
[058460c] | 180 | if not len(s_info.parameters.kernel_parameters) >= 2: |
---|
[6e7ba14] | 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)) |
---|
[f88e248] | 186 | if not s_info.parameters.magnetism_index == []: |
---|
| 187 | raise TypeError("S should not have SLD parameters") |
---|
[b2f0e5d] | 188 | if RADIUS_ID in p_info.parameters: |
---|
| 189 | raise TypeError("P should not have {}".format(RADIUS_ID)) |
---|
[f619de7] | 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 |
---|
[b2f0e5d] | 192 | p_has_volfrac = VOLFRAC_ID in p_info.parameters |
---|
[6d6508e] | 193 | |
---|
[b2f0e5d] | 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. |
---|
[f88e248] | 199 | p_set = set(p.id for p in p_pars.kernel_parameters) |
---|
[b2f0e5d] | 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] |
---|
[f88e248] | 205 | |
---|
[b2f0e5d] | 206 | # Build combined parameter table |
---|
[6e7ba14] | 207 | combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) |
---|
[f619de7] | 208 | parameters = ParameterTable(combined_pars) |
---|
[b2f0e5d] | 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} |
---|
[765eb0e] | 220 | def random(): |
---|
[b297ba9] | 221 | """Random set of model parameters for product model""" |
---|
[765eb0e] | 222 | combined_pars = p_info.random() |
---|
[b2f0e5d] | 223 | combined_pars.update((s_translate[k], v) |
---|
[2d81cfe] | 224 | for k, v in s_info.random().items() |
---|
[b2f0e5d] | 225 | if k in s_translate) |
---|
[765eb0e] | 226 | return combined_pars |
---|
[6d6508e] | 227 | |
---|
| 228 | model_info = ModelInfo() |
---|
[6a5ccfb] | 229 | model_info.id = '@'.join((p_id, s_id)) |
---|
| 230 | model_info.name = '@'.join((p_name, s_name)) |
---|
[6d6508e] | 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" |
---|
[f619de7] | 236 | model_info.parameters = parameters |
---|
[765eb0e] | 237 | model_info.random = random |
---|
[6d6508e] | 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 = [] |
---|
[6dc78e4] | 243 | # Remember the component info blocks so we can build the model |
---|
[6d6508e] | 244 | model_info.composition = ('product', [p_info, s_info]) |
---|
[439ffcd] | 245 | model_info.hidden = p_info.hidden |
---|
[ee95012] | 246 | if getattr(p_info, 'profile', None) is not None: |
---|
[edb0f85] | 247 | profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) |
---|
[ee95012] | 248 | def profile(**kwargs): |
---|
[b297ba9] | 249 | """Return SLD profile of the form factor as a function of radius.""" |
---|
[edb0f85] | 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) |
---|
[ee95012] | 253 | else: |
---|
| 254 | profile = None |
---|
| 255 | model_info.profile = profile |
---|
[439ffcd] | 256 | model_info.profile_axes = p_info.profile_axes |
---|
[edb0f85] | 257 | |
---|
[8f04da4] | 258 | # TODO: delegate random to p_info, s_info |
---|
| 259 | #model_info.random = lambda: {} |
---|
[f88e248] | 260 | |
---|
[765eb0e] | 261 | ## Show the parameter table |
---|
[f88e248] | 262 | #from .compare import get_pars, parlist |
---|
| 263 | #print("==== %s ====="%model_info.name) |
---|
[765eb0e] | 264 | #values = get_pars(model_info) |
---|
[f88e248] | 265 | #print(parlist(model_info, values, is2d=True)) |
---|
[17bbadd] | 266 | return model_info |
---|
| 267 | |
---|
[f88e248] | 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 | """ |
---|
[6dc78e4] | 275 | par = copy(par) |
---|
[f88e248] | 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):] |
---|
[6dc78e4] | 280 | par.id = par.id + "_S" |
---|
[f88e248] | 281 | par.name = par.id + vector_length |
---|
[6dc78e4] | 282 | return par |
---|
| 283 | |
---|
[e44432d] | 284 | def _intermediates( |
---|
[b2f0e5d] | 285 | F, # type: np.ndarray |
---|
| 286 | Fsq, # type: np.ndarray |
---|
[e44432d] | 287 | S, # type: np.ndarray |
---|
| 288 | scale, # type: float |
---|
[a34b811] | 289 | radius_effective, # type: float |
---|
[e44432d] | 290 | beta_mode, # type: bool |
---|
[b297ba9] | 291 | ): |
---|
[2773c66] | 292 | # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] |
---|
[d32de68] | 293 | """ |
---|
[e44432d] | 294 | Returns intermediate results for beta approximation-enabled product. |
---|
| 295 | The result may be an array or a float. |
---|
[d32de68] | 296 | """ |
---|
[065d77d] | 297 | # CRUFT: remove effective_radius once SasView 5.0 is released. |
---|
[e44432d] | 298 | if beta_mode: |
---|
| 299 | # TODO: 1. include calculated Q vector |
---|
| 300 | # TODO: 2. consider implications if there are intermediate results in P(Q) |
---|
| 301 | parts = OrderedDict(( |
---|
[b2f0e5d] | 302 | ("P(Q)", scale*Fsq), |
---|
[e44432d] | 303 | ("S(Q)", S), |
---|
[b2f0e5d] | 304 | ("beta(Q)", F**2 / Fsq), |
---|
| 305 | ("S_eff(Q)", 1 + (F**2 / Fsq)*(S-1)), |
---|
[a34b811] | 306 | ("effective_radius", radius_effective), |
---|
[065d77d] | 307 | ("radius_effective", radius_effective), |
---|
[b2f0e5d] | 308 | # ("I(Q)", scale*(Fsq + (F**2)*(S-1)) + bg), |
---|
[e44432d] | 309 | )) |
---|
| 310 | else: |
---|
| 311 | parts = OrderedDict(( |
---|
[b2f0e5d] | 312 | ("P(Q)", scale*Fsq), |
---|
[e44432d] | 313 | ("S(Q)", S), |
---|
[a34b811] | 314 | ("effective_radius", radius_effective), |
---|
[065d77d] | 315 | ("radius_effective", radius_effective), |
---|
[e44432d] | 316 | )) |
---|
| 317 | return parts |
---|
[d32de68] | 318 | |
---|
[f619de7] | 319 | class ProductModel(KernelModel): |
---|
[b297ba9] | 320 | """ |
---|
| 321 | Model definition for product model. |
---|
| 322 | """ |
---|
[72a081d] | 323 | def __init__(self, model_info, P, S): |
---|
[f619de7] | 324 | # type: (ModelInfo, KernelModel, KernelModel) -> None |
---|
[146793b] | 325 | #: Combined info plock for the product model |
---|
[72a081d] | 326 | self.info = model_info |
---|
[146793b] | 327 | #: Form factor modelling individual particles. |
---|
[17bbadd] | 328 | self.P = P |
---|
[146793b] | 329 | #: Structure factor modelling interaction between particles. |
---|
[17bbadd] | 330 | self.S = S |
---|
[c036ddb] | 331 | |
---|
[146793b] | 332 | #: Model precision. This is not really relevant, since it is the |
---|
| 333 | #: individual P and S models that control the effective dtype, |
---|
| 334 | #: converting the q-vectors to the correct type when the kernels |
---|
| 335 | #: for each are created. Ideally this should be set to the more |
---|
| 336 | #: precise type to avoid loss of precision, but precision in q is |
---|
| 337 | #: not critical (single is good enough for our purposes), so it just |
---|
| 338 | #: uses the precision of the form factor. |
---|
| 339 | self.dtype = P.dtype # type: np.dtype |
---|
[17bbadd] | 340 | |
---|
[6dc78e4] | 341 | def make_kernel(self, q_vectors): |
---|
[f619de7] | 342 | # type: (List[np.ndarray]) -> Kernel |
---|
[17bbadd] | 343 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
| 344 | # are only needed once. It would mess up modularity quite a bit to |
---|
| 345 | # handle this optimally, especially since there are many cases where |
---|
| 346 | # separate q vectors are needed (e.g., form in python and structure |
---|
| 347 | # in opencl; or both in opencl, but one in single precision and the |
---|
| 348 | # other in double precision). |
---|
[c036ddb] | 349 | |
---|
[f619de7] | 350 | p_kernel = self.P.make_kernel(q_vectors) |
---|
| 351 | s_kernel = self.S.make_kernel(q_vectors) |
---|
[17bbadd] | 352 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
[b297ba9] | 353 | make_kernel.__doc__ = KernelModel.make_kernel.__doc__ |
---|
[17bbadd] | 354 | |
---|
| 355 | def release(self): |
---|
[f619de7] | 356 | # type: (None) -> None |
---|
[17bbadd] | 357 | """ |
---|
| 358 | Free resources associated with the model. |
---|
| 359 | """ |
---|
| 360 | self.P.release() |
---|
| 361 | self.S.release() |
---|
| 362 | |
---|
| 363 | |
---|
[f619de7] | 364 | class ProductKernel(Kernel): |
---|
[b297ba9] | 365 | """ |
---|
| 366 | Instantiated kernel for product model. |
---|
| 367 | """ |
---|
[17bbadd] | 368 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
[f619de7] | 369 | # type: (ModelInfo, Kernel, Kernel) -> None |
---|
[17bbadd] | 370 | self.info = model_info |
---|
| 371 | self.p_kernel = p_kernel |
---|
| 372 | self.s_kernel = s_kernel |
---|
[6dc78e4] | 373 | self.dtype = p_kernel.dtype |
---|
| 374 | self.results = [] # type: List[np.ndarray] |
---|
| 375 | |
---|
[b2f0e5d] | 376 | # Find index of volfraction parameter in parameter list |
---|
| 377 | for k, p in enumerate(model_info.parameters.call_parameters): |
---|
| 378 | if p.id == VOLFRAC_ID: |
---|
| 379 | self._volfrac_index = k |
---|
| 380 | break |
---|
| 381 | else: |
---|
| 382 | raise RuntimeError("no %s parameter in %s"%(VOLFRAC_ID, self)) |
---|
[e44432d] | 383 | |
---|
[6dc78e4] | 384 | p_info, s_info = self.info.composition[1] |
---|
[e44432d] | 385 | p_npars = p_info.parameters.npars |
---|
| 386 | s_npars = s_info.parameters.npars |
---|
| 387 | |
---|
| 388 | have_beta_mode = p_info.have_Fq |
---|
[b2f0e5d] | 389 | have_er_mode = p_info.radius_effective_modes is not None |
---|
| 390 | volfrac_in_p = self._volfrac_index < p_npars + 2 # scale & background |
---|
| 391 | |
---|
| 392 | # Slices into the details length/offset structure for P@S. |
---|
| 393 | # Made complicated by the possibly missing volfraction in S. |
---|
| 394 | self._p_detail_slice = slice(0, p_npars) |
---|
| 395 | self._s_detail_slice = slice(p_npars, p_npars+s_npars-volfrac_in_p) |
---|
| 396 | self._volfrac_in_p = volfrac_in_p |
---|
| 397 | |
---|
| 398 | # P block from data vector, without scale and background |
---|
| 399 | first_p = 2 |
---|
| 400 | last_p = p_npars + 2 |
---|
| 401 | self._p_value_slice = slice(first_p, last_p) |
---|
| 402 | |
---|
| 403 | # radius_effective is the first parameter in S from the data vector. |
---|
| 404 | self._er_index = last_p |
---|
| 405 | |
---|
| 406 | # S block from data vector, without scale, background, volfrac or er. |
---|
| 407 | first_s = last_p + 2 - volfrac_in_p |
---|
| 408 | last_s = first_s + s_npars - 2 |
---|
| 409 | self._s_value_slice = slice(first_s, last_s) |
---|
| 410 | |
---|
| 411 | # S distribution block in S data vector starts after all S values |
---|
| 412 | self._s_dist_slice = slice(2 + s_npars, None) |
---|
| 413 | |
---|
| 414 | # structure_factor_mode is the first parameter after P and S. Skip |
---|
| 415 | # 2 for scale and background, and subtract 1 in case there is no |
---|
| 416 | # volfraction in S. |
---|
| 417 | self._beta_mode_index = last_s if have_beta_mode else 0 |
---|
| 418 | |
---|
| 419 | # radius_effective_mode is the second parameter after P and S |
---|
| 420 | # unless structure_factor_mode isn't available, in which case it |
---|
| 421 | # is first. |
---|
| 422 | self._er_mode_index = last_s + have_beta_mode if have_er_mode else 0 |
---|
| 423 | |
---|
| 424 | # Magnetic parameters are after everything else. If they exist, |
---|
| 425 | # they will only be for form factor P, not structure factor S. |
---|
| 426 | first_mag = last_s + have_beta_mode + have_er_mode |
---|
| 427 | mag_pars = 3*p_info.parameters.nmagnetic |
---|
| 428 | last_mag = first_mag + (mag_pars + 3 if mag_pars else 0) |
---|
| 429 | self._magentic_slice = slice(first_mag, last_mag) |
---|
[e44432d] | 430 | |
---|
[b2f0e5d] | 431 | def Iq(self, call_details, values, cutoff, magnetic): |
---|
| 432 | # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray |
---|
| 433 | p_info, s_info = self.info.composition[1] |
---|
| 434 | |
---|
| 435 | # Retrieve values from the data vector |
---|
[e44432d] | 436 | scale, background = values[0], values[1] |
---|
[b2f0e5d] | 437 | volfrac = values[self._volfrac_index] |
---|
| 438 | er_mode = (int(values[self._er_mode_index]) |
---|
| 439 | if self._er_mode_index > 0 else 0) |
---|
| 440 | beta_mode = (values[self._beta_mode_index] > 0 |
---|
| 441 | if self._beta_mode_index > 0 else False) |
---|
[6dc78e4] | 442 | |
---|
| 443 | nvalues = self.info.parameters.nvalues |
---|
| 444 | nweights = call_details.num_weights |
---|
| 445 | weights = values[nvalues:nvalues + 2*nweights] |
---|
| 446 | |
---|
[b2f0e5d] | 447 | # Can't do 2d and beta_mode just yet |
---|
| 448 | if beta_mode and self.p_kernel.dim == '2d': |
---|
| 449 | raise NotImplementedError("beta not yet supported for 2D") |
---|
| 450 | |
---|
[6dc78e4] | 451 | # Construct the calling parameters for P. |
---|
[b2f0e5d] | 452 | p_length = call_details.length[self._p_detail_slice] |
---|
| 453 | p_offset = call_details.offset[self._p_detail_slice] |
---|
[6dc78e4] | 454 | p_details = make_details(p_info, p_length, p_offset, nweights) |
---|
[e44432d] | 455 | p_values = [ |
---|
| 456 | [1., 0.], # scale=1, background=0, |
---|
[b2f0e5d] | 457 | values[self._p_value_slice], |
---|
| 458 | values[self._magentic_slice], |
---|
[e44432d] | 459 | weights] |
---|
[6dc78e4] | 460 | spacer = (32 - sum(len(v) for v in p_values)%32)%32 |
---|
| 461 | p_values.append([0.]*spacer) |
---|
| 462 | p_values = np.hstack(p_values).astype(self.p_kernel.dtype) |
---|
| 463 | |
---|
[b2f0e5d] | 464 | # Call the form factor kernel to compute <F> and <F^2>. |
---|
| 465 | # If the model doesn't support Fq the returned <F> will be None. |
---|
| 466 | F, Fsq, radius_effective, shell_volume, volume_ratio \ |
---|
| 467 | = self.p_kernel.Fq(p_details, p_values, cutoff, magnetic, er_mode) |
---|
| 468 | |
---|
| 469 | # TODO: async call to the GPU |
---|
| 470 | |
---|
[6dc78e4] | 471 | # Construct the calling parameters for S. |
---|
[b2f0e5d] | 472 | s_length = call_details.length[self._s_detail_slice] |
---|
| 473 | s_offset = call_details.offset[self._s_detail_slice] |
---|
| 474 | if self._volfrac_in_p: |
---|
| 475 | # Volfrac is in P and missing from S so insert a slot for it. Say |
---|
| 476 | # the distribution is length 1 and use the slot for volfraction |
---|
| 477 | # from the P distribution. |
---|
| 478 | s_length = np.insert(s_length, 1, 1) |
---|
| 479 | s_offset = np.insert(s_offset, 1, p_offset[self._volfrac_index - 2]) |
---|
| 480 | if er_mode > 0: |
---|
| 481 | # If effective_radius comes from P, make sure it is monodisperse. |
---|
| 482 | # Weight is set to 1 later, after the value array is created |
---|
[e44432d] | 483 | s_length[0] = 1 |
---|
| 484 | s_details = make_details(s_info, s_length, s_offset, nweights) |
---|
[9951a86] | 485 | s_values = [ |
---|
[b2f0e5d] | 486 | [1., # scale=1 |
---|
| 487 | 0., # background=0, |
---|
| 488 | values[self._er_index], # S.radius_effective; may be replaced by P |
---|
| 489 | 0.], # volfraction; will be replaced by volfrac * volume_ratio |
---|
| 490 | # followed by S parameters after effective_radius and volfraction |
---|
| 491 | values[self._s_value_slice], |
---|
[6e7ba14] | 492 | weights, |
---|
[9951a86] | 493 | ] |
---|
[6dc78e4] | 494 | spacer = (32 - sum(len(v) for v in s_values)%32)%32 |
---|
| 495 | s_values.append([0.]*spacer) |
---|
| 496 | s_values = np.hstack(s_values).astype(self.s_kernel.dtype) |
---|
[17bbadd] | 497 | |
---|
[e44432d] | 498 | # Plug R_eff from the form factor into structure factor parameters |
---|
| 499 | # and scale volume fraction by form:shell volume ratio. These changes |
---|
| 500 | # needs to be both in the initial value slot as well as the |
---|
| 501 | # polydispersity distribution slot in the values array due to |
---|
| 502 | # implementation details in kernel_iq.c. |
---|
[8795b6f] | 503 | #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g" |
---|
| 504 | # % (radius_type, radius_effective, volfrac, volume_ratio)) |
---|
[b2f0e5d] | 505 | s_dist = s_values[self._s_dist_slice] |
---|
| 506 | if er_mode > 0: |
---|
[39a06c9] | 507 | # set the value to the model R_eff and set the weight to 1 |
---|
[b2f0e5d] | 508 | s_values[2] = s_dist[s_offset[0]] = radius_effective |
---|
| 509 | s_dist[s_offset[0]+nweights] = 1.0 |
---|
| 510 | s_values[3] = s_dist[s_offset[1]] = volfrac*volume_ratio |
---|
| 511 | s_dist[s_offset[1]+nweights] = 1.0 |
---|
[e44432d] | 512 | |
---|
[b2f0e5d] | 513 | # Call the structure factor kernel to compute S. |
---|
| 514 | S = self.s_kernel.Iq(s_details, s_values, cutoff, False) |
---|
| 515 | #print("P", Fsq[:10]) |
---|
| 516 | #print("S", S[:10]) |
---|
| 517 | #print(radius_effective, volfrac*volume_ratio) |
---|
[e44432d] | 518 | |
---|
| 519 | # Combine form factor and structure factor |
---|
[b2f0e5d] | 520 | #print("beta", beta_mode, F, Fsq, S) |
---|
| 521 | PS = Fsq + F**2*(S-1) if beta_mode else Fsq*S |
---|
| 522 | |
---|
| 523 | # Determine overall scale factor. Hollow shapes are weighted by |
---|
| 524 | # shell_volume, so that is needed for number density estimation. |
---|
| 525 | # For solid shapes we can use shell_volume as well since it is |
---|
| 526 | # equal to form volume. If P already has a volfraction parameter, |
---|
| 527 | # then assume that it is already on absolute scale, and don't |
---|
| 528 | # include volfrac in the combined_scale. |
---|
| 529 | combined_scale = scale*(volfrac if not self._volfrac_in_p else 1.0) |
---|
| 530 | final_result = combined_scale/shell_volume*PS + background |
---|
[e44432d] | 531 | |
---|
| 532 | # Capture intermediate values so user can see them. These are |
---|
| 533 | # returned as a lazy evaluator since they are only needed in the |
---|
| 534 | # GUI, and not for each evaluation during a fit. |
---|
| 535 | # TODO: return the results structure with the final results |
---|
| 536 | # That way the model calcs are idempotent. Further, we can |
---|
| 537 | # generalize intermediates to various other model types if we put it |
---|
| 538 | # kernel calling interface. Could do this as an "optional" |
---|
| 539 | # return value in the caller, though in that case we could return |
---|
| 540 | # the results directly rather than through a lazy evaluator. |
---|
| 541 | self.results = lambda: _intermediates( |
---|
[b2f0e5d] | 542 | F, Fsq, S, combined_scale, radius_effective, beta_mode) |
---|
[d3ffeb7] | 543 | |
---|
[01c8d9e] | 544 | return final_result |
---|
[17bbadd] | 545 | |
---|
[b297ba9] | 546 | Iq.__doc__ = Kernel.Iq.__doc__ |
---|
| 547 | __call__ = Iq |
---|
| 548 | |
---|
[17bbadd] | 549 | def release(self): |
---|
[f619de7] | 550 | # type: () -> None |
---|
[b297ba9] | 551 | """Free resources associated with the kernel.""" |
---|
[17bbadd] | 552 | self.p_kernel.release() |
---|
[f619de7] | 553 | self.s_kernel.release() |
---|