source: sasmodels/sasmodels/kernelpy.py @ 7ae2b7f

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

still more linting; ignore numpy types

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