source: sasmodels/sasmodels/details.py @ b297ba9

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b297ba9 was b297ba9, checked in by Paul Kienzle <pkienzle@…>, 7 months 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    value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ())
241    #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight)
242    length = np.array([len(w) for w in weight])
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 boundary
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 + weight + ZEROS[:extra])
249    data = data.astype(kernel.dtype)
250    is_magnetic = convert_magnetism(kernel.info.parameters, data)
251    #call_details.show()
252    #print("data", data)
253    return call_details, data, is_magnetic
254
255def correct_theta_weights(parameters, # type: ParameterTable
256                          dispersity, # type: Sequence[np.ndarray]
257                          weights     # type: Sequence[np.ndarray]
258                         ):
259    # type: (...) -> Sequence[np.ndarray]
260    """
261    **Deprecated** Theta weights will be computed in the kernel wrapper if
262    they are needed.
263
264    If there is a theta parameter, update the weights of that parameter so that
265    the cosine weighting required for polar integration is preserved.
266
267    Avoid evaluation strictly at the pole, which would otherwise send the
268    weight to zero.  This is probably not a problem in practice (if dispersity
269    is +/- 90, then you probably should be using a 1-D model of the circular
270    average).
271
272    Note: scale and background parameters are not include in the tuples for
273    dispersity and weights, so index is parameters.theta_offset, not
274    parameters.theta_offset+2
275
276    Returns updated weights vectors
277    """
278    # Apparently the parameters.theta_offset similarly skips scale and
279    # and background, so the indexing works out, but they are still shipped
280    # to the kernel, so we need to add two there.
281    if parameters.theta_offset >= 0:
282        index = parameters.theta_offset
283        theta = dispersity[index]
284        theta_weight = abs(cos(radians(theta)))
285        weights = tuple(theta_weight*w if k == index else w
286                        for k, w in enumerate(weights))
287    return weights
288
289
290def convert_magnetism(parameters, values):
291    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
292    """
293    Convert magnetism values from polar to rectangular coordinates.
294
295    Returns True if any magnetism is present.
296    """
297    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
298    mag = mag.reshape(-1, 3)
299    if np.any(mag[:, 0] != 0.0):
300        M0 = mag[:, 0].copy()
301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
302        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx
303        mag[:, 1] = +M0*sin(theta) # my
304        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz
305        return True
306    else:
307        return False
308
309
310def dispersion_mesh(model_info, mesh):
311    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]]
312    """
313    Create a mesh grid of dispersion parameters and weights.
314
315    *mesh* is a list of (value, dispersity, weights), where the values
316    are the individual parameter values, and (dispersity, weights) is
317    the distribution of parameter values.
318
319    Only the volume parameters should be included in this list.  Orientation
320    parameters do not affect the calculation of effective radius or volume
321    ratio.  This is convenient since it avoids the distinction between
322    value and dispersity that is present in orientation parameters but not
323    shape parameters.
324
325    Returns [p1,p2,...],w where pj is a vector of values for parameter j
326    and w is a vector containing the products for weights for each
327    parameter set in the vector.
328    """
329    _, dispersity, weight = zip(*mesh)
330    #weight = [w if len(w)>0 else [1.] for w in weight]
331    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
332    weight = np.prod(weight, axis=0)
333    dispersity = [v.flatten() for v in meshgrid(*dispersity)]
334    lengths = [par.length for par in model_info.parameters.kernel_parameters
335               if par.type == 'volume']
336    if any(n > 1 for n in lengths):
337        pars = []
338        offset = 0
339        for n in lengths:
340            pars.append(np.vstack(dispersity[offset:offset+n])
341                        if n > 1 else dispersity[offset])
342            offset += n
343        dispersity = pars
344    return dispersity, weight
Note: See TracBrowser for help on using the repository browser.