[5378e40] | 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 |
---|
[8a20be5] | 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 |
---|
[5378e40] | 23 | |
---|
| 24 | class GpuCoreShellCylinder(object): |
---|
[8a20be5] | 25 | PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':-1e-6, 'solvent_sld':0, |
---|
[5378e40] | 26 | 'background':0, 'axis_theta':0, 'axis_phi':0} |
---|
| 27 | PD_PARS = ['radius', 'length', 'thickness', 'axis_phi', 'axis_theta'] |
---|
| 28 | |
---|
[8a20be5] | 29 | def __init__(self, qx, qy, dtype='float32'): |
---|
[5378e40] | 30 | #create context, queue, and build program |
---|
[8a20be5] | 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 | |
---|
[5378e40] | 36 | |
---|
| 37 | #buffers |
---|
| 38 | mf = cl.mem_flags |
---|
[8a20be5] | 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) |
---|
[5378e40] | 43 | |
---|
| 44 | def eval(self, pars): |
---|
| 45 | |
---|
[8a20be5] | 46 | _ctx,queue = card() |
---|
[5378e40] | 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, 1000, True) |
---|
| 51 | length.value, length.weight = length.get_weights(pars['length'], 0, 1000, True) |
---|
| 52 | thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 1000, 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 | size = len(axis_theta.weight) |
---|
| 58 | |
---|
| 59 | for i in xrange(len(radius.weight)): |
---|
| 60 | for j in xrange(len(length.weight)): |
---|
| 61 | for k in xrange(len(axis_theta.weight)): |
---|
| 62 | for l in xrange(len(axis_phi.weight)): |
---|
| 63 | for f in xrange(len(thickness.weight)): |
---|
| 64 | |
---|
[8a20be5] | 65 | self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, |
---|
[5378e40] | 66 | np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), np.float32(thickness.value[f]), |
---|
| 67 | np.float32(length.value[j]), np.float32(radius.value[i]), np.float32(pars['scale']), |
---|
| 68 | np.float32(radius.weight[i]), np.float32(length.weight[j]), np.float32(thickness.weight[f]), |
---|
| 69 | np.float32(axis_theta.weight[k]), np.float32(axis_phi.weight[l]), np.float32(pars['core_sld']), |
---|
| 70 | np.float32(pars['shell_sld']), np.float32(pars['solvent_sld']),np.uint32(size), |
---|
| 71 | np.uint32(self.qx.size)) |
---|
[8a20be5] | 72 | cl.enqueue_copy(queue, self.res, self.res_b) |
---|
[5378e40] | 73 | |
---|
| 74 | sum += self.res |
---|
| 75 | vol += radius.weight[i]*length.weight[j]*thickness.weight[f]*pow(radius.value[i]+thickness.value[f],2)\ |
---|
| 76 | *(length.value[j]+2.0*thickness.value[f]) |
---|
| 77 | norm_vol += radius.weight[i]*length.weight[j]*thickness.weight[k] |
---|
| 78 | norm += radius.weight[i]*length.weight[j]*thickness.weight[f]*axis_theta.weight[k]\ |
---|
| 79 | *axis_phi.weight[l] |
---|
| 80 | |
---|
| 81 | if size>1: |
---|
| 82 | norm /= math.asin(1.0) |
---|
| 83 | if vol != 0.0 and norm_vol != 0.0: |
---|
| 84 | sum *= norm_vol/vol |
---|
| 85 | |
---|
| 86 | return sum/norm + pars['background'] |
---|