[dbb0048] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | |
---|
| 4 | import numpy as np |
---|
| 5 | import math |
---|
| 6 | import pyopencl as cl |
---|
| 7 | from weights import GaussianDispersion |
---|
| 8 | from sasmodel import card |
---|
| 9 | |
---|
| 10 | def set_precision(src, qx, qy, dtype): |
---|
| 11 | qx = np.ascontiguousarray(qx, dtype=dtype) |
---|
| 12 | qy = np.ascontiguousarray(qy, dtype=dtype) |
---|
| 13 | if np.dtype(dtype) == np.dtype('float32'): |
---|
| 14 | header = """\ |
---|
| 15 | #define real float |
---|
| 16 | """ |
---|
| 17 | else: |
---|
| 18 | header = """\ |
---|
| 19 | #pragma OPENCL EXTENSION cl_khr_fp64: enable |
---|
| 20 | #define real double |
---|
| 21 | """ |
---|
| 22 | return header+src, qx, qy |
---|
| 23 | |
---|
| 24 | class GpuCoreShellCylinder(object): |
---|
| 25 | PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':-1e-6, 'solvent_sld':0, |
---|
| 26 | 'background':0, 'axis_theta':0, 'axis_phi':0} |
---|
| 27 | PD_PARS = ['radius', 'length', 'thickness', 'axis_phi', 'axis_theta'] |
---|
| 28 | |
---|
| 29 | def __init__(self, qx, qy, dtype='float32'): |
---|
| 30 | #create context, queue, and build program |
---|
| 31 | ctx,_queue = card() |
---|
| 32 | src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-CoreShellCylinder.cpp').read(), qx, qy, dtype=dtype) |
---|
| 33 | self.prg = cl.Program(ctx, src).build() |
---|
| 34 | self.qx, self.qy = qx, qy |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | #buffers |
---|
| 38 | mf = cl.mem_flags |
---|
| 39 | self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) |
---|
| 40 | self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) |
---|
| 41 | self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes) |
---|
| 42 | self.res = np.empty_like(qx) |
---|
| 43 | |
---|
| 44 | def eval(self, pars): |
---|
| 45 | |
---|
| 46 | _ctx,queue = card() |
---|
| 47 | radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) |
---|
| 48 | for base in GpuCoreShellCylinder.PD_PARS] |
---|
| 49 | |
---|
| 50 | radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True) |
---|
| 51 | length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) |
---|
| 52 | thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 10000, True) |
---|
| 53 | axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False) |
---|
| 54 | axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False) |
---|
| 55 | |
---|
| 56 | sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0 |
---|
| 57 | |
---|
| 58 | print radius.value |
---|
| 59 | print thickness.weight |
---|
| 60 | print axis_phi.weight |
---|
| 61 | print axis_theta.weight |
---|
| 62 | print length.value |
---|
| 63 | |
---|
| 64 | size = len(axis_theta.weight) |
---|
| 65 | |
---|
| 66 | real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 |
---|
| 67 | for r in xrange(len(radius.weight)): |
---|
| 68 | for l in xrange(len(length.weight)): |
---|
| 69 | for at in xrange(len(axis_theta.weight)): |
---|
| 70 | for p in xrange(len(axis_phi.weight)): |
---|
| 71 | for th in xrange(len(thickness.weight)): |
---|
| 72 | |
---|
| 73 | self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, |
---|
| 74 | real(axis_theta.value[at]), real(axis_phi.value[p]), real(thickness.value[th]), |
---|
| 75 | real(length.value[l]), real(radius.value[r]), real(pars['scale']), |
---|
| 76 | real(radius.weight[r]), real(length.weight[l]), real(thickness.weight[th]), |
---|
| 77 | real(axis_theta.weight[at]), real(axis_phi.weight[p]), real(pars['core_sld']), |
---|
| 78 | real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size), |
---|
| 79 | np.uint32(self.qx.size)) |
---|
| 80 | cl.enqueue_copy(queue, self.res, self.res_b) |
---|
| 81 | |
---|
| 82 | sum += self.res |
---|
| 83 | vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\ |
---|
| 84 | *(length.value[l]+2.0*thickness.value[th]) |
---|
| 85 | norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th] |
---|
| 86 | norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\ |
---|
| 87 | *axis_phi.weight[p] |
---|
| 88 | |
---|
| 89 | #if size>1: |
---|
| 90 | # norm /= math.asin(1.0) |
---|
| 91 | if vol != 0.0 and norm_vol != 0.0: |
---|
| 92 | sum *= norm_vol/vol |
---|
| 93 | |
---|
| 94 | return sum/norm + pars['background'] |
---|