source: sasmodels/sasmodels/kernelpy.py @ 6592f56

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 6592f56 was 14a15a3, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

treat normalization volume of 0. as 1., which mitigates n_shells=0 problem for spherical sld. Refs #635.

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