source: sasmodels/sasmodels/details.py @ ef07e95

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

lint

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