source: sasmodels/sasmodels/kernelpy.py @ aa8c6e0

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since aa8c6e0 was aa8c6e0, checked in by Paul Kienzle <pkienzle@…>, 8 months ago

Merge branch 'master' into beta_approx

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