source: sasmodels/sasmodels/kernelpy.py @ 62d7601

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

Merge branch 'master' into ticket-776-orientation

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