1 | """ |
---|
2 | Python driver for python kernels |
---|
3 | |
---|
4 | Calls the kernel with a vector of $q$ values for a single parameter set. |
---|
5 | Polydispersity is supported by looping over different parameter sets and |
---|
6 | summing the results. The interface to :class:`PyModel` matches those for |
---|
7 | :class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. |
---|
8 | """ |
---|
9 | from __future__ import division, print_function |
---|
10 | |
---|
11 | import logging |
---|
12 | |
---|
13 | import numpy as np # type: ignore |
---|
14 | |
---|
15 | from .generate import F64 |
---|
16 | from .kernel import KernelModel, Kernel |
---|
17 | |
---|
18 | # pylint: disable=unused-import |
---|
19 | try: |
---|
20 | from typing import Union, Callable |
---|
21 | except ImportError: |
---|
22 | pass |
---|
23 | else: |
---|
24 | from . import details |
---|
25 | DType = Union[None, str, np.dtype] |
---|
26 | # pylint: enable=unused-import |
---|
27 | |
---|
28 | logger = logging.getLogger(__name__) |
---|
29 | |
---|
30 | class PyModel(KernelModel): |
---|
31 | """ |
---|
32 | Wrapper for pure python models. |
---|
33 | """ |
---|
34 | def __init__(self, model_info): |
---|
35 | # Make sure Iq is available and vectorized |
---|
36 | _create_default_functions(model_info) |
---|
37 | self.info = model_info |
---|
38 | self.dtype = np.dtype('d') |
---|
39 | |
---|
40 | def make_kernel(self, q_vectors): |
---|
41 | q_input = PyInput(q_vectors, dtype=F64) |
---|
42 | return PyKernel(self.info, q_input) |
---|
43 | |
---|
44 | def release(self): |
---|
45 | """ |
---|
46 | Free resources associated with the model. |
---|
47 | """ |
---|
48 | pass |
---|
49 | |
---|
50 | class PyInput(object): |
---|
51 | """ |
---|
52 | Make q data available to the gpu. |
---|
53 | |
---|
54 | *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data, |
---|
55 | and *[qx, qy]* for 2-D data. Internally, the vectors will be reallocated |
---|
56 | to get the best performance on OpenCL, which may involve shifting and |
---|
57 | stretching the array to better match the memory architecture. Additional |
---|
58 | points will be evaluated with *q=1e-3*. |
---|
59 | |
---|
60 | *dtype* is the data type for the q vectors. The data type should be |
---|
61 | set to match that of the kernel, which is an attribute of |
---|
62 | :class:`GpuProgram`. Note that not all kernels support double |
---|
63 | precision, so even if the program was created for double precision, |
---|
64 | the *GpuProgram.dtype* may be single precision. |
---|
65 | |
---|
66 | Call :meth:`release` when complete. Even if not called directly, the |
---|
67 | buffer will be released when the data object is freed. |
---|
68 | """ |
---|
69 | def __init__(self, q_vectors, dtype): |
---|
70 | self.nq = q_vectors[0].size |
---|
71 | self.dtype = dtype |
---|
72 | self.is_2d = (len(q_vectors) == 2) |
---|
73 | if self.is_2d: |
---|
74 | self.q = np.empty((self.nq, 2), dtype=dtype) |
---|
75 | self.q[:, 0] = q_vectors[0] |
---|
76 | self.q[:, 1] = q_vectors[1] |
---|
77 | else: |
---|
78 | self.q = np.empty(self.nq, dtype=dtype) |
---|
79 | self.q[:self.nq] = q_vectors[0] |
---|
80 | |
---|
81 | def release(self): |
---|
82 | """ |
---|
83 | Free resources associated with the model inputs. |
---|
84 | """ |
---|
85 | self.q = None |
---|
86 | |
---|
87 | class PyKernel(Kernel): |
---|
88 | """ |
---|
89 | Callable SAS kernel. |
---|
90 | |
---|
91 | *kernel* is the kernel object to call. |
---|
92 | |
---|
93 | *model_info* is the module information |
---|
94 | |
---|
95 | *q_input* is the DllInput q vectors at which the kernel should be |
---|
96 | evaluated. |
---|
97 | |
---|
98 | The resulting call method takes the *pars*, a list of values for |
---|
99 | the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) |
---|
100 | vectors for the polydisperse parameters. *cutoff* determines the |
---|
101 | integration limits: any points with combined weight less than *cutoff* |
---|
102 | will not be calculated. |
---|
103 | |
---|
104 | Call :meth:`release` when done with the kernel instance. |
---|
105 | """ |
---|
106 | def __init__(self, model_info, q_input): |
---|
107 | # type: (callable, ModelInfo, List[np.ndarray]) -> None |
---|
108 | self.dtype = np.dtype('d') |
---|
109 | self.info = model_info |
---|
110 | self.q_input = q_input |
---|
111 | self.res = np.empty(q_input.nq, q_input.dtype) |
---|
112 | self.dim = '2d' if q_input.is_2d else '1d' |
---|
113 | |
---|
114 | partable = model_info.parameters |
---|
115 | #kernel_parameters = (partable.iqxy_parameters if q_input.is_2d |
---|
116 | # else partable.iq_parameters) |
---|
117 | kernel_parameters = partable.iq_parameters |
---|
118 | volume_parameters = partable.form_volume_parameters |
---|
119 | |
---|
120 | # Create an array to hold the parameter values. There will be a |
---|
121 | # single array whose values are updated as the calculator goes |
---|
122 | # through the loop. Arguments to the kernel and volume functions |
---|
123 | # will use views into this vector, relying on the fact that a |
---|
124 | # an array of no dimensions acts like a scalar. |
---|
125 | parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') |
---|
126 | |
---|
127 | # Create views into the array to hold the arguments |
---|
128 | offset = 0 |
---|
129 | kernel_args, volume_args = [], [] |
---|
130 | for p in partable.kernel_parameters: |
---|
131 | if p.length == 1: |
---|
132 | # Scalar values are length 1 vectors with no dimensions. |
---|
133 | v = parameter_vector[offset:offset+1].reshape(()) |
---|
134 | else: |
---|
135 | # Vector values are simple views. |
---|
136 | v = parameter_vector[offset:offset+p.length] |
---|
137 | offset += p.length |
---|
138 | if p in kernel_parameters: |
---|
139 | kernel_args.append(v) |
---|
140 | if p in volume_parameters: |
---|
141 | volume_args.append(v) |
---|
142 | |
---|
143 | # Hold on to the parameter vector so we can use it to call kernel later. |
---|
144 | # This may also be required to preserve the views into the vector. |
---|
145 | self._parameter_vector = parameter_vector |
---|
146 | |
---|
147 | # Generate a closure which calls the kernel with the views into the |
---|
148 | # parameter array. |
---|
149 | if q_input.is_2d: |
---|
150 | form = model_info.Iqxy |
---|
151 | qx, qy = q_input.q[:, 0], q_input.q[:, 1] |
---|
152 | self._form = lambda: form(qx, qy, *kernel_args) |
---|
153 | else: |
---|
154 | form = model_info.Iq |
---|
155 | q = q_input.q |
---|
156 | self._form = lambda: form(q, *kernel_args) |
---|
157 | |
---|
158 | # Generate a closure which calls the form_volume if it exists. |
---|
159 | form_volume = model_info.form_volume |
---|
160 | self._volume = ((lambda: form_volume(*volume_args)) if form_volume else |
---|
161 | (lambda: 1.0)) |
---|
162 | |
---|
163 | def __call__(self, call_details, values, cutoff, magnetic): |
---|
164 | # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray |
---|
165 | if magnetic: |
---|
166 | raise NotImplementedError("Magnetism not implemented for pure python models") |
---|
167 | #print("Calling python kernel") |
---|
168 | #call_details.show(values) |
---|
169 | res = _loops(self._parameter_vector, self._form, self._volume, |
---|
170 | self.q_input.nq, call_details, values, cutoff) |
---|
171 | return res |
---|
172 | |
---|
173 | def release(self): |
---|
174 | # type: () -> None |
---|
175 | """ |
---|
176 | Free resources associated with the kernel. |
---|
177 | """ |
---|
178 | self.q_input.release() |
---|
179 | self.q_input = None |
---|
180 | |
---|
181 | def _loops(parameters, # type: np.ndarray |
---|
182 | form, # type: Callable[[], np.ndarray] |
---|
183 | form_volume, # type: Callable[[], float] |
---|
184 | nq, # type: int |
---|
185 | call_details, # type: details.CallDetails |
---|
186 | values, # type: np.ndarray |
---|
187 | cutoff # type: float |
---|
188 | ): |
---|
189 | # type: (...) -> None |
---|
190 | ################################################################ |
---|
191 | # # |
---|
192 | # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # |
---|
193 | # !! !! # |
---|
194 | # !! KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C !! # |
---|
195 | # !! !! # |
---|
196 | # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # |
---|
197 | # # |
---|
198 | ################################################################ |
---|
199 | n_pars = len(parameters) |
---|
200 | parameters[:] = values[2:n_pars+2] |
---|
201 | if call_details.num_active == 0: |
---|
202 | pd_norm = float(form_volume()) |
---|
203 | scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) |
---|
204 | background = values[1] |
---|
205 | return scale*form() + background |
---|
206 | |
---|
207 | pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] |
---|
208 | pd_weight = values[2+n_pars + call_details.num_weights:] |
---|
209 | |
---|
210 | pd_norm = 0.0 |
---|
211 | partial_weight = np.NaN |
---|
212 | weight = np.NaN |
---|
213 | |
---|
214 | p0_par = call_details.pd_par[0] |
---|
215 | p0_length = call_details.pd_length[0] |
---|
216 | p0_index = p0_length |
---|
217 | p0_offset = call_details.pd_offset[0] |
---|
218 | |
---|
219 | pd_par = call_details.pd_par[:call_details.num_active] |
---|
220 | pd_offset = call_details.pd_offset[:call_details.num_active] |
---|
221 | pd_stride = call_details.pd_stride[:call_details.num_active] |
---|
222 | pd_length = call_details.pd_length[:call_details.num_active] |
---|
223 | |
---|
224 | total = np.zeros(nq, 'd') |
---|
225 | for loop_index in range(call_details.num_eval): |
---|
226 | # update polydispersity parameter values |
---|
227 | if p0_index == p0_length: |
---|
228 | pd_index = (loop_index//pd_stride)%pd_length |
---|
229 | parameters[pd_par] = pd_value[pd_offset+pd_index] |
---|
230 | partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) |
---|
231 | p0_index = loop_index%p0_length |
---|
232 | |
---|
233 | weight = partial_weight * pd_weight[p0_offset + p0_index] |
---|
234 | parameters[p0_par] = pd_value[p0_offset + p0_index] |
---|
235 | p0_index += 1 |
---|
236 | if weight > cutoff: |
---|
237 | # Call the scattering function |
---|
238 | # Assume that NaNs are only generated if the parameters are bad; |
---|
239 | # exclude all q for that NaN. Even better would be to have an |
---|
240 | # INVALID expression like the C models, but that is too expensive. |
---|
241 | Iq = np.asarray(form(), 'd') |
---|
242 | if np.isnan(Iq).any(): |
---|
243 | continue |
---|
244 | |
---|
245 | # update value and norm |
---|
246 | total += weight * Iq |
---|
247 | pd_norm += weight * form_volume() |
---|
248 | |
---|
249 | scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) |
---|
250 | background = values[1] |
---|
251 | return scale*total + background |
---|
252 | |
---|
253 | |
---|
254 | def _create_default_functions(model_info): |
---|
255 | """ |
---|
256 | Autogenerate missing functions, such as Iqxy from Iq. |
---|
257 | |
---|
258 | This only works for Iqxy when Iq is written in python. :func:`make_source` |
---|
259 | performs a similar role for Iq written in C. This also vectorizes |
---|
260 | any functions that are not already marked as vectorized. |
---|
261 | """ |
---|
262 | # Note: must call create_vector_Iq before create_vector_Iqxy |
---|
263 | _create_vector_Iq(model_info) |
---|
264 | _create_vector_Iqxy(model_info) |
---|
265 | |
---|
266 | |
---|
267 | def _create_vector_Iq(model_info): |
---|
268 | """ |
---|
269 | Define Iq as a vector function if it exists. |
---|
270 | """ |
---|
271 | Iq = model_info.Iq |
---|
272 | if callable(Iq) and not getattr(Iq, 'vectorized', False): |
---|
273 | #print("vectorizing Iq") |
---|
274 | def vector_Iq(q, *args): |
---|
275 | """ |
---|
276 | Vectorized 1D kernel. |
---|
277 | """ |
---|
278 | return np.array([Iq(qi, *args) for qi in q]) |
---|
279 | vector_Iq.vectorized = True |
---|
280 | model_info.Iq = vector_Iq |
---|
281 | |
---|
282 | |
---|
283 | def _create_vector_Iqxy(model_info): |
---|
284 | """ |
---|
285 | Define Iqxy as a vector function if it exists, or default it from Iq(). |
---|
286 | """ |
---|
287 | Iqxy = getattr(model_info, 'Iqxy', None) |
---|
288 | if callable(Iqxy): |
---|
289 | if not getattr(Iqxy, 'vectorized', False): |
---|
290 | #print("vectorizing Iqxy") |
---|
291 | def vector_Iqxy(qx, qy, *args): |
---|
292 | """ |
---|
293 | Vectorized 2D kernel. |
---|
294 | """ |
---|
295 | return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) |
---|
296 | vector_Iqxy.vectorized = True |
---|
297 | model_info.Iqxy = vector_Iqxy |
---|
298 | else: |
---|
299 | #print("defaulting Iqxy") |
---|
300 | # Iq is vectorized because create_vector_Iq was already called. |
---|
301 | Iq = model_info.Iq |
---|
302 | def default_Iqxy(qx, qy, *args): |
---|
303 | """ |
---|
304 | Default 2D kernel. |
---|
305 | """ |
---|
306 | return Iq(np.sqrt(qx**2 + qy**2), *args) |
---|
307 | default_Iqxy.vectorized = True |
---|
308 | model_info.Iqxy = default_Iqxy |
---|