source: sasmodels/sasmodels/kernelpy.py @ e44432d

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

support hollow models in structure factor calculations

  • 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
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        shell = model_info.shell_volume
167        radius = model_info.effective_radius
168        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume
169                        else (lambda: [volume(*volume_args)]*2) if volume
170                        else (lambda: (1.0, 1.0)))
171        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius
172                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume
173                        else (lambda mode: 1.0))
174
175
176
177    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
178        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
179        if magnetic:
180            raise NotImplementedError("Magnetism not implemented for pure python models")
181        #print("Calling python kernel")
182        #call_details.show(values)
183        radius = ((lambda: 0.0) if effective_radius_type == 0
184                  else (lambda: self._radius(effective_radius_type)))
185        self.result = _loops(
186            self._parameter_vector, self._form, self._volume, radius,
187            self.q_input.nq, call_details, values, cutoff)
188
189    def release(self):
190        # type: () -> None
191        """
192        Free resources associated with the kernel.
193        """
194        self.q_input.release()
195        self.q_input = None
196
197def _loops(parameters,    # type: np.ndarray
198           form,          # type: Callable[[], np.ndarray]
199           form_volume,   # type: Callable[[], float]
200           form_radius,   # type: Callable[[], float]
201           nq,            # type: int
202           call_details,  # type: details.CallDetails
203           values,        # type: np.ndarray
204           cutoff,        # type: float
205          ):
206    # type: (...) -> None
207    ################################################################
208    #                                                              #
209    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
210    #   !!                                                    !!   #
211    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
212    #   !!                                                    !!   #
213    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
214    #                                                              #
215    ################################################################
216
217    # WARNING: Trickery ahead
218    # The parameters[] vector is embedded in the closures for form(),
219    # form_volume() and form_radius().  We set the initial vector from
220    # the values for the model parameters. As we loop through the polydispesity
221    # mesh, we update the components with the polydispersity values before
222    # calling the respective functions.
223    n_pars = len(parameters)
224    parameters[:] = values[2:n_pars+2]
225
226    if call_details.num_active == 0:
227        total = form()
228        weight_norm = 1.0
229        weighted_shell, weighted_form = form_volume()
230        weighted_radius = form_radius()
231
232    else:
233        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
234        pd_weight = values[2+n_pars + call_details.num_weights:]
235
236        weight_norm = 0.0
237        weighted_form = 0.0
238        weighted_shell = 0.0
239        weighted_radius = 0.0
240        partial_weight = np.NaN
241        weight = np.NaN
242
243        p0_par = call_details.pd_par[0]
244        p0_length = call_details.pd_length[0]
245        p0_index = p0_length
246        p0_offset = call_details.pd_offset[0]
247
248        pd_par = call_details.pd_par[:call_details.num_active]
249        pd_offset = call_details.pd_offset[:call_details.num_active]
250        pd_stride = call_details.pd_stride[:call_details.num_active]
251        pd_length = call_details.pd_length[:call_details.num_active]
252
253        total = np.zeros(nq, 'd')
254        for loop_index in range(call_details.num_eval):
255            # update polydispersity parameter values
256            if p0_index == p0_length:
257                pd_index = (loop_index//pd_stride)%pd_length
258                parameters[pd_par] = pd_value[pd_offset+pd_index]
259                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
260                p0_index = loop_index%p0_length
261
262            weight = partial_weight * pd_weight[p0_offset + p0_index]
263            parameters[p0_par] = pd_value[p0_offset + p0_index]
264            p0_index += 1
265            if weight > cutoff:
266                # Call the scattering function
267                # Assume that NaNs are only generated if the parameters are bad;
268                # exclude all q for that NaN.  Even better would be to have an
269                # INVALID expression like the C models, but that is expensive.
270                Iq = np.asarray(form(), 'd')
271                if np.isnan(Iq).any():
272                    continue
273
274                # update value and norm
275                total += weight * Iq
276                weight_norm += weight
277                shell, form = form_volume()
278                weighted_form += weight * form
279                weighted_shell += weight * shell
280                weighted_radius += weight * form_radius()
281
282    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius))
283    return result
284
285
286def _create_default_functions(model_info):
287    """
288    Autogenerate missing functions, such as Iqxy from Iq.
289
290    This only works for Iqxy when Iq is written in python. :func:`make_source`
291    performs a similar role for Iq written in C.  This also vectorizes
292    any functions that are not already marked as vectorized.
293    """
294    # Note: must call create_vector_Iq before create_vector_Iqxy
295    _create_vector_Iq(model_info)
296    _create_vector_Iqxy(model_info)
297
298
299def _create_vector_Iq(model_info):
300    """
301    Define Iq as a vector function if it exists.
302    """
303    Iq = model_info.Iq
304    if callable(Iq) and not getattr(Iq, 'vectorized', False):
305        #print("vectorizing Iq")
306        def vector_Iq(q, *args):
307            """
308            Vectorized 1D kernel.
309            """
310            return np.array([Iq(qi, *args) for qi in q])
311        vector_Iq.vectorized = True
312        model_info.Iq = vector_Iq
313
314
315def _create_vector_Iqxy(model_info):
316    """
317    Define Iqxy as a vector function if it exists, or default it from Iq().
318    """
319    Iqxy = getattr(model_info, 'Iqxy', None)
320    if callable(Iqxy):
321        if not getattr(Iqxy, 'vectorized', False):
322            #print("vectorizing Iqxy")
323            def vector_Iqxy(qx, qy, *args):
324                """
325                Vectorized 2D kernel.
326                """
327                return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)])
328            vector_Iqxy.vectorized = True
329            model_info.Iqxy = vector_Iqxy
330    else:
331        #print("defaulting Iqxy")
332        # Iq is vectorized because create_vector_Iq was already called.
333        Iq = model_info.Iq
334        def default_Iqxy(qx, qy, *args):
335            """
336            Default 2D kernel.
337            """
338            return Iq(np.sqrt(qx**2 + qy**2), *args)
339        default_Iqxy.vectorized = True
340        model_info.Iqxy = default_Iqxy
Note: See TracBrowser for help on using the repository browser.