source: sasmodels/sasmodels/product.py @ ce99754

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

make sure that nominal values get into the weight vector even when there is no pd so that pd loops are simpler

  • Property mode set to 100644
File size: 12.5 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 copy import copy
16import numpy as np  # type: ignore
17
18from .modelinfo import Parameter, ParameterTable, ModelInfo
19from .kernel import KernelModel, Kernel
20from .details import make_details, dispersion_mesh
21
22try:
23    from typing import Tuple
24except ImportError:
25    pass
26
27# TODO: make estimates available to constraints
28#ESTIMATED_PARAMETERS = [
29#    ["est_radius_effective", "A", 0.0, [0, np.inf], "", "Estimated effective radius"],
30#    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"],
31#]
32
33ER_ID = "radius_effective"
34VF_ID = "volfraction"
35
36# TODO: core_shell_sphere model has suppressed the volume ratio calculation
37# revert it after making VR and ER available at run time as constraints.
38def make_product_info(p_info, s_info):
39    # type: (ModelInfo, ModelInfo) -> ModelInfo
40    """
41    Create info block for product model.
42    """
43    # Make sure effective radius is the first parameter and
44    # make sure volume fraction is the second parameter of the
45    # structure factor calculator.  Structure factors should not
46    # have any magnetic parameters
47    if not len(s_info.parameters.kernel_parameters) >= 2:
48        raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID))
49    if not s_info.parameters.kernel_parameters[0].id == ER_ID:
50        raise TypeError("S needs {} as first parameter".format(ER_ID))
51    if not s_info.parameters.kernel_parameters[1].id == VF_ID:
52        raise TypeError("S needs {} as second parameter".format(VF_ID))
53    if not s_info.parameters.magnetism_index == []:
54        raise TypeError("S should not have SLD parameters")
55    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters
56    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
57
58    # Create list of parameters for the combined model.  Skip the first
59    # parameter of S, which we verified above is effective radius.  If there
60    # are any names in P that overlap with those in S, modify the name in S
61    # to distinguish it.
62    p_set = set(p.id for p in p_pars.kernel_parameters)
63    s_list = [(_tag_parameter(par) if par.id in p_set else par)
64              for par in s_pars.kernel_parameters[1:]]
65    # Check if still a collision after renaming.  This could happen if for
66    # example S has volfrac and P has both volfrac and volfrac_S.
67    if any(p.id in p_set for p in s_list):
68        raise TypeError("name collision: P has P.name and P.name_S while S has S.name")
69
70    translate_name = dict((old.id, new.id) for old, new
71                          in zip(s_pars.kernel_parameters[1:], s_list))
72    combined_pars = p_pars.kernel_parameters + s_list
73    parameters = ParameterTable(combined_pars)
74    parameters.max_pd = p_pars.max_pd + s_pars.max_pd
75    def random():
76        combined_pars = p_info.random()
77        s_names = set(par.id for par in s_pars.kernel_parameters[1:])
78        s = s_info.random()
79        combined_pars.update((translate_name[k], v)
80                    for k, v in s_info.random().items()
81                    if k in s_names)
82        return combined_pars
83
84    model_info = ModelInfo()
85    model_info.id = '@'.join((p_id, s_id))
86    model_info.name = '@'.join((p_name, s_name))
87    model_info.filename = None
88    model_info.title = 'Product of %s and %s'%(p_name, s_name)
89    model_info.description = model_info.title
90    model_info.docs = model_info.title
91    model_info.category = "custom"
92    model_info.parameters = parameters
93    model_info.random = random
94    #model_info.single = p_info.single and s_info.single
95    model_info.structure_factor = False
96    model_info.variant_info = None
97    #model_info.tests = []
98    #model_info.source = []
99    # Iq, Iqxy, form_volume, ER, VR and sesans
100    # Remember the component info blocks so we can build the model
101    model_info.composition = ('product', [p_info, s_info])
102    model_info.control = p_info.control
103    model_info.hidden = p_info.hidden
104    if getattr(p_info, 'profile', None) is not None:
105        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
106        def profile(**kwargs):
107            # extract the profile args
108            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
109            return p_info.profile(**kwargs)
110    else:
111        profile = None
112    model_info.profile = profile
113    model_info.profile_axes = p_info.profile_axes
114
115    # TODO: delegate random to p_info, s_info
116    #model_info.random = lambda: {}
117
118    ## Show the parameter table
119    #from .compare import get_pars, parlist
120    #print("==== %s ====="%model_info.name)
121    #values = get_pars(model_info)
122    #print(parlist(model_info, values, is2d=True))
123    return model_info
124
125def _tag_parameter(par):
126    """
127    Tag the parameter name with _S to indicate that the parameter comes from
128    the structure factor parameter set.  This is only necessary if the
129    form factor model includes a parameter of the same name as a parameter
130    in the structure factor.
131    """
132    par = copy(par)
133    # Protect against a vector parameter in S by appending the vector length
134    # to the renamed parameter.  Note: haven't tested this since no existing
135    # structure factor models contain vector parameters.
136    vector_length = par.name[len(par.id):]
137    par.id = par.id + "_S"
138    par.name = par.id + vector_length
139    return par
140
141class ProductModel(KernelModel):
142    def __init__(self, model_info, P, S):
143        # type: (ModelInfo, KernelModel, KernelModel) -> None
144        self.info = model_info
145        self.P = P
146        self.S = S
147        self.dtype = P.dtype
148
149    def make_kernel(self, q_vectors):
150        # type: (List[np.ndarray]) -> Kernel
151        # Note: may be sending the q_vectors to the GPU twice even though they
152        # are only needed once.  It would mess up modularity quite a bit to
153        # handle this optimally, especially since there are many cases where
154        # separate q vectors are needed (e.g., form in python and structure
155        # in opencl; or both in opencl, but one in single precision and the
156        # other in double precision).
157        p_kernel = self.P.make_kernel(q_vectors)
158        s_kernel = self.S.make_kernel(q_vectors)
159        return ProductKernel(self.info, p_kernel, s_kernel)
160
161    def release(self):
162        # type: (None) -> None
163        """
164        Free resources associated with the model.
165        """
166        self.P.release()
167        self.S.release()
168
169
170class ProductKernel(Kernel):
171    def __init__(self, model_info, p_kernel, s_kernel):
172        # type: (ModelInfo, Kernel, Kernel) -> None
173        self.info = model_info
174        self.p_kernel = p_kernel
175        self.s_kernel = s_kernel
176        self.dtype = p_kernel.dtype
177        self.results = []  # type: List[np.ndarray]
178
179    def __call__(self, call_details, values, cutoff, magnetic):
180        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
181        p_info, s_info = self.info.composition[1]
182
183        # if there are magnetic parameters, they will only be on the
184        # form factor P, not the structure factor S.
185        nmagnetic = len(self.info.parameters.magnetism_index)
186        if nmagnetic:
187            spin_index = self.info.parameters.npars + 2
188            magnetism = values[spin_index: spin_index+3+3*nmagnetic]
189        else:
190            magnetism = []
191        nvalues = self.info.parameters.nvalues
192        nweights = call_details.num_weights
193        weights = values[nvalues:nvalues + 2*nweights]
194
195        # Construct the calling parameters for P.
196        p_npars = p_info.parameters.npars
197        p_length = call_details.length[:p_npars]
198        p_offset = call_details.offset[:p_npars]
199        p_details = make_details(p_info, p_length, p_offset, nweights)
200        # Set p scale to the volume fraction in s, which is the first of the
201        # 'S' parameters in the parameter list, or 2+np in 0-origin.
202        volfrac = values[2+p_npars]
203        p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights]
204        spacer = (32 - sum(len(v) for v in p_values)%32)%32
205        p_values.append([0.]*spacer)
206        p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
207
208        # Call ER and VR for P since these are needed for S.
209        p_er, p_vr = calc_er_vr(p_info, p_details, p_values)
210        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac)
211        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr))
212
213        # Construct the calling parameters for S.
214        # The  effective radius is not in the combined parameter list, so
215        # the number of 'S' parameters is one less than expected.  The
216        # computed effective radius needs to be added into the weights
217        # vector, especially since it is a polydisperse parameter in the
218        # stand-alone structure factor models.  We will added it at the
219        # end so the remaining offsets don't need to change.
220        s_npars = s_info.parameters.npars-1
221        s_length = call_details.length[p_npars:p_npars+s_npars]
222        s_offset = call_details.offset[p_npars:p_npars+s_npars]
223        s_length = np.hstack((1, s_length))
224        s_offset = np.hstack((nweights, s_offset))
225        s_details = make_details(s_info, s_length, s_offset, nweights+1)
226        v, w = weights[:nweights], weights[nweights:]
227        s_values = [
228            # scale=1, background=0, radius_effective=p_er, volfraction=s_vr
229            [1., 0., p_er, s_vr],
230            # structure factor parameters start after scale, background and
231            # all the form factor parameters.  Skip the volfraction parameter
232            # as well, since it is computed elsewhere, and go to the end of the
233            # parameter list.
234            values[2+p_npars+1:2+p_npars+s_npars],
235            # no magnetism parameters to include for S
236            # add er into the (value, weights) pairs
237            v, [p_er], w, [1.0]
238        ]
239        spacer = (32 - sum(len(v) for v in s_values)%32)%32
240        s_values.append([0.]*spacer)
241        s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
242
243        # Call the kernels
244        p_result = self.p_kernel(p_details, p_values, cutoff, magnetic)
245        s_result = self.s_kernel(s_details, s_values, cutoff, False)
246
247        #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars])
248        #call_details.show(values)
249        #print("values", values)
250        #p_details.show(p_values)
251        #print("=>", p_result)
252        #s_details.show(s_values)
253        #print("=>", s_result)
254
255        # remember the parts for plotting later
256        self.results = [p_result, s_result]
257
258        #import pylab as plt
259        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-')
260        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
261        #plt.figure()
262
263        return values[0]*(p_result*s_result) + values[1]
264
265    def release(self):
266        # type: () -> None
267        self.p_kernel.release()
268        self.s_kernel.release()
269
270
271def calc_er_vr(model_info, call_details, values):
272    # type: (ModelInfo, ParameterSet) -> Tuple[float, float]
273
274    if model_info.ER is None and model_info.VR is None:
275        return 1.0, 1.0
276
277    nvalues = model_info.parameters.nvalues
278    value = values[nvalues:nvalues + call_details.num_weights]
279    weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
280    npars = model_info.parameters.npars
281    # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the
282    # dispersion mesh code doesn't actually care about the nominal parameter
283    # value, p, so set it to None.
284    pairs = [(None, value[offset:offset+length], weight[offset:offset+length])
285             for p, offset, length
286             in zip(model_info.parameters.call_parameters[2:2+npars],
287                    call_details.offset,
288                    call_details.length)
289             if p.type == 'volume']
290    value, weight = dispersion_mesh(model_info, pairs)
291
292    if model_info.ER is not None:
293        individual_radii = model_info.ER(*value)
294        radius_effective = np.sum(weight*individual_radii) / np.sum(weight)
295    else:
296        radius_effective = 1.0
297
298    if model_info.VR is not None:
299        whole, part = model_info.VR(*value)
300        volume_ratio = np.sum(weight*part)/np.sum(weight*whole)
301    else:
302        volume_ratio = 1.0
303
304    return radius_effective, volume_ratio
Note: See TracBrowser for help on using the repository browser.