source: sasmodels/sasmodels/details.py @ 9e771a3

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

sort out weights (seems to be correct this time)

  • Property mode set to 100644
File size: 13.1 KB
Line 
1"""
2Kernel Call Details
3===================
4
5When calling sas computational kernels with polydispersity there are a
6number of details that need to be sent to the caller.  This includes the
7list of polydisperse parameters, the number of points in the polydispersity
8weight distribution, and which parameter is the "theta" parameter for
9polar coordinate integration.  The :class:`CallDetails` object maintains
10this data.  Use :func:`build_details` to build a *details* object which
11can be passed to one of the computational kernels.
12"""
13
14from __future__ import print_function
15
16import numpy as np  # type: ignore
17from numpy import pi, cos, sin, radians
18
19try:
20    np.meshgrid([])
21    meshgrid = np.meshgrid
22except Exception:
23    # CRUFT: np.meshgrid requires multiple vectors
24    def meshgrid(*args):
25        if len(args) > 1:
26            return np.meshgrid(*args)
27        else:
28            return [np.asarray(v) for v in args]
29
30try:
31    from typing import List, Tuple, Sequence
32except ImportError:
33    pass
34else:
35    from .modelinfo import ModelInfo
36    from .modelinfo import ParameterTable
37
38
39class CallDetails(object):
40    """
41    Manage the polydispersity information for the kernel call.
42
43    Conceptually, a polydispersity calculation is an integral over a mesh
44    in n-D space where n is the number of polydisperse parameters.  In order
45    to keep the program responsive, and not crash the GPU, only a portion
46    of the mesh is computed at a time.  Meshes with a large number of points
47    will therefore require many calls to the polydispersity loop.  Restarting
48    a nested loop in the middle requires that the indices of the individual
49    mesh dimensions can be computed for the current loop location.  This
50    is handled by the *pd_stride* vector, with n//stride giving the loop
51    index and n%stride giving the position in the sub loops.
52
53    One of the parameters may be the latitude.  When integrating in polar
54    coordinates, the total circumference decreases as latitude varies from
55    pi r^2 at the equator to 0 at the pole, and the weight associated
56    with a range of latitude values needs to be scaled by this circumference.
57    This scale factor needs to be updated each time the theta value
58    changes.  *theta_par* indicates which of the values in the parameter
59    vector is the latitude parameter, or -1 if there is no latitude
60    parameter in the model.  In practice, the normalization term cancels
61    if the latitude is not a polydisperse parameter.
62    """
63    parts = None  # type: List["CallDetails"]
64    def __init__(self, model_info):
65        # type: (ModelInfo) -> None
66        parameters = model_info.parameters
67        max_pd = parameters.max_pd
68
69        # Structure of the call details buffer:
70        #   pd_par[max_pd]     pd params in order of length
71        #   pd_length[max_pd]  length of each pd param
72        #   pd_offset[max_pd]  offset of pd values in parameter array
73        #   pd_stride[max_pd]  index of pd value in loop = n//stride[k]
74        #   num_eval           total length of pd loop
75        #   num_weights        total length of the weight vector
76        #   num_active         number of pd params
77        #   theta_par          parameter number for theta parameter
78        self.buffer = np.empty(4*max_pd + 4, 'i4')
79
80        # generate views on different parts of the array
81        self._pd_par = self.buffer[0 * max_pd:1 * max_pd]
82        self._pd_length = self.buffer[1 * max_pd:2 * max_pd]
83        self._pd_offset = self.buffer[2 * max_pd:3 * max_pd]
84        self._pd_stride = self.buffer[3 * max_pd:4 * max_pd]
85
86        # theta_par is fixed
87        self.theta_par = parameters.theta_offset
88
89        # offset and length are for all parameters, not just pd parameters
90        # They are not sent to the kernel function, though they could be.
91        # They are used by the composite models (sum and product) to
92        # figure out offsets into the combined value list.
93        self.offset = None  # type: np.ndarray
94        self.length = None  # type: np.ndarray
95
96        # keep hold of ifno show() so we can break a values vector
97        # into the individual components
98        self.info = model_info
99
100    @property
101    def pd_par(self):
102        """List of polydisperse parameters"""
103        return self._pd_par
104
105    @property
106    def pd_length(self):
107        """Number of weights for each polydisperse parameter"""
108        return self._pd_length
109
110    @property
111    def pd_offset(self):
112        """Offsets for the individual weight vectors in the set of weights"""
113        return self._pd_offset
114
115    @property
116    def pd_stride(self):
117        """Stride in the pd mesh for each pd dimension"""
118        return self._pd_stride
119
120    @property
121    def num_eval(self):
122        """Total size of the pd mesh"""
123        return self.buffer[-4]
124
125    @num_eval.setter
126    def num_eval(self, v):
127        """Total size of the pd mesh"""
128        self.buffer[-4] = v
129
130    @property
131    def num_weights(self):
132        """Total length of all the weight vectors"""
133        return self.buffer[-3]
134
135    @num_weights.setter
136    def num_weights(self, v):
137        """Total length of all the weight vectors"""
138        self.buffer[-3] = v
139
140    @property
141    def num_active(self):
142        """Number of active polydispersity loops"""
143        return self.buffer[-2]
144
145    @num_active.setter
146    def num_active(self, v):
147        """Number of active polydispersity loops"""
148        self.buffer[-2] = v
149
150    @property
151    def theta_par(self):
152        """Location of the theta parameter in the parameter vector"""
153        return self.buffer[-1]
154
155    @theta_par.setter
156    def theta_par(self, v):
157        """Location of the theta parameter in the parameter vector"""
158        self.buffer[-1] = v
159
160    def show(self, values=None):
161        """Print the polydispersity call details to the console"""
162        print("===== %s details ===="%self.info.name)
163        print("num_active:%d  num_eval:%d  num_weights:%d  theta=%d"
164              % (self.num_active, self.num_eval, self.num_weights, self.theta_par))
165        if self.pd_par.size:
166            print("pd_par", self.pd_par)
167            print("pd_length", self.pd_length)
168            print("pd_offset", self.pd_offset)
169            print("pd_stride", self.pd_stride)
170        if values is not None:
171            nvalues = self.info.parameters.nvalues
172            print("scale, background", values[:2])
173            print("val", values[2:nvalues])
174            print("pd", values[nvalues:nvalues+self.num_weights])
175            print("wt", values[nvalues+self.num_weights:nvalues+2*self.num_weights])
176            print("offsets", self.offset)
177
178
179def make_details(model_info, length, offset, num_weights):
180    # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails
181    """
182    Return a :class:`CallDetails` object for a polydisperse calculation
183    of the model defined by *model_info*.  Polydispersity is defined by
184    the *length* of the polydispersity distribution for each parameter
185    and the *offset* of the distribution in the polydispersity array.
186    Monodisperse parameters should use a polydispersity length of one
187    with weight 1.0. *num_weights* is the total length of the polydispersity
188    array.
189    """
190    #pars = model_info.parameters.call_parameters[2:model_info.parameters.npars+2]
191    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(pars)))
192    #print("len:",length)
193    #print("off:",offset)
194
195    # Check that we arn't using too many polydispersity loops
196    num_active = np.sum(length > 1)
197    max_pd = model_info.parameters.max_pd
198    if num_active > max_pd:
199        raise ValueError("Too many polydisperse parameters")
200
201    # Decreasing list of polydpersity lengths
202    # Note: the reversing view, x[::-1], does not require a copy
203    idx = np.argsort(length)[::-1][:max_pd]
204    pd_stride = np.cumprod(np.hstack((1, length[idx])))
205
206    call_details = CallDetails(model_info)
207    call_details.pd_par[:max_pd] = idx
208    call_details.pd_length[:max_pd] = length[idx]
209    call_details.pd_offset[:max_pd] = offset[idx]
210    call_details.pd_stride[:max_pd] = pd_stride[:-1]
211    call_details.num_eval = pd_stride[-1]
212    call_details.num_weights = num_weights
213    call_details.num_active = num_active
214    call_details.length = length
215    call_details.offset = offset
216    #call_details.show()
217    return call_details
218
219
220ZEROS = tuple([0.]*31)
221def make_kernel_args(kernel, mesh):
222    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
223    """
224    Converts (value, dispersity, weight) for each parameter into kernel pars.
225
226    Returns a CallDetails object indicating the polydispersity, a data object
227    containing the different values, and the magnetic flag indicating whether
228    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are
229    converted to rectangular coordinates (mx, my, mz).
230    """
231    npars = kernel.info.parameters.npars
232    nvalues = kernel.info.parameters.nvalues
233    scalars = [value for value, dispersity, weight in mesh]
234    # skipping scale and background when building values and weights
235    values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ())
236    weights = correct_theta_weights(kernel.info.parameters, dispersity, weights)
237    length = np.array([len(w) for w in weights])
238    offset = np.cumsum(np.hstack((0, length)))
239    call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
240    # Pad value array to a 32 value boundaryd
241    data_len = nvalues + 2*sum(len(v) for v in dispersity)
242    extra = (32 - data_len%32)%32
243    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra])
244    data = data.astype(kernel.dtype)
245    is_magnetic = convert_magnetism(kernel.info.parameters, data)
246    #call_details.show()
247    return call_details, data, is_magnetic
248
249def correct_theta_weights(parameters, dispersity, weights):
250    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray]
251    """
252    If there is a theta parameter, update the weights of that parameter so that
253    the cosine weighting required for polar integration is preserved.
254
255    Avoid evaluation strictly at the pole, which would otherwise send the
256    weight to zero.  This is probably not a problem in practice (if dispersity
257    is +/- 90, then you probably should be using a 1-D model of the circular
258    average).
259
260    Note: scale and background parameters are not include in the tuples for
261    dispersity and weights, so index is parameters.theta_offset, not
262    parameters.theta_offset+2
263
264    Returns updated weights vectors
265    """
266    # TODO: explain in a comment why scale and background are missing
267    # Apparently the parameters.theta_offset similarly skips scale and
268    # and background, so the indexing works out, but they are still shipped
269    # to the kernel, so we need to add two there.
270    if parameters.theta_offset >= 0:
271        index = parameters.theta_offset
272        theta = dispersity[index]
273        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,...
274        theta_weight = abs(cos(radians(theta)))
275        weights = tuple(theta_weight*v if k == index else v
276                        for k, v in enumerate(weights))
277    return weights
278
279
280def convert_magnetism(parameters, values):
281    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
282    """
283    Convert magnetism values from polar to rectangular coordinates.
284
285    Returns True if any magnetism is present.
286    """
287    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
288    mag = mag.reshape(-1, 3)
289    scale = mag[:,0]
290    if np.any(scale):
291        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
292        cos_theta = cos(theta)
293        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
294        mag[:, 1] = scale*sin(theta)  # my
295        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
296        return True
297    else:
298        return False
299
300
301def dispersion_mesh(model_info, mesh):
302    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]]
303    """
304    Create a mesh grid of dispersion parameters and weights.
305
306    *mesh* is a list of (value, dispersity, weights), where the values
307    are the individual parameter values, and (dispersity, weights) is
308    the distribution of parameter values.
309
310    Only the volume parameters should be included in this list.  Orientation
311    parameters do not affect the calculation of effective radius or volume
312    ratio.  This is convenient since it avoids the distinction between
313    value and dispersity that is present in orientation parameters but not
314    shape parameters.
315
316    Returns [p1,p2,...],w where pj is a vector of values for parameter j
317    and w is a vector containing the products for weights for each
318    parameter set in the vector.
319    """
320    value, dispersity, weight = zip(*mesh)
321    #weight = [w if len(w)>0 else [1.] for w in weight]
322    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
323    weight = np.prod(weight, axis=0)
324    dispersity = [v.flatten() for v in meshgrid(*dispersity)]
325    lengths = [par.length for par in model_info.parameters.kernel_parameters
326               if par.type == 'volume']
327    if any(n > 1 for n in lengths):
328        pars = []
329        offset = 0
330        for n in lengths:
331            pars.append(np.vstack(dispersity[offset:offset+n])
332                        if n > 1 else dispersity[offset])
333            offset += n
334        dispersity = pars
335    return dispersity, weight
Note: See TracBrowser for help on using the repository browser.