source: sasmodels/sasmodels/product.py @ e44432d

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e44432d was e44432d, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

support hollow models in structure factor calculations

  • Property mode set to 100644
File size: 16.6 KB
Line 
1"""
2Product model
3-------------
4
5The product model multiplies the structure factor by the form factor,
6modulated by the effective radius of the form.  The resulting model
7has a attributes of both the model description (with parameters, etc.)
8and the module evaluator (with call, release, etc.).
9
10To use it, first load form factor P and structure factor S, then create
11*make_product_info(P, S)*.
12"""
13from __future__ import print_function, division
14
15from collections import OrderedDict
16
17from copy import copy
18import numpy as np  # type: ignore
19
20from .modelinfo import ParameterTable, ModelInfo, Parameter, parse_parameter
21from .kernel import KernelModel, Kernel
22from .details import make_details, dispersion_mesh
23
24# pylint: disable=unused-import
25try:
26    from typing import Tuple, Callable, Union
27except ImportError:
28    pass
29else:
30    from .modelinfo import ParameterSet
31# pylint: enable=unused-import
32
33# TODO: make estimates available to constraints
34#ESTIMATED_PARAMETERS = [
35#    ["est_radius_effective", "A", 0.0, [0, np.inf], "", "Estimated effective radius"],
36#    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"],
37#]
38
39RADIUS_ID = "radius_effective"
40VOLFRAC_ID = "volfraction"
41def make_extra_pars(p_info):
42    pars = []
43    if p_info.have_Fq:
44        par = parse_parameter(
45                "structure_factor_mode",
46                "",
47                0,
48                [["P*S","P*(1+beta*(S-1))"]],
49                "",
50                "Structure factor calculation")
51        pars.append(par)
52    if p_info.effective_radius_type is not None:
53        par = parse_parameter(
54                "radius_effective_mode",
55                "",
56                0,
57                [["unconstrained"] + p_info.effective_radius_type],
58                "",
59                "Effective radius calculation")
60        pars.append(par)
61    return pars
62
63# TODO: core_shell_sphere model has suppressed the volume ratio calculation
64# revert it after making VR and ER available at run time as constraints.
65def make_product_info(p_info, s_info):
66    # type: (ModelInfo, ModelInfo) -> ModelInfo
67    """
68    Create info block for product model.
69    """
70    # Make sure effective radius is the first parameter and
71    # make sure volume fraction is the second parameter of the
72    # structure factor calculator.  Structure factors should not
73    # have any magnetic parameters
74    if not len(s_info.parameters.kernel_parameters) >= 2:
75        raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID))
76    if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID:
77        raise TypeError("S needs {} as first parameter".format(RADIUS_ID))
78    if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID:
79        raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID))
80    if not s_info.parameters.magnetism_index == []:
81        raise TypeError("S should not have SLD parameters")
82    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters
83    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
84
85    # Create list of parameters for the combined model.  If there
86    # are any names in P that overlap with those in S, modify the name in S
87    # to distinguish it.
88    p_set = set(p.id for p in p_pars.kernel_parameters)
89    s_list = [(_tag_parameter(par) if par.id in p_set else par)
90              for par in s_pars.kernel_parameters]
91    # Check if still a collision after renaming.  This could happen if for
92    # example S has volfrac and P has both volfrac and volfrac_S.
93    if any(p.id in p_set for p in s_list):
94        raise TypeError("name collision: P has P.name and P.name_S while S has S.name")
95
96    # make sure effective radius is not a polydisperse parameter in product
97    s_list[0] = copy(s_list[0])
98    s_list[0].polydisperse = False
99
100    translate_name = dict((old.id, new.id) for old, new
101                          in zip(s_pars.kernel_parameters, s_list))
102    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info)
103    parameters = ParameterTable(combined_pars)
104    parameters.max_pd = p_pars.max_pd + s_pars.max_pd
105    def random():
106        combined_pars = p_info.random()
107        s_names = set(par.id for par in s_pars.kernel_parameters)
108        combined_pars.update((translate_name[k], v)
109                             for k, v in s_info.random().items()
110                             if k in s_names)
111        return combined_pars
112
113    model_info = ModelInfo()
114    model_info.id = '@'.join((p_id, s_id))
115    model_info.name = '@'.join((p_name, s_name))
116    model_info.filename = None
117    model_info.title = 'Product of %s and %s'%(p_name, s_name)
118    model_info.description = model_info.title
119    model_info.docs = model_info.title
120    model_info.category = "custom"
121    model_info.parameters = parameters
122    model_info.random = random
123    #model_info.single = p_info.single and s_info.single
124    model_info.structure_factor = False
125    model_info.variant_info = None
126    #model_info.tests = []
127    #model_info.source = []
128    # Iq, Iqxy, form_volume, ER, VR and sesans
129    # Remember the component info blocks so we can build the model
130    model_info.composition = ('product', [p_info, s_info])
131    model_info.control = p_info.control
132    model_info.hidden = p_info.hidden
133    if getattr(p_info, 'profile', None) is not None:
134        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
135        def profile(**kwargs):
136            # extract the profile args
137            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
138            return p_info.profile(**kwargs)
139    else:
140        profile = None
141    model_info.profile = profile
142    model_info.profile_axes = p_info.profile_axes
143
144    # TODO: delegate random to p_info, s_info
145    #model_info.random = lambda: {}
146
147    ## Show the parameter table
148    #from .compare import get_pars, parlist
149    #print("==== %s ====="%model_info.name)
150    #values = get_pars(model_info)
151    #print(parlist(model_info, values, is2d=True))
152    return model_info
153
154def _tag_parameter(par):
155    """
156    Tag the parameter name with _S to indicate that the parameter comes from
157    the structure factor parameter set.  This is only necessary if the
158    form factor model includes a parameter of the same name as a parameter
159    in the structure factor.
160    """
161    par = copy(par)
162    # Protect against a vector parameter in S by appending the vector length
163    # to the renamed parameter.  Note: haven't tested this since no existing
164    # structure factor models contain vector parameters.
165    vector_length = par.name[len(par.id):]
166    par.id = par.id + "_S"
167    par.name = par.id + vector_length
168    return par
169
170def _intermediates(
171        F1,               # type: np.ndarray
172        F2,               # type: np.ndarray
173        S,                # type: np.ndarray
174        scale,            # type: float
175        effective_radius, # type: float
176        beta_mode,        # type: bool
177        ):
178    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]]
179    """
180    Returns intermediate results for beta approximation-enabled product.
181    The result may be an array or a float.
182    """
183    if beta_mode:
184        # TODO: 1. include calculated Q vector
185        # TODO: 2. consider implications if there are intermediate results in P(Q)
186        parts = OrderedDict((
187            ("P(Q)", scale*F2),
188            ("S(Q)", S),
189            ("beta(Q)", F1**2 / F2),
190            ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),
191            ("effective_radius", effective_radius),
192            # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),
193        ))
194    else:
195        parts = OrderedDict((
196            ("P(Q)", scale*F2),
197            ("S(Q)", S),
198            ("effective_radius", effective_radius),
199        ))
200    return parts
201
202class ProductModel(KernelModel):
203    def __init__(self, model_info, P, S):
204        # type: (ModelInfo, KernelModel, KernelModel) -> None
205        #: Combined info plock for the product model
206        self.info = model_info
207        #: Form factor modelling individual particles.
208        self.P = P
209        #: Structure factor modelling interaction between particles.
210        self.S = S
211
212        #: Model precision. This is not really relevant, since it is the
213        #: individual P and S models that control the effective dtype,
214        #: converting the q-vectors to the correct type when the kernels
215        #: for each are created. Ideally this should be set to the more
216        #: precise type to avoid loss of precision, but precision in q is
217        #: not critical (single is good enough for our purposes), so it just
218        #: uses the precision of the form factor.
219        self.dtype = P.dtype  # type: np.dtype
220
221    def make_kernel(self, q_vectors):
222        # type: (List[np.ndarray]) -> Kernel
223        # Note: may be sending the q_vectors to the GPU twice even though they
224        # are only needed once.  It would mess up modularity quite a bit to
225        # handle this optimally, especially since there are many cases where
226        # separate q vectors are needed (e.g., form in python and structure
227        # in opencl; or both in opencl, but one in single precision and the
228        # other in double precision).
229
230        p_kernel = self.P.make_kernel(q_vectors)
231        s_kernel = self.S.make_kernel(q_vectors)
232        return ProductKernel(self.info, p_kernel, s_kernel)
233
234    def release(self):
235        # type: (None) -> None
236        """
237        Free resources associated with the model.
238        """
239        self.P.release()
240        self.S.release()
241
242
243class ProductKernel(Kernel):
244    def __init__(self, model_info, p_kernel, s_kernel):
245        # type: (ModelInfo, Kernel, Kernel) -> None
246        self.info = model_info
247        self.p_kernel = p_kernel
248        self.s_kernel = s_kernel
249        self.dtype = p_kernel.dtype
250        self.results = []  # type: List[np.ndarray]
251
252    def __call__(self, call_details, values, cutoff, magnetic):
253        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
254
255        p_info, s_info = self.info.composition[1]
256        p_npars = p_info.parameters.npars
257        p_length = call_details.length[:p_npars]
258        p_offset = call_details.offset[:p_npars]
259        s_npars = s_info.parameters.npars
260        s_length = call_details.length[p_npars:p_npars+s_npars]
261        s_offset = call_details.offset[p_npars:p_npars+s_npars]
262
263        # Beta mode parameter is the first parameter after P and S parameters
264        have_beta_mode = p_info.have_Fq
265        beta_mode_offset = 2+p_npars+s_npars
266        beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False
267        if beta_mode and self.p_kernel.dim== '2d':
268            raise NotImplementedError("beta not yet supported for 2D")
269
270        # R_eff type parameter is the second parameter after P and S parameters
271        # unless the model doesn't support beta mode, in which case it is first
272        have_radius_type = p_info.effective_radius_type is not None
273        radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0)
274        radius_type = int(values[radius_type_offset]) if have_radius_type else 0
275
276        # Retrieve the volume fraction, which is the second of the
277        # 'S' parameters in the parameter list, or 2+np in 0-origin,
278        # as well as the scale and background.
279        volfrac = values[3+p_npars]
280        scale, background = values[0], values[1]
281
282        # if there are magnetic parameters, they will only be on the
283        # form factor P, not the structure factor S.
284        nmagnetic = len(self.info.parameters.magnetism_index)
285        if nmagnetic:
286            spin_index = self.info.parameters.npars + 2
287            magnetism = values[spin_index: spin_index+3+3*nmagnetic]
288        else:
289            magnetism = []
290        nvalues = self.info.parameters.nvalues
291        nweights = call_details.num_weights
292        weights = values[nvalues:nvalues + 2*nweights]
293
294        # Construct the calling parameters for P.
295        p_details = make_details(p_info, p_length, p_offset, nweights)
296        p_values = [
297            [1., 0.], # scale=1, background=0,
298            values[2:2+p_npars],
299            magnetism,
300            weights]
301        spacer = (32 - sum(len(v) for v in p_values)%32)%32
302        p_values.append([0.]*spacer)
303        p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
304
305        # Construct the calling parameters for S.
306        if radius_type > 0:
307            # If R_eff comes from form factor, make sure it is monodisperse.
308            # weight is set to 1 later, after the value array is created
309            s_length[0] = 1
310        s_details = make_details(s_info, s_length, s_offset, nweights)
311        s_values = [
312            [1., 0.], # scale=1, background=0,
313            values[2+p_npars:2+p_npars+s_npars],
314            weights,
315        ]
316        spacer = (32 - sum(len(v) for v in s_values)%32)%32
317        s_values.append([0.]*spacer)
318        s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
319
320        # Call the form factor kernel to compute <F> and <F^2>.
321        # If the model doesn't support Fq the returned <F> will be None.
322        F1, F2, effective_radius, shell_volume, volume_ratio = self.p_kernel.Fq(
323            p_details, p_values, cutoff, magnetic, radius_type)
324
325        # Call the structure factor kernel to compute S.
326        # Plug R_eff from the form factor into structure factor parameters
327        # and scale volume fraction by form:shell volume ratio. These changes
328        # needs to be both in the initial value slot as well as the
329        # polydispersity distribution slot in the values array due to
330        # implementation details in kernel_iq.c.
331        #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g"%(radius_type, effective_radius, volfrac, volume_ratio))
332        if radius_type > 0:
333            # set the value to the model ER and set the weight to 1
334            s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
335            s_values[2+s_npars+s_offset[0]+nweights] = 1.0
336        s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio
337        S = self.s_kernel.Iq(s_details, s_values, cutoff, False)
338
339        # Determine overall scale factor. Hollow shapes are weighted by
340        # shell_volume, so that is needed for volume normalization.  For
341        # solid shapes we can use shell_volume as well since it is equal
342        # to form volume.
343        combined_scale = scale*volfrac/shell_volume
344
345        # Combine form factor and structure factor
346        #print("beta", beta_mode, F1, F2, S)
347        PS = F2 + F1**2*(S-1) if beta_mode else F2*S
348        final_result = combined_scale*PS + background
349
350        # Capture intermediate values so user can see them.  These are
351        # returned as a lazy evaluator since they are only needed in the
352        # GUI, and not for each evaluation during a fit.
353        # TODO: return the results structure with the final results
354        # That way the model calcs are idempotent. Further, we can
355        # generalize intermediates to various other model types if we put it
356        # kernel calling interface.  Could do this as an "optional"
357        # return value in the caller, though in that case we could return
358        # the results directly rather than through a lazy evaluator.
359        self.results = lambda: _intermediates(
360            F1, F2, S, combined_scale, effective_radius, beta_mode)
361
362        return final_result
363
364    def release(self):
365        # type: () -> None
366        self.p_kernel.release()
367        self.s_kernel.release()
368
369
370def calc_er_vr(model_info, call_details, values):
371    # type: (ModelInfo, ParameterSet) -> Tuple[float, float]
372
373    if model_info.ER is None and model_info.VR is None:
374        return 1.0, 1.0
375
376    nvalues = model_info.parameters.nvalues
377    value = values[nvalues:nvalues + call_details.num_weights]
378    weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
379    npars = model_info.parameters.npars
380    # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the
381    # dispersion mesh code doesn't actually care about the nominal parameter
382    # value, p, so set it to None.
383    pairs = [(None, value[offset:offset+length], weight[offset:offset+length])
384             for p, offset, length
385             in zip(model_info.parameters.call_parameters[2:2+npars],
386                    call_details.offset,
387                    call_details.length)
388             if p.type == 'volume']
389    value, weight = dispersion_mesh(model_info, pairs)
390
391    if model_info.ER is not None:
392        individual_radii = model_info.ER(*value)
393        radius_effective = np.sum(weight*individual_radii) / np.sum(weight)
394    else:
395        radius_effective = 1.0
396
397    if model_info.VR is not None:
398        whole, part = model_info.VR(*value)
399        volume_ratio = np.sum(weight*part)/np.sum(weight*whole)
400    else:
401        volume_ratio = 1.0
402
403    return radius_effective, volume_ratio
Note: See TracBrowser for help on using the repository browser.