| 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, 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 | real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 |
|---|
| 60 | for i in xrange(len(radius.weight)): |
|---|
| 61 | for j in xrange(len(length.weight)): |
|---|
| 62 | for k in xrange(len(axis_theta.weight)): |
|---|
| 63 | for l in xrange(len(axis_phi.weight)): |
|---|
| 64 | for f in xrange(len(thickness.weight)): |
|---|
| 65 | |
|---|
| 66 | self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, |
|---|
| 67 | real(axis_theta.value[k]), real(axis_phi.value[l]), real(thickness.value[f]), |
|---|
| 68 | real(length.value[j]), real(radius.value[i]), real(pars['scale']), |
|---|
| 69 | real(radius.weight[i]), real(length.weight[j]), real(thickness.weight[f]), |
|---|
| 70 | real(axis_theta.weight[k]), real(axis_phi.weight[l]), real(pars['core_sld']), |
|---|
| 71 | real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size), |
|---|
| 72 | np.uint32(self.qx.size)) |
|---|
| 73 | cl.enqueue_copy(queue, self.res, self.res_b) |
|---|
| 74 | |
|---|
| 75 | sum += self.res |
|---|
| 76 | vol += radius.weight[i]*length.weight[j]*thickness.weight[f]*pow(radius.value[i]+thickness.value[f],2)\ |
|---|
| 77 | *(length.value[j]+2.0*thickness.value[f]) |
|---|
| 78 | norm_vol += radius.weight[i]*length.weight[j]*thickness.weight[k] |
|---|
| 79 | norm += radius.weight[i]*length.weight[j]*thickness.weight[f]*axis_theta.weight[k]\ |
|---|
| 80 | *axis_phi.weight[l] |
|---|
| 81 | |
|---|
| 82 | if size>1: |
|---|
| 83 | norm /= math.asin(1.0) |
|---|
| 84 | if vol != 0.0 and norm_vol != 0.0: |
|---|
| 85 | sum *= norm_vol/vol |
|---|
| 86 | |
|---|
| 87 | return sum/norm + pars['background'] |
|---|
