source: sasmodels/sasmodels/details.py @ 885753a

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

correct magnetic calculations; fractional polarization no longer matches 3.1.2

  • Property mode set to 100644
File size: 13.4 KB
RevLine 
[baa79c2]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
[a738209]14from __future__ import print_function
15
[7ae2b7f]16import numpy as np  # type: ignore
[2d81cfe]17from numpy import cos, sin, radians
[9eb3632]18
19try:
20    np.meshgrid([])
21    meshgrid = np.meshgrid
[725ee36]22except Exception:
[9eb3632]23    # CRUFT: np.meshgrid requires multiple vectors
24    def meshgrid(*args):
[110f69c]25        """See docs from a recent version of numpy"""
[9eb3632]26        if len(args) > 1:
27            return np.meshgrid(*args)
28        else:
29            return [np.asarray(v) for v in args]
[7ae2b7f]30
[2d81cfe]31# pylint: disable=unused-import
[7ae2b7f]32try:
[d4c33d6]33    from typing import List, Tuple, Sequence
[7ae2b7f]34except ImportError:
35    pass
[8d62008]36else:
37    from .modelinfo import ModelInfo
[d4c33d6]38    from .modelinfo import ParameterTable
[2d81cfe]39# pylint: enable=unused-import
[7ae2b7f]40
[d2fc9a4]41
42class CallDetails(object):
[baa79c2]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
[d4c33d6]59    with a range of latitude values needs to be scaled by this circumference.
[baa79c2]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    """
[7ae2b7f]66    parts = None  # type: List["CallDetails"]
[d2fc9a4]67    def __init__(self, model_info):
[8d62008]68        # type: (ModelInfo) -> None
[d2fc9a4]69        parameters = model_info.parameters
70        max_pd = parameters.max_pd
[a738209]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]
[bde38b5]77        #   num_eval           total length of pd loop
78        #   num_weights        total length of the weight vector
[a738209]79        #   num_active         number of pd params
80        #   theta_par          parameter number for theta parameter
[bde38b5]81        self.buffer = np.empty(4*max_pd + 4, 'i4')
[d2fc9a4]82
83        # generate views on different parts of the array
[baa79c2]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]
[d2fc9a4]88
89        # theta_par is fixed
[a738209]90        self.theta_par = parameters.theta_offset
[d2fc9a4]91
[bde38b5]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
[d2fc9a4]103    @property
[baa79c2]104    def pd_par(self):
105        """List of polydisperse parameters"""
106        return self._pd_par
[d2fc9a4]107
108    @property
[baa79c2]109    def pd_length(self):
110        """Number of weights for each polydisperse parameter"""
111        return self._pd_length
[d2fc9a4]112
113    @property
[baa79c2]114    def pd_offset(self):
115        """Offsets for the individual weight vectors in the set of weights"""
116        return self._pd_offset
[d2fc9a4]117
118    @property
[baa79c2]119    def pd_stride(self):
120        """Stride in the pd mesh for each pd dimension"""
121        return self._pd_stride
[d2fc9a4]122
123    @property
[bde38b5]124    def num_eval(self):
[baa79c2]125        """Total size of the pd mesh"""
126        return self.buffer[-4]
127
[bde38b5]128    @num_eval.setter
129    def num_eval(self, v):
[40a87fa]130        """Total size of the pd mesh"""
[baa79c2]131        self.buffer[-4] = v
[d2fc9a4]132
133    @property
[bde38b5]134    def num_weights(self):
[baa79c2]135        """Total length of all the weight vectors"""
136        return self.buffer[-3]
137
[bde38b5]138    @num_weights.setter
139    def num_weights(self, v):
[40a87fa]140        """Total length of all the weight vectors"""
[baa79c2]141        self.buffer[-3] = v
[d2fc9a4]142
143    @property
[baa79c2]144    def num_active(self):
145        """Number of active polydispersity loops"""
146        return self.buffer[-2]
147
[d2fc9a4]148    @num_active.setter
[baa79c2]149    def num_active(self, v):
[40a87fa]150        """Number of active polydispersity loops"""
[baa79c2]151        self.buffer[-2] = v
[d2fc9a4]152
153    @property
[baa79c2]154    def theta_par(self):
155        """Location of the theta parameter in the parameter vector"""
156        return self.buffer[-1]
157
[a738209]158    @theta_par.setter
[baa79c2]159    def theta_par(self, v):
[40a87fa]160        """Location of the theta parameter in the parameter vector"""
[baa79c2]161        self.buffer[-1] = v
[d2fc9a4]162
[bde38b5]163    def show(self, values=None):
[baa79c2]164        """Print the polydispersity call details to the console"""
[bde38b5]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):
[6dc78e4]183    # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails
[baa79c2]184    """
185    Return a :class:`CallDetails` object for a polydisperse calculation
[bde38b5]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.
[baa79c2]192    """
[bde38b5]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
[6aee3ab]198    # Check that we aren't using too many polydispersity loops
[bde38b5]199    num_active = np.sum(length > 1)
[9eb3632]200    max_pd = model_info.parameters.max_pd
201    if num_active > max_pd:
[d2fc9a4]202        raise ValueError("Too many polydisperse parameters")
203
[bde38b5]204    # Decreasing list of polydpersity lengths
[a738209]205    # Note: the reversing view, x[::-1], does not require a copy
[bde38b5]206    idx = np.argsort(length)[::-1][:max_pd]
207    pd_stride = np.cumprod(np.hstack((1, length[idx])))
[d2fc9a4]208
209    call_details = CallDetails(model_info)
[9eb3632]210    call_details.pd_par[:max_pd] = idx
[bde38b5]211    call_details.pd_length[:max_pd] = length[idx]
212    call_details.pd_offset[:max_pd] = offset[idx]
[9eb3632]213    call_details.pd_stride[:max_pd] = pd_stride[:-1]
[bde38b5]214    call_details.num_eval = pd_stride[-1]
215    call_details.num_weights = num_weights
[d2fc9a4]216    call_details.num_active = num_active
[bde38b5]217    call_details.length = length
218    call_details.offset = offset
[d2fc9a4]219    #call_details.show()
220    return call_details
[9eb3632]221
222
[bde38b5]223ZEROS = tuple([0.]*31)
[2d81cfe]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]
[9eb3632]228    """
[8698a0d]229    Converts (value, dispersity, weight) for each parameter into kernel pars.
[9eb3632]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    """
[bde38b5]236    npars = kernel.info.parameters.npars
237    nvalues = kernel.info.parameters.nvalues
[110f69c]238    scalars = [value for value, _dispersity, _weight in mesh]
[77bfb5f]239    # skipping scale and background when building values and weights
[110f69c]240    _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ())
[767dca8]241    #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights)
[bde38b5]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])
[885753a]245    # Pad value array to a 32 value boundary
[8698a0d]246    data_len = nvalues + 2*sum(len(v) for v in dispersity)
[bde38b5]247    extra = (32 - data_len%32)%32
[8698a0d]248    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra])
[bde38b5]249    data = data.astype(kernel.dtype)
[9eb3632]250    is_magnetic = convert_magnetism(kernel.info.parameters, data)
251    #call_details.show()
[885753a]252    #print("data", data)
[9eb3632]253    return call_details, data, is_magnetic
254
[2d81cfe]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]
[d4c33d6]260    """
[108e70e]261    **Deprecated** Theta weights will be computed in the kernel wrapper if
262    they are needed.
263
[d4c33d6]264    If there is a theta parameter, update the weights of that parameter so that
[9e771a3]265    the cosine weighting required for polar integration is preserved.
[77bfb5f]266
[9e771a3]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
[d4c33d6]277    """
[77bfb5f]278    # Apparently the parameters.theta_offset similarly skips scale and
[9e771a3]279    # and background, so the indexing works out, but they are still shipped
280    # to the kernel, so we need to add two there.
[d4c33d6]281    if parameters.theta_offset >= 0:
[77bfb5f]282        index = parameters.theta_offset
[9e771a3]283        theta = dispersity[index]
284        theta_weight = abs(cos(radians(theta)))
[2bccb5a]285        weights = tuple(theta_weight*w if k == index else w
286                        for k, w in enumerate(weights))
[d4c33d6]287    return weights
288
[bde38b5]289
[9eb3632]290def convert_magnetism(parameters, values):
[d4c33d6]291    # type: (ParameterTable, Sequence[np.ndarray]) -> bool
[9eb3632]292    """
[baa79c2]293    Convert magnetism values from polar to rectangular coordinates.
[9eb3632]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)
[885753a]299    if np.any(mag[:, 0] != 0.0):
300        M0 = mag[:, 0].copy()
[d4c33d6]301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2])
[885753a]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
[9eb3632]305        return True
306    else:
307        return False
[bde38b5]308
309
[8698a0d]310def dispersion_mesh(model_info, mesh):
[bde38b5]311    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]]
312    """
313    Create a mesh grid of dispersion parameters and weights.
314
[8698a0d]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.
[d4c33d6]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
[8698a0d]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.
[d4c33d6]324
[bde38b5]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    """
[ce99754]329    _, dispersity, weight = zip(*mesh)
[6dc78e4]330    #weight = [w if len(w)>0 else [1.] for w in weight]
[bde38b5]331    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
332    weight = np.prod(weight, axis=0)
[8698a0d]333    dispersity = [v.flatten() for v in meshgrid(*dispersity)]
[bde38b5]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:
[8698a0d]340            pars.append(np.vstack(dispersity[offset:offset+n])
341                        if n > 1 else dispersity[offset])
[bde38b5]342            offset += n
[8698a0d]343        dispersity = pars
344    return dispersity, weight
Note: See TracBrowser for help on using the repository browser.