source: sasmodels/sasmodels/details.py @ 40a87fa

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

lint and latex cleanup

  • Property mode set to 100644
File size: 10.5 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 ValueError:
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        #   pd_prod            total length of pd loop
74        #   pd_sum             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.zeros(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    @property
89    def pd_par(self):
90        """List of polydisperse parameters"""
91        return self._pd_par
92
93    @property
94    def pd_length(self):
95        """Number of weights for each polydisperse parameter"""
96        return self._pd_length
97
98    @property
99    def pd_offset(self):
100        """Offsets for the individual weight vectors in the set of weights"""
101        return self._pd_offset
102
103    @property
104    def pd_stride(self):
105        """Stride in the pd mesh for each pd dimension"""
106        return self._pd_stride
107
108    @property
109    def pd_prod(self):
110        """Total size of the pd mesh"""
111        return self.buffer[-4]
112
113    @pd_prod.setter
114    def pd_prod(self, v):
115        """Total size of the pd mesh"""
116        self.buffer[-4] = v
117
118    @property
119    def pd_sum(self):
120        """Total length of all the weight vectors"""
121        return self.buffer[-3]
122
123    @pd_sum.setter
124    def pd_sum(self, v):
125        """Total length of all the weight vectors"""
126        self.buffer[-3] = v
127
128    @property
129    def num_active(self):
130        """Number of active polydispersity loops"""
131        return self.buffer[-2]
132
133    @num_active.setter
134    def num_active(self, v):
135        """Number of active polydispersity loops"""
136        self.buffer[-2] = v
137
138    @property
139    def theta_par(self):
140        """Location of the theta parameter in the parameter vector"""
141        return self.buffer[-1]
142
143    @theta_par.setter
144    def theta_par(self, v):
145        """Location of the theta parameter in the parameter vector"""
146        self.buffer[-1] = v
147
148    def show(self):
149        """Print the polydispersity call details to the console"""
150        print("num_active", self.num_active)
151        print("pd_prod", self.pd_prod)
152        print("pd_sum", self.pd_sum)
153        print("theta par", self.theta_par)
154        print("pd_par", self.pd_par)
155        print("pd_length", self.pd_length)
156        print("pd_offset", self.pd_offset)
157        print("pd_stride", self.pd_stride)
158
159
160def mono_details(model_info):
161    # type: (ModelInfo) -> CallDetails
162    """
163    Return a :class:`CallDetails` object for a monodisperse calculation
164    of the model defined by *model_info*.
165    """
166    call_details = CallDetails(model_info)
167    call_details.pd_prod = 1
168    call_details.pd_sum = model_info.parameters.nvalues
169    call_details.pd_par[:] = np.arange(0, model_info.parameters.max_pd)
170    call_details.pd_length[:] = 1
171    call_details.pd_offset[:] = np.arange(0, model_info.parameters.max_pd)
172    call_details.pd_stride[:] = 1
173    return call_details
174
175
176def poly_details(model_info, weights):
177    """
178    Return a :class:`CallDetails` object for a polydisperse calculation
179    of the model defined by *model_info* for the given set of *weights*.
180    *weights* is a list of vectors, one for each parameter.  Monodisperse
181    parameters should provide a weight vector containing [1.0].
182    """
183    # type: (ModelInfo) -> CallDetails
184    #print("weights",weights)
185    #weights = weights[2:] # Skip scale and background
186
187    # Decreasing list of polydispersity lengths
188    #print([p.id for p in model_info.parameters.call_parameters])
189    pd_length = np.array([len(w)
190                          for w in weights[2:2+model_info.parameters.npars]])
191    num_active = np.sum(pd_length > 1)
192    max_pd = model_info.parameters.max_pd
193    if num_active > max_pd:
194        raise ValueError("Too many polydisperse parameters")
195
196    pd_offset = np.cumsum(np.hstack((0, pd_length)))
197    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters)))
198    #print("len:",pd_length)
199    #print("off:",pd_offset)
200    # Note: the reversing view, x[::-1], does not require a copy
201    idx = np.argsort(pd_length)[::-1][:max_pd]
202    pd_stride = np.cumprod(np.hstack((1, pd_length[idx])))
203
204    call_details = CallDetails(model_info)
205    call_details.pd_par[:max_pd] = idx
206    call_details.pd_length[:max_pd] = pd_length[idx]
207    call_details.pd_offset[:max_pd] = pd_offset[idx]
208    call_details.pd_stride[:max_pd] = pd_stride[:-1]
209    call_details.pd_prod = pd_stride[-1]
210    call_details.pd_sum = sum(len(w) for w in weights)
211    call_details.num_active = num_active
212    #call_details.show()
213    return call_details
214
215
216def dispersion_mesh(model_info, pars):
217    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]]
218    """
219    Create a mesh grid of dispersion parameters and weights.
220
221    Returns [p1,p2,...],w where pj is a vector of values for parameter j
222    and w is a vector containing the products for weights for each
223    parameter set in the vector.
224    """
225    value, weight = zip(*pars)
226    weight = [w if w else [1.] for w in weight]
227    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
228    weight = np.prod(weight, axis=0)
229    value = [v.flatten() for v in meshgrid(*value)]
230    lengths = [par.length for par in model_info.parameters.kernel_parameters
231               if par.type == 'volume']
232    if any(n > 1 for n in lengths):
233        pars = []
234        offset = 0
235        for n in lengths:
236            pars.append(np.vstack(value[offset:offset+n])
237                        if n > 1 else value[offset])
238            offset += n
239        value = pars
240    return value, weight
241
242
243
244def build_details(kernel, pairs):
245    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
246    """
247    Converts (value, weight) pairs into parameters for the kernel call.
248
249    Returns a CallDetails object indicating the polydispersity, a data object
250    containing the different values, and the magnetic flag indicating whether
251    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are
252    converted to rectangular coordinates (mx, my, mz).
253    """
254    values, weights = zip(*pairs)
255    scalars = [v[0] for v in values]
256    if all(len(w) == 1 for w in weights):
257        call_details = mono_details(kernel.info)
258        # Pad value array to a 32 value boundary
259        data_len = 3*len(scalars)
260        extra = ((data_len+31)//32)*32 - data_len
261        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra,
262                        dtype=kernel.dtype)
263    else:
264        call_details = poly_details(kernel.info, weights)
265        # Pad value array to a 32 value boundary
266        data_len = len(scalars) + 2*sum(len(v) for v in values)
267        extra = ((data_len+31)//32)*32 - data_len
268        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra)
269        data = data.astype(kernel.dtype)
270    is_magnetic = convert_magnetism(kernel.info.parameters, data)
271    #call_details.show()
272    return call_details, data, is_magnetic
273
274def convert_magnetism(parameters, values):
275    """
276    Convert magnetism values from polar to rectangular coordinates.
277
278    Returns True if any magnetism is present.
279    """
280    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
281    mag = mag.reshape(-1, 3)
282    scale = mag[:,0]
283    if np.any(scale):
284        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180.
285        cos_theta = cos(theta)
286        mag[:, 0] = scale*cos_theta*cos(phi)  # mx
287        mag[:, 1] = scale*sin(theta)  # my
288        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz
289        return True
290    else:
291        return False
Note: See TracBrowser for help on using the repository browser.