source: sasmodels/sasmodels/product.py @ d277229

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

fix R_eff support infrastructure so more tests pass

  • Property mode set to 100644
File size: 15.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"""
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
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 = Parameter("structure_factor_mode", ["P*S","P*(1+beta*(S-1))"], 0, None, "",
45                        "Structure factor calculation")
46        pars.append(par)
47    if p_info.effective_radius_type is not None:
48        par = Parameter("radius_effective_mode", ["unconstrained"] + p_info.effective_radius_type, 0, None, "",
49                        "Effective radius calculation")
50        pars.append(par)
51    return pars
52
53# TODO: core_shell_sphere model has suppressed the volume ratio calculation
54# revert it after making VR and ER available at run time as constraints.
55def make_product_info(p_info, s_info):
56    # type: (ModelInfo, ModelInfo) -> ModelInfo
57    """
58    Create info block for product model.
59    """
60    # Make sure effective radius is the first parameter and
61    # make sure volume fraction is the second parameter of the
62    # structure factor calculator.  Structure factors should not
63    # have any magnetic parameters
64    if not len(s_info.parameters.kernel_parameters) >= 2:
65        raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID))
66    if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID:
67        raise TypeError("S needs {} as first parameter".format(RADIUS_ID))
68    if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID:
69        raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID))
70    if not s_info.parameters.magnetism_index == []:
71        raise TypeError("S should not have SLD parameters")
72    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters
73    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
74
75    # Create list of parameters for the combined model.  If there
76    # are any names in P that overlap with those in S, modify the name in S
77    # to distinguish it.
78    p_set = set(p.id for p in p_pars.kernel_parameters)
79    s_list = [(_tag_parameter(par) if par.id in p_set else par)
80              for par in s_pars.kernel_parameters]
81    # Check if still a collision after renaming.  This could happen if for
82    # example S has volfrac and P has both volfrac and volfrac_S.
83    if any(p.id in p_set for p in s_list):
84        raise TypeError("name collision: P has P.name and P.name_S while S has S.name")
85
86    # make sure effective radius is not a polydisperse parameter in product
87    s_list[0] = copy(s_list[0])
88    s_list[0].polydisperse = False
89
90    translate_name = dict((old.id, new.id) for old, new
91                          in zip(s_pars.kernel_parameters, s_list))
92    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info)
93    parameters = ParameterTable(combined_pars)
94    parameters.max_pd = p_pars.max_pd + s_pars.max_pd
95    def random():
96        combined_pars = p_info.random()
97        s_names = set(par.id for par in s_pars.kernel_parameters)
98        combined_pars.update((translate_name[k], v)
99                             for k, v in s_info.random().items()
100                             if k in s_names)
101        return combined_pars
102
103    model_info = ModelInfo()
104    model_info.id = '@'.join((p_id, s_id))
105    model_info.name = '@'.join((p_name, s_name))
106    model_info.filename = None
107    model_info.title = 'Product of %s and %s'%(p_name, s_name)
108    model_info.description = model_info.title
109    model_info.docs = model_info.title
110    model_info.category = "custom"
111    model_info.parameters = parameters
112    model_info.random = random
113    #model_info.single = p_info.single and s_info.single
114    model_info.structure_factor = False
115    model_info.variant_info = None
116    #model_info.tests = []
117    #model_info.source = []
118    # Iq, Iqxy, form_volume, ER, VR and sesans
119    # Remember the component info blocks so we can build the model
120    model_info.composition = ('product', [p_info, s_info])
121    model_info.control = p_info.control
122    model_info.hidden = p_info.hidden
123    if getattr(p_info, 'profile', None) is not None:
124        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
125        def profile(**kwargs):
126            # extract the profile args
127            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
128            return p_info.profile(**kwargs)
129    else:
130        profile = None
131    model_info.profile = profile
132    model_info.profile_axes = p_info.profile_axes
133
134    # TODO: delegate random to p_info, s_info
135    #model_info.random = lambda: {}
136
137    ## Show the parameter table
138    #from .compare import get_pars, parlist
139    #print("==== %s ====="%model_info.name)
140    #values = get_pars(model_info)
141    #print(parlist(model_info, values, is2d=True))
142    return model_info
143
144def _tag_parameter(par):
145    """
146    Tag the parameter name with _S to indicate that the parameter comes from
147    the structure factor parameter set.  This is only necessary if the
148    form factor model includes a parameter of the same name as a parameter
149    in the structure factor.
150    """
151    par = copy(par)
152    # Protect against a vector parameter in S by appending the vector length
153    # to the renamed parameter.  Note: haven't tested this since no existing
154    # structure factor models contain vector parameters.
155    vector_length = par.name[len(par.id):]
156    par.id = par.id + "_S"
157    par.name = par.id + vector_length
158    return par
159
160def _intermediates(P, S, effective_radius):
161    # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray]
162    """
163    Returns intermediate results for standard product (P(Q)*S(Q))
164    """
165    return OrderedDict((
166        ("P(Q)", P),
167        ("S(Q)", S),
168        ("effective_radius", effective_radius),
169    ))
170
171def _intermediates_beta(F1, F2, S, scale, bg, effective_radius):
172    # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float) -> OrderedDict[str, Union[np.ndarray, float]]
173    """
174    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float.
175    """
176    # TODO: 1. include calculated Q vector
177    # TODO: 2. consider implications if there are intermediate results in P(Q)
178    return OrderedDict((
179        ("P(Q)", scale*F2),
180        ("S(Q)", S),
181        ("beta(Q)", F1**2 / F2),
182        ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),
183        ("effective_radius", effective_radius),
184        # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),
185    ))
186
187class ProductModel(KernelModel):
188    def __init__(self, model_info, P, S):
189        # type: (ModelInfo, KernelModel, KernelModel) -> None
190        #: Combined info plock for the product model
191        self.info = model_info
192        #: Form factor modelling individual particles.
193        self.P = P
194        #: Structure factor modelling interaction between particles.
195        self.S = S
196
197        #: Model precision. This is not really relevant, since it is the
198        #: individual P and S models that control the effective dtype,
199        #: converting the q-vectors to the correct type when the kernels
200        #: for each are created. Ideally this should be set to the more
201        #: precise type to avoid loss of precision, but precision in q is
202        #: not critical (single is good enough for our purposes), so it just
203        #: uses the precision of the form factor.
204        self.dtype = P.dtype  # type: np.dtype
205
206    def make_kernel(self, q_vectors):
207        # type: (List[np.ndarray]) -> Kernel
208        # Note: may be sending the q_vectors to the GPU twice even though they
209        # are only needed once.  It would mess up modularity quite a bit to
210        # handle this optimally, especially since there are many cases where
211        # separate q vectors are needed (e.g., form in python and structure
212        # in opencl; or both in opencl, but one in single precision and the
213        # other in double precision).
214
215        p_kernel = self.P.make_kernel(q_vectors)
216        s_kernel = self.S.make_kernel(q_vectors)
217        return ProductKernel(self.info, p_kernel, s_kernel)
218
219    def release(self):
220        # type: (None) -> None
221        """
222        Free resources associated with the model.
223        """
224        self.P.release()
225        self.S.release()
226
227
228class ProductKernel(Kernel):
229    def __init__(self, model_info, p_kernel, s_kernel):
230        # type: (ModelInfo, Kernel, Kernel) -> None
231        self.info = model_info
232        self.p_kernel = p_kernel
233        self.s_kernel = s_kernel
234        self.dtype = p_kernel.dtype
235        self.results = []  # type: List[np.ndarray]
236
237    def __call__(self, call_details, values, cutoff, magnetic):
238        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
239        p_info, s_info = self.info.composition[1]
240
241        # if there are magnetic parameters, they will only be on the
242        # form factor P, not the structure factor S.
243        nmagnetic = len(self.info.parameters.magnetism_index)
244        if nmagnetic:
245            spin_index = self.info.parameters.npars + 2
246            magnetism = values[spin_index: spin_index+3+3*nmagnetic]
247        else:
248            magnetism = []
249        nvalues = self.info.parameters.nvalues
250        nweights = call_details.num_weights
251        weights = values[nvalues:nvalues + 2*nweights]
252
253        # Construct the calling parameters for P.
254        p_npars = p_info.parameters.npars
255        p_length = call_details.length[:p_npars]
256        p_offset = call_details.offset[:p_npars]
257        p_details = make_details(p_info, p_length, p_offset, nweights)
258
259        # Set p scale to the volume fraction in s, which is the first of the
260        # 'S' parameters in the parameter list, or 2+np in 0-origin.
261        volfrac = values[2+p_npars]
262        p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights]
263        spacer = (32 - sum(len(v) for v in p_values)%32)%32
264        p_values.append([0.]*spacer)
265        p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
266
267        # Construct the calling parameters for S.
268        s_npars = s_info.parameters.npars
269        s_length = call_details.length[p_npars:p_npars+s_npars]
270        s_offset = call_details.offset[p_npars:p_npars+s_npars]
271        s_length = np.hstack((1, s_length))
272        s_offset = np.hstack((nweights, s_offset))
273        s_details = make_details(s_info, s_length, s_offset, nweights+1)
274        s_values = [
275            # scale=1, background=0,
276            [1., 0.],
277            values[2+p_npars:2+p_npars+s_npars],
278            weights,
279        ]
280        spacer = (32 - sum(len(v) for v in s_values)%32)%32
281        s_values.append([0.]*spacer)
282        s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
283
284        # beta mode is the first parameter after the structure factor pars
285        extra_offset = 2+p_npars+s_npars
286        if p_info.have_Fq:
287            beta_mode = values[extra_offset]
288            extra_offset += 1
289        else:
290            beta_mode = 0
291        if p_info.effective_radius_type is not None:
292            effective_radius_type = int(values[extra_offset])
293            extra_offset += 1
294        else:
295            effective_radius_type = 0
296
297        # Call the kernels
298        scale, background = values[0], values[1]
299        if beta_mode:
300            F1, F2, volume_avg, effective_radius = self.p_kernel.beta(
301                p_details, p_values, cutoff, magnetic, effective_radius_type)
302            if effective_radius_type > 0:
303                # Plug R_eff from p into S model (initial value and pd value)
304                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
305            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
306            combined_scale = scale*volfrac/volume_avg
307
308            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius)
309            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background
310
311        else:
312            p_result, effective_radius = self.p_kernel.Pq_Reff(
313                p_details, p_values, cutoff, magnetic, effective_radius_type)
314            if effective_radius_type > 0:
315                # Plug R_eff from p into S model (initial value and pd value)
316                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
317            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
318            # remember the parts for plotting later
319            self.results = lambda: _intermediates(p_result, s_result, effective_radius)
320            final_result = scale*(p_result*s_result) + background
321
322        #call_details.show(values)
323        #print("values", values)
324        #p_details.show(p_values)
325        #print("=>", p_result)
326        #s_details.show(s_values)
327        #print("=>", s_result)
328        #import pylab as plt
329        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-')
330        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
331        #plt.figure()
332        return final_result
333
334    def release(self):
335        # type: () -> None
336        self.p_kernel.release()
337        self.s_kernel.release()
338
339
340def calc_er_vr(model_info, call_details, values):
341    # type: (ModelInfo, ParameterSet) -> Tuple[float, float]
342
343    if model_info.ER is None and model_info.VR is None:
344        return 1.0, 1.0
345
346    nvalues = model_info.parameters.nvalues
347    value = values[nvalues:nvalues + call_details.num_weights]
348    weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
349    npars = model_info.parameters.npars
350    # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the
351    # dispersion mesh code doesn't actually care about the nominal parameter
352    # value, p, so set it to None.
353    pairs = [(None, value[offset:offset+length], weight[offset:offset+length])
354             for p, offset, length
355             in zip(model_info.parameters.call_parameters[2:2+npars],
356                    call_details.offset,
357                    call_details.length)
358             if p.type == 'volume']
359    value, weight = dispersion_mesh(model_info, pairs)
360
361    if model_info.ER is not None:
362        individual_radii = model_info.ER(*value)
363        radius_effective = np.sum(weight*individual_radii) / np.sum(weight)
364    else:
365        radius_effective = 1.0
366
367    if model_info.VR is not None:
368        whole, part = model_info.VR(*value)
369        volume_ratio = np.sum(weight*part)/np.sum(weight*whole)
370    else:
371        volume_ratio = 1.0
372
373    return radius_effective, volume_ratio
Note: See TracBrowser for help on using the repository browser.