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