source: sasmodels/sasmodels/product.py @ 6dc78e4

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 6dc78e4 was 6dc78e4, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

re-enable structure factor product

  • Property mode set to 100644
File size: 10.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 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    assert(s_info.parameters.kernel_parameters[0].id == ER_ID)
48    assert(s_info.parameters.kernel_parameters[1].id == VF_ID)
49    assert(s_info.parameters.magnetism_index == [])
50    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters
51    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
52    p_set = set(p.id for p in p_pars.call_parameters)
53    s_set = set(p.id for p in s_pars.call_parameters)
54
55    if p_set & s_set:
56        # there is some overlap between the parameter names; tag the
57        # overlapping S parameters with name_S.
58        # Skip the first parameter of s, which is effective radius
59        s_list = [(suffix_parameter(par) if par.id in p_set else par)
60                  for par in s_pars.kernel_parameters[1:]]
61    else:
62        # Skip the first parameter of s, which is effective radius
63        s_list = s_pars.kernel_parameters[1:]
64    translate_name = dict((old.id, new.id) for old, new
65                          in zip(s_pars.kernel_parameters[1:], s_list))
66    demo = {}
67    demo.update((k, v) for k, v in p_info.demo.items()
68                if k not in ("background", "scale"))
69    demo.update((translate_name[k], v) for k, v in s_info.demo.items()
70                if k not in ("background", "scale") and not k.startswith(ER_ID))
71    combined_pars = p_pars.kernel_parameters + s_list
72    parameters = ParameterTable(combined_pars)
73    parameters.max_pd = p_pars.max_pd + s_pars.max_pd
74
75    model_info = ModelInfo()
76    model_info.id = '*'.join((p_id, s_id))
77    model_info.name = '*'.join((p_name, s_name))
78    model_info.filename = None
79    model_info.title = 'Product of %s and %s'%(p_name, s_name)
80    model_info.description = model_info.title
81    model_info.docs = model_info.title
82    model_info.category = "custom"
83    model_info.parameters = parameters
84    #model_info.single = p_info.single and s_info.single
85    model_info.structure_factor = False
86    model_info.variant_info = None
87    #model_info.tests = []
88    #model_info.source = []
89    # Iq, Iqxy, form_volume, ER, VR and sesans
90    # Remember the component info blocks so we can build the model
91    model_info.composition = ('product', [p_info, s_info])
92    model_info.demo = {}
93    return model_info
94
95def suffix_parameter(par, suffix):
96    par = copy(par)
97    par.name = par.name + " S"
98    par.id = par.id + "_S"
99    return par
100
101class ProductModel(KernelModel):
102    def __init__(self, model_info, P, S):
103        # type: (ModelInfo, KernelModel, KernelModel) -> None
104        self.info = model_info
105        self.P = P
106        self.S = S
107
108    def make_kernel(self, q_vectors):
109        # type: (List[np.ndarray]) -> Kernel
110        # Note: may be sending the q_vectors to the GPU twice even though they
111        # are only needed once.  It would mess up modularity quite a bit to
112        # handle this optimally, especially since there are many cases where
113        # separate q vectors are needed (e.g., form in python and structure
114        # in opencl; or both in opencl, but one in single precision and the
115        # other in double precision).
116        p_kernel = self.P.make_kernel(q_vectors)
117        s_kernel = self.S.make_kernel(q_vectors)
118        return ProductKernel(self.info, p_kernel, s_kernel)
119
120    def release(self):
121        # type: (None) -> None
122        """
123        Free resources associated with the model.
124        """
125        self.P.release()
126        self.S.release()
127
128
129class ProductKernel(Kernel):
130    def __init__(self, model_info, p_kernel, s_kernel):
131        # type: (ModelInfo, Kernel, Kernel) -> None
132        self.info = model_info
133        self.p_kernel = p_kernel
134        self.s_kernel = s_kernel
135        self.dtype = p_kernel.dtype
136        self.results = []  # type: List[np.ndarray]
137
138    def __call__(self, call_details, values, cutoff, magnetic):
139        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
140        p_info, s_info = self.info.composition[1]
141
142        # if there are magnetic parameters, they will only be on the
143        # form factor P, not the structure factor S.
144        nmagnetic = len(self.info.parameters.magnetism_index)
145        if nmagnetic:
146            spin_index = self.info.parameters.npars + 2
147            magnetism = values[spin_index: spin_index+3+3*nmagnetic]
148        else:
149            magnetism = []
150        nvalues = self.info.parameters.nvalues
151        nweights = call_details.num_weights
152        weights = values[nvalues:nvalues + 2*nweights]
153
154        # Construct the calling parameters for P.
155        p_npars = p_info.parameters.npars
156        p_length = call_details.length[:p_npars]
157        p_offset = call_details.offset[:p_npars]
158        p_details = make_details(p_info, p_length, p_offset, nweights)
159        # Set p scale to the volume fraction in s, which is the first of the
160        # 'S' parameters in the parameter list, or 2+np in 0-origin.
161        volfrac = values[2+p_npars]
162        p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights]
163        spacer = (32 - sum(len(v) for v in p_values)%32)%32
164        p_values.append([0.]*spacer)
165        p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
166
167        # Call ER and VR for P since these are needed for S.
168        p_er, p_vr = calc_er_vr(p_info, p_details, p_values)
169        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac)
170        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr))
171
172        # Construct the calling parameters for S.
173        # The  effective radius is not in the combined parameter list, so
174        # the number of 'S' parameters is one less than expected.  The
175        # computed effective radius needs to be added into the weights
176        # vector, especially since it is a polydisperse parameter in the
177        # stand-alone structure factor models.  We will added it at the
178        # end so the remaining offsets don't need to change.
179        s_npars = s_info.parameters.npars-1
180        s_length = call_details.length[p_npars:p_npars+s_npars]
181        s_offset = call_details.offset[p_npars:p_npars+s_npars]
182        s_length = np.hstack((1, s_length))
183        s_offset = np.hstack((nweights, s_offset))
184        s_details = make_details(s_info, s_length, s_offset, nweights+1)
185        v, w = weights[:nweights], weights[nweights:]
186        s_values = [[1., 0., p_er, s_vr],
187                    # er and vf already included, so start one past the
188                    # volfrac parameter, and go two less than the number
189                    # of S parameters.
190                    values[2+p_npars+2:2+p_npars+s_npars-1],
191                    # no magnetism parameters to include for S
192                    # add er into the (value, weights) pairs
193                    v, [p_er], w, [1.0]]
194        spacer = (32 - sum(len(v) for v in s_values)%32)%32
195        s_values.append([0.]*spacer)
196        s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
197
198        # Call the kernels
199        p_result = self.p_kernel(p_details, p_values, cutoff, magnetic)
200        s_result = self.s_kernel(s_details, s_values, cutoff, False)
201
202        #call_details.show(values)
203        #print("values", values)
204        #p_details.show(p_values)
205        #print("=>", p_result)
206        #print("p val", s_values)
207        #s_details.show(s_values)
208        #print("=>", s_result)
209
210        # remember the parts for plotting later
211        self.results = [p_result, s_result]
212
213        #import pylab as plt
214        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-')
215        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
216        #plt.figure()
217
218        return values[0]*(p_result*s_result) + values[1]
219
220    def release(self):
221        # type: () -> None
222        self.p_kernel.release()
223        self.s_kernel.release()
224
225
226def calc_er_vr(model_info, call_details, values):
227    # type: (ModelInfo, ParameterSet) -> Tuple[float, float]
228
229    if model_info.ER is None and model_info.VR is None:
230        return 1.0, 1.0
231
232    nvalues = model_info.parameters.nvalues
233    value = values[nvalues:nvalues + call_details.num_weights]
234    weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
235    npars = model_info.parameters.npars
236    pairs = [(value[offset:offset+length], weight[offset:offset+length])
237             for p, offset, length
238             in zip(model_info.parameters.call_parameters[2:2+npars],
239                    call_details.offset,
240                    call_details.length)
241             if p.type == 'volume']
242    value, weight = dispersion_mesh(model_info, pairs)
243
244    if model_info.ER is not None:
245        individual_radii = model_info.ER(*value)
246        radius_effective = np.sum(weight*individual_radii) / np.sum(weight)
247    else:
248        radius_effective = 1.0
249
250    if model_info.VR is not None:
251        whole, part = model_info.VR(*value)
252        volume_ratio = np.sum(weight*part)/np.sum(weight*whole)
253    else:
254        volume_ratio = 1.0
255
256    return radius_effective, volume_ratio
Note: See TracBrowser for help on using the repository browser.