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