source: sasmodels/sasmodels/details.py @ 4f1f876

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

Intel GPU wants data vectors to follow cache alignment

  • Property mode set to 100644
File size: 7.2 KB
Line 
1from __future__ import print_function
2
3import numpy as np  # type: ignore
4from numpy import pi, cos, sin
5
6try:
7    np.meshgrid([])
8    meshgrid = np.meshgrid
9except ValueError:
10    # CRUFT: np.meshgrid requires multiple vectors
11    def meshgrid(*args):
12        if len(args) > 1:
13            return np.meshgrid(*args)
14        else:
15            return [np.asarray(v) for v in args]
16
17try:
18    from typing import List
19except ImportError:
20    pass
21else:
22    from .modelinfo import ModelInfo
23
24
25class CallDetails(object):
26    parts = None  # type: List["CallDetails"]
27    def __init__(self, model_info):
28        # type: (ModelInfo) -> None
29        parameters = model_info.parameters
30        max_pd = parameters.max_pd
31
32        # Structure of the call details buffer:
33        #   pd_par[max_pd]     pd params in order of length
34        #   pd_length[max_pd]  length of each pd param
35        #   pd_offset[max_pd]  offset of pd values in parameter array
36        #   pd_stride[max_pd]  index of pd value in loop = n//stride[k]
37        #   pd_prod            total length of pd loop
38        #   pd_sum             total length of the weight vector
39        #   num_active         number of pd params
40        #   theta_par          parameter number for theta parameter
41        self.buffer = np.zeros(4*max_pd + 4, 'i4')
42
43        # generate views on different parts of the array
44        self._pd_par     = self.buffer[0 * max_pd:1 * max_pd]
45        self._pd_length  = self.buffer[1 * max_pd:2 * max_pd]
46        self._pd_offset  = self.buffer[2 * max_pd:3 * max_pd]
47        self._pd_stride  = self.buffer[3 * max_pd:4 * max_pd]
48
49        # theta_par is fixed
50        self.theta_par = parameters.theta_offset
51
52    @property
53    def pd_par(self): return self._pd_par
54
55    @property
56    def pd_length(self): return self._pd_length
57
58    @property
59    def pd_offset(self): return self._pd_offset
60
61    @property
62    def pd_stride(self): return self._pd_stride
63
64    @property
65    def pd_prod(self): return self.buffer[-4]
66    @pd_prod.setter
67    def pd_prod(self, v): self.buffer[-4] = v
68
69    @property
70    def pd_sum(self): return self.buffer[-3]
71    @pd_sum.setter
72    def pd_sum(self, v): self.buffer[-3] = v
73
74    @property
75    def num_active(self): return self.buffer[-2]
76    @num_active.setter
77    def num_active(self, v): self.buffer[-2] = v
78
79    @property
80    def theta_par(self): return self.buffer[-1]
81    @theta_par.setter
82    def theta_par(self, v): self.buffer[-1] = v
83
84    def show(self):
85        print("num_active", self.num_active)
86        print("pd_prod", self.pd_prod)
87        print("pd_sum", self.pd_sum)
88        print("theta par", self.theta_par)
89        print("pd_par", self.pd_par)
90        print("pd_length", self.pd_length)
91        print("pd_offset", self.pd_offset)
92        print("pd_stride", self.pd_stride)
93
94
95def mono_details(model_info):
96    call_details = CallDetails(model_info)
97    call_details.pd_prod = 1
98    call_details.pd_sum = model_info.parameters.nvalues
99    call_details.pd_par[:] = np.arange(0, model_info.parameters.max_pd)
100    call_details.pd_length[:] = 1
101    call_details.pd_offset[:] = np.arange(0, model_info.parameters.max_pd)
102    call_details.pd_stride[:] = 1
103    return call_details
104
105
106def poly_details(model_info, weights):
107    #print("weights",weights)
108    #weights = weights[2:] # Skip scale and background
109
110    # Decreasing list of polydispersity lengths
111    #print([p.id for p in model_info.parameters.call_parameters])
112    pd_length = np.array([len(w) for w in weights[2:2+model_info.parameters.npars]])
113    num_active = np.sum(pd_length>1)
114    max_pd = model_info.parameters.max_pd
115    if num_active > max_pd:
116        raise ValueError("Too many polydisperse parameters")
117
118    pd_offset = np.cumsum(np.hstack((0, pd_length)))
119    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters)))
120    #print("len:",pd_length)
121    #print("off:",pd_offset)
122    # Note: the reversing view, x[::-1], does not require a copy
123    idx = np.argsort(pd_length)[::-1][:max_pd]
124    pd_stride = np.cumprod(np.hstack((1, pd_length[idx])))
125
126    call_details = CallDetails(model_info)
127    call_details.pd_par[:max_pd] = idx
128    call_details.pd_length[:max_pd] = pd_length[idx]
129    call_details.pd_offset[:max_pd] = pd_offset[idx]
130    call_details.pd_stride[:max_pd] = pd_stride[:-1]
131    call_details.pd_prod = pd_stride[-1]
132    call_details.pd_sum = sum(len(w) for w in weights)
133    call_details.num_active = num_active
134    #call_details.show()
135    return call_details
136
137
138def dispersion_mesh(model_info, pars):
139    """
140    Create a mesh grid of dispersion parameters and weights.
141
142    Returns [p1,p2,...],w where pj is a vector of values for parameter j
143    and w is a vector containing the products for weights for each
144    parameter set in the vector.
145    """
146    value, weight = zip(*pars)
147    weight = [w if w else [1.] for w in weight]
148    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
149    weight = np.prod(weight, axis=0)
150    value = [v.flatten() for v in meshgrid(*value)]
151    lengths = [par.length for par in model_info.parameters.kernel_parameters
152               if par.type == 'volume']
153    if any(n > 1 for n in lengths):
154        pars = []
155        offset = 0
156        for n in lengths:
157            pars.append(np.vstack(value[offset:offset+n]) if n > 1 else value[offset])
158            offset += n
159        value = pars
160    return value, weight
161
162
163
164def build_details(kernel, pairs):
165    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool]
166    """
167    Converts (value, weight) pairs into parameters for the kernel call.
168
169    Returns a CallDetails object indicating the polydispersity, a data object
170    containing the different values, and the magnetic flag indicating whether
171    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are
172    converted to rectangular coordinates (mx, my, mz).
173    """
174    values, weights = zip(*pairs)
175    scalars = [v[0] for v in values]
176    if all(len(w)==1 for w in weights):
177        call_details = mono_details(kernel.info)
178        # Pad value array to a 32 value boundary
179        data_len = 3*len(scalars)
180        extra = ((data_len+31)//32)*32 - data_len
181        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra, dtype=kernel.dtype)
182    else:
183        call_details = poly_details(kernel.info, weights)
184        # Pad value array to a 32 value boundary
185        data_len = len(scalars) + 2*sum(len(v) for v in values)
186        extra = ((data_len+31)//32)*32 - data_len
187        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra).astype(kernel.dtype)
188    is_magnetic = convert_magnetism(kernel.info.parameters, data)
189    #call_details.show()
190    return call_details, data, is_magnetic
191
192def convert_magnetism(parameters, values):
193    """
194    Convert magnetism in value table from polar to rectangular coordinates.
195
196    Returns True if any magnetism is present.
197    """
198    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues]
199    mag = mag.reshape(-1, 3)
200    M0 = mag[:,0]
201    if np.any(M0):
202        theta, phi = mag[:,1]*pi/180., mag[:,2]*pi/180.
203        cos_theta = cos(theta)
204        mx = M0*cos_theta*cos(phi)
205        my = M0*sin(theta)
206        mz = -M0*cos_theta*sin(phi)
207        mag[:,0], mag[:,1], mag[:,2] = mx, my, mz
208        return True
209    else:
210        return False
Note: See TracBrowser for help on using the repository browser.