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

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

re-enable structure factor product

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