source: sasmodels/sasmodels/kernelpy.py @ 2586ab72

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2586ab72 was 9ae0d2e, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Merge branch 'master' into beta_approx

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