source: sasmodels/sasmodels/product.py @ d8e81f7

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since d8e81f7 was d8e81f7, checked in by ajj, 10 months ago

Merge branch 'ticket_822_v5_unit_tests' of https://github.com/SasView/sasmodels into ticket_822_v5_unit_tests

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