source: sasmodels/sasmodels/kernelpy.py @ 352964b

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 352964b was 352964b, checked in by wojciech, 8 years ago

kernel.py a few checks

  • Property mode set to 100644
File size: 20.8 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
10from numpy import pi, cos
11
12from .generate import F64
13
14class PyModel(object):
15    """
16    Wrapper for pure python models.
17    """
18    def __init__(self, model_info):
19        self.info = model_info
20
21    def __call__(self, q_vectors):
22        q_input = PyInput(q_vectors, dtype=F64)
23        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']
24        return PyKernel(kernel, self.info, q_input)
25
26    def release(self):
27        """
28        Free resources associated with the model.
29        """
30        pass
31
32class PyInput(object):
33    """
34    Make q data available to the gpu.
35
36    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
37    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
38    to get the best performance on OpenCL, which may involve shifting and
39    stretching the array to better match the memory architecture.  Additional
40    points will be evaluated with *q=1e-3*.
41
42    *dtype* is the data type for the q vectors. The data type should be
43    set to match that of the kernel, which is an attribute of
44    :class:`GpuProgram`.  Note that not all kernels support double
45    precision, so even if the program was created for double precision,
46    the *GpuProgram.dtype* may be single precision.
47
48    Call :meth:`release` when complete.  Even if not called directly, the
49    buffer will be released when the data object is freed.
50    """
51    def __init__(self, q_vectors, dtype):
52        self.nq = q_vectors[0].size
53        self.dtype = dtype
54        self.is_2d = (len(q_vectors) == 2)
55        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
56        self.q_pointers = [q.ctypes.data for q in self.q_vectors]
57
58    def release(self):
59        """
60        Free resources associated with the model inputs.
61        """
62        self.q_vectors = []
63
64class PyKernel(object):
65    """
66    Callable SAS kernel.
67
68    *kernel* is the DllKernel object to call.
69
70    *model_info* is the module information
71
72    *q_input* is the DllInput q vectors at which the kernel should be
73    evaluated.
74
75    The resulting call method takes the *pars*, a list of values for
76    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
77    vectors for the polydisperse parameters.  *cutoff* determines the
78    integration limits: any points with combined weight less than *cutoff*
79    will not be calculated.
80
81    Call :meth:`release` when done with the kernel instance.
82    """
83    def __init__(self, kernel, model_info, q_input):
84        self.info = model_info
85        self.q_input = q_input
86        self.res = np.empty(q_input.nq, q_input.dtype)
87        dim = '2d' if q_input.is_2d else '1d'
88        # Loop over q unless user promises that the kernel is vectorized by
89        # taggining it with vectorized=True
90        if not getattr(kernel, 'vectorized', False):
91            if dim == '2d':
92                def vector_kernel(qx, qy, *args):
93                    """
94                    Vectorized 2D kernel.
95                    """
96                    return np.array([kernel(qxi, qyi, *args)
97                                     for qxi, qyi in zip(qx, qy)])
98            else:
99                def vector_kernel(q, *args):
100                    """
101                    Vectorized 1D kernel.
102                    """
103                    return np.array([kernel(qi, *args)
104                                     for qi in q])
105            self.kernel = vector_kernel
106        else:
107            self.kernel = kernel
108        fixed_pars = model_info['partype']['fixed-' + dim]
109        pd_pars = model_info['partype']['pd-' + dim]
110        vol_pars = model_info['partype']['volume']
111
112        # First two fixed pars are scale and background
113        pars = [p.name for p in model_info['parameters'][2:]]
114        offset = len(self.q_input.q_vectors)
115        self.args = self.q_input.q_vectors + [None] * len(pars)
116        self.fixed_index = np.array([pars.index(p) + offset
117                                     for p in fixed_pars[2:]])
118        self.pd_index = np.array([pars.index(p) + offset
119                                  for p in pd_pars])
120        self.vol_index = np.array([pars.index(p) + offset
121                                   for p in vol_pars])
122        try: self.theta_index = pars.index('theta') + offset
123        except ValueError: self.theta_index = -1
124
125        # Caller needs fixed_pars and pd_pars
126        self.fixed_pars = fixed_pars
127        self.pd_pars = pd_pars
128
129    def __call__(self, fixed, pd, cutoff=1e-5):
130        print("fixed",fixed)
131        print("pd", pd)
132        args = self.args[:]  # grab a copy of the args
133        form, form_volume = self.kernel, self.info['form_volume']
134        # First two fixed
135        scale, background = fixed[:2]
136        for index, value in zip(self.fixed_index, fixed[2:]):
137            args[index] = float(value)
138        res = _loops(form, form_volume, cutoff, scale, background, args,
139                     pd, self.pd_index, self.vol_index, self.theta_index)
140
141        return res
142
143    def release(self):
144        """
145        Free resources associated with the kernel.
146        """
147        self.q_input = None
148
149def _loops(form, form_volume, cutoff, scale, background,
150           args, pd, pd_index, vol_index, theta_index):
151    """
152    *form* is the name of the form function, which should be vectorized over
153    q, but otherwise have an interface like the opencl kernels, with the
154    q parameters first followed by the individual parameters in the order
155    given in model.parameters (see :mod:`sasmodels.generate`).
156
157    *form_volume* calculates the volume of the shape.  *vol_index* gives
158    the list of volume parameters
159
160    *cutoff* ignores the corners of the dispersion hypercube
161
162    *scale*, *background* multiplies the resulting form and adds an offset
163
164    *args* is the prepopulated set of arguments to the form function, starting
165    with the q vectors, and including slots for all the parameters.  The
166    values for the parameters will be substituted with values from the
167    dispersion functions.
168
169    *pd* is the list of dispersion parameters
170
171    *pd_index* are the indices of the dispersion parameters in *args*
172
173    *vol_index* are the indices of the volume parameters in *args*
174
175    *theta_index* is the index of the theta parameter for the sasview
176    spherical correction, or -1 if there is no angular dispersion
177    """
178
179    ################################################################
180    #                                                              #
181    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
182    #   !!                                                    !!   #
183    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
184    #   !!                                                    !!   #
185    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
186    #                                                              #
187    ################################################################
188
189    #TODO: Wojtek's notes
190    #TODO: The goal is to restructure polydispersity loop
191    #so it allows fitting arbitrary polydispersity parameters
192    #and it accepts cases like coupled parameters
193    weight = np.empty(len(pd), 'd')
194    if weight.size > 0:
195        # weight vector, to be populated by polydispersity loops
196
197        # identify which pd parameters are volume parameters
198        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
199
200        # Sort parameters in decreasing order of pd length
201        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
202        order = np.argsort(Npd)[::-1]
203        stride = np.cumprod(Npd[order])
204        pd = [pd[index] for index in order]
205        pd_index = pd_index[order]
206        vol_weight_index = vol_weight_index[order]
207
208        fast_value = pd[0][0]
209        fast_weight = pd[0][1]
210    else:
211        stride = np.array([1])
212        vol_weight_index = slice(None, None)
213        # keep lint happy
214        fast_value = [None]
215        fast_weight = [None]
216
217    ret = np.zeros_like(args[0])
218    norm = np.zeros_like(ret)
219    vol = np.zeros_like(ret)
220    vol_norm = np.zeros_like(ret)
221    for k in range(stride[-1]):
222        # update polydispersity parameter values
223        fast_index = k % stride[0]
224        if fast_index == 0:  # bottom loop complete ... check all other loops
225            if weight.size > 0:
226                for i, index, in enumerate(k % stride):
227                    args[pd_index[i]] = pd[i][0][index]
228                    weight[i] = pd[i][1][index]
229        else:
230            args[pd_index[0]] = fast_value[fast_index]
231            weight[0] = fast_weight[fast_index]
232
233        # Computes the weight, and if it is not sufficient then ignore this
234        # parameter set.
235        # Note: could precompute w1*...*wn so we only need to multiply by w0
236        w = np.prod(weight)
237        if w > cutoff:
238            # Note: can precompute spherical correction if theta_index is not
239            # the fast index. Correction factor for spherical integration
240            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
241            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
242                                    if theta_index >= 0 else 1.0)
243            #spherical_correction = 1.0
244
245            # Call the scattering function and adds it to the total.
246            I = form(*args)
247            if np.isnan(I).any(): continue
248            ret += w * I * spherical_correction
249            norm += w
250
251            # Volume normalization.
252            # If there are "volume" polydispersity parameters, then these
253            # will be used to call the form_volume function from the user
254            # supplied kernel, and accumulate a normalized weight.
255            # Note: can precompute volume norm if fast index is not a volume
256            if form_volume:
257                vol_args = [args[index] for index in vol_index]
258                vol_weight = np.prod(weight[vol_weight_index])
259                vol += vol_weight * form_volume(*vol_args)
260                vol_norm += vol_weight
261
262    positive = (vol * vol_norm != 0.0)
263    ret[positive] *= vol_norm[positive] / vol[positive]
264    result = scale * ret / norm + background
265    return result
266=======
267"""
268Python driver for python kernels
269
270Calls the kernel with a vector of $q$ values for a single parameter set.
271Polydispersity is supported by looping over different parameter sets and
272summing the results.  The interface to :class:`PyModel` matches those for
273:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`.
274"""
275import numpy as np
276from numpy import pi, cos
277
278from .generate import F64
279
280class PyModel(object):
281    """
282    Wrapper for pure python models.
283    """
284    def __init__(self, model_info):
285        self.info = model_info
286
287    def __call__(self, q_vectors):
288        q_input = PyInput(q_vectors, dtype=F64)
289        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']
290        return PyKernel(kernel, self.info, q_input)
291
292    def release(self):
293        """
294        Free resources associated with the model.
295        """
296        pass
297
298class PyInput(object):
299    """
300    Make q data available to the gpu.
301
302    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
303    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
304    to get the best performance on OpenCL, which may involve shifting and
305    stretching the array to better match the memory architecture.  Additional
306    points will be evaluated with *q=1e-3*.
307
308    *dtype* is the data type for the q vectors. The data type should be
309    set to match that of the kernel, which is an attribute of
310    :class:`GpuProgram`.  Note that not all kernels support double
311    precision, so even if the program was created for double precision,
312    the *GpuProgram.dtype* may be single precision.
313
314    Call :meth:`release` when complete.  Even if not called directly, the
315    buffer will be released when the data object is freed.
316    """
317    def __init__(self, q_vectors, dtype):
318        self.nq = q_vectors[0].size
319        self.dtype = dtype
320        self.is_2d = (len(q_vectors) == 2)
321        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
322        self.q_pointers = [q.ctypes.data for q in self.q_vectors]
323
324    def release(self):
325        """
326        Free resources associated with the model inputs.
327        """
328        self.q_vectors = []
329
330class PyKernel(object):
331    """
332    Callable SAS kernel.
333
334    *kernel* is the DllKernel object to call.
335
336    *model_info* is the module information
337
338    *q_input* is the DllInput q vectors at which the kernel should be
339    evaluated.
340
341    The resulting call method takes the *pars*, a list of values for
342    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
343    vectors for the polydisperse parameters.  *cutoff* determines the
344    integration limits: any points with combined weight less than *cutoff*
345    will not be calculated.
346
347    Call :meth:`release` when done with the kernel instance.
348    """
349    def __init__(self, kernel, model_info, q_input):
350        self.info = model_info
351        self.q_input = q_input
352        self.res = np.empty(q_input.nq, q_input.dtype)
353        dim = '2d' if q_input.is_2d else '1d'
354        # Loop over q unless user promises that the kernel is vectorized by
355        # taggining it with vectorized=True
356        if not getattr(kernel, 'vectorized', False):
357            if dim == '2d':
358                def vector_kernel(qx, qy, *args):
359                    """
360                    Vectorized 2D kernel.
361                    """
362                    return np.array([kernel(qxi, qyi, *args)
363                                     for qxi, qyi in zip(qx, qy)])
364            else:
365                def vector_kernel(q, *args):
366                    """
367                    Vectorized 1D kernel.
368                    """
369                    return np.array([kernel(qi, *args)
370                                     for qi in q])
371            self.kernel = vector_kernel
372        else:
373            self.kernel = kernel
374        fixed_pars = model_info['partype']['fixed-' + dim]
375        pd_pars = model_info['partype']['pd-' + dim]
376        vol_pars = model_info['partype']['volume']
377
378        # First two fixed pars are scale and background
379        pars = [p.name for p in model_info['parameters'][2:]]
380        offset = len(self.q_input.q_vectors)
381        self.args = self.q_input.q_vectors + [None] * len(pars)
382        self.fixed_index = np.array([pars.index(p) + offset
383                                     for p in fixed_pars[2:]])
384        self.pd_index = np.array([pars.index(p) + offset
385                                  for p in pd_pars])
386        self.vol_index = np.array([pars.index(p) + offset
387                                   for p in vol_pars])
388        try: self.theta_index = pars.index('theta') + offset
389        except ValueError: self.theta_index = -1
390
391        # Caller needs fixed_pars and pd_pars
392        self.fixed_pars = fixed_pars
393        self.pd_pars = pd_pars
394
395    def __call__(self, fixed, pd, cutoff=1e-5):
396        #print("fixed",fixed)
397        #print("pd", pd)
398        args = self.args[:]  # grab a copy of the args
399        form, form_volume = self.kernel, self.info['form_volume']
400        # First two fixed
401        scale, background = fixed[:2]
402        for index, value in zip(self.fixed_index, fixed[2:]):
403            args[index] = float(value)
404        res = _loops(form, form_volume, cutoff, scale, background, args,
405                     pd, self.pd_index, self.vol_index, self.theta_index)
406
407        return res
408
409    def release(self):
410        """
411        Free resources associated with the kernel.
412        """
413        self.q_input = None
414
415def _loops(form, form_volume, cutoff, scale, background,
416           args, pd, pd_index, vol_index, theta_index):
417    """
418    *form* is the name of the form function, which should be vectorized over
419    q, but otherwise have an interface like the opencl kernels, with the
420    q parameters first followed by the individual parameters in the order
421    given in model.parameters (see :mod:`sasmodels.generate`).
422
423    *form_volume* calculates the volume of the shape.  *vol_index* gives
424    the list of volume parameters
425
426    *cutoff* ignores the corners of the dispersion hypercube
427
428    *scale*, *background* multiplies the resulting form and adds an offset
429
430    *args* is the prepopulated set of arguments to the form function, starting
431    with the q vectors, and including slots for all the parameters.  The
432    values for the parameters will be substituted with values from the
433    dispersion functions.
434
435    *pd* is the list of dispersion parameters
436
437    *pd_index* are the indices of the dispersion parameters in *args*
438
439    *vol_index* are the indices of the volume parameters in *args*
440
441    *theta_index* is the index of the theta parameter for the sasview
442    spherical correction, or -1 if there is no angular dispersion
443    """
444
445    ################################################################
446    #                                                              #
447    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
448    #   !!                                                    !!   #
449    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
450    #   !!                                                    !!   #
451    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
452    #                                                              #
453    ################################################################
454
455    weight = np.empty(len(pd), 'd')
456    if weight.size > 0:
457        # weight vector, to be populated by polydispersity loops
458
459        # identify which pd parameters are volume parameters
460        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
461
462        # Sort parameters in decreasing order of pd length
463        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
464        order = np.argsort(Npd)[::-1]
465        stride = np.cumprod(Npd[order])
466        pd = [pd[index] for index in order]
467        pd_index = pd_index[order]
468        vol_weight_index = vol_weight_index[order]
469
470        fast_value = pd[0][0]
471        fast_weight = pd[0][1]
472    else:
473        stride = np.array([1])
474        vol_weight_index = slice(None, None)
475        # keep lint happy
476        fast_value = [None]
477        fast_weight = [None]
478
479    ret = np.zeros_like(args[0])
480    norm = np.zeros_like(ret)
481    vol = np.zeros_like(ret)
482    vol_norm = np.zeros_like(ret)
483    for k in range(stride[-1]):
484        # update polydispersity parameter values
485        fast_index = k % stride[0]
486        if fast_index == 0:  # bottom loop complete ... check all other loops
487            if weight.size > 0:
488                for i, index, in enumerate(k % stride):
489                    args[pd_index[i]] = pd[i][0][index]
490                    weight[i] = pd[i][1][index]
491        else:
492            args[pd_index[0]] = fast_value[fast_index]
493            weight[0] = fast_weight[fast_index]
494
495        # Computes the weight, and if it is not sufficient then ignore this
496        # parameter set.
497        # Note: could precompute w1*...*wn so we only need to multiply by w0
498        w = np.prod(weight)
499        if w > cutoff:
500            # Note: can precompute spherical correction if theta_index is not
501            # the fast index. Correction factor for spherical integration
502            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
503            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
504                                    if theta_index >= 0 else 1.0)
505            #spherical_correction = 1.0
506
507            # Call the scattering function and adds it to the total.
508            I = form(*args)
509            if np.isnan(I).any(): continue
510            ret += w * I * spherical_correction
511            norm += w
512
513            # Volume normalization.
514            # If there are "volume" polydispersity parameters, then these
515            # will be used to call the form_volume function from the user
516            # supplied kernel, and accumulate a normalized weight.
517            # Note: can precompute volume norm if fast index is not a volume
518            if form_volume:
519                vol_args = [args[index] for index in vol_index]
520                vol_weight = np.prod(weight[vol_weight_index])
521                vol += vol_weight * form_volume(*vol_args)
522                vol_norm += vol_weight
523
524    positive = (vol * vol_norm != 0.0)
525    ret[positive] *= vol_norm[positive] / vol[positive]
526    result = scale * ret / norm + background
527    return result
Note: See TracBrowser for help on using the repository browser.