source: sasmodels/sasmodels/kernelpy.py @ 3199b17

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

PEP 8 changes and improved consistency between OpenCL/CUDA/DLL code

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