source: sasmodels/sasmodels/kernelpy.py @ 96153e4

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 96153e4 was 108e70e, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

use Iqac/Iqabc? for the new orientation interface but Iqxy for the old

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