source: sasmodels/sasmodels/product.py @ b2f0e5d

ticket-1257-vesicle-product
Last change on this file since b2f0e5d was b2f0e5d, checked in by Paul Kienzle <pkienzle@…>, 5 months ago

Fix handling of volfraction in vesicle@hardsphere. Refs 1257.

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