[dbb0048] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | |
---|
| 4 | import numpy as np |
---|
[a42fec0] | 5 | from math import sqrt, fabs, atan |
---|
[dbb0048] | 6 | import pyopencl as cl |
---|
[473183c] | 7 | |
---|
[dbb0048] | 8 | from weights import GaussianDispersion |
---|
[a42fec0] | 9 | from sasmodel import card, set_precision |
---|
| 10 | |
---|
| 11 | |
---|
[dbb0048] | 12 | |
---|
| 13 | class GpuCapCylinder(object): |
---|
| 14 | PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0, |
---|
| 15 | 'theta':0, 'phi':0} |
---|
| 16 | |
---|
| 17 | PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi'] |
---|
| 18 | |
---|
| 19 | def __init__(self, qx, qy, dtype='float32'): |
---|
| 20 | |
---|
| 21 | #create context, queue, and build program |
---|
| 22 | ctx,_queue = card() |
---|
[ca6c007] | 23 | trala = open('Kernel/NR_BessJ1.cpp').read()+"\n"+open('Kernel/Capcyl_Kfun.cpp').read()+"\n"+open('Kernel/Kernel-CapCyl.cpp').read() |
---|
[dbb0048] | 24 | src, qx, qy = set_precision(trala, qx, qy, dtype=dtype) |
---|
| 25 | self.prg = cl.Program(ctx, src).build() |
---|
| 26 | self.qx, self.qy = qx, qy |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | #buffers |
---|
| 30 | mf = cl.mem_flags |
---|
| 31 | self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx) |
---|
| 32 | self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy) |
---|
| 33 | self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes) |
---|
| 34 | self.res = np.empty_like(self.qx) |
---|
| 35 | |
---|
| 36 | def eval(self, pars): |
---|
| 37 | |
---|
| 38 | _ctx,queue = card() |
---|
[a42fec0] | 39 | self.res[:] = 0 |
---|
| 40 | cl.enqueue_copy(queue, self.res_b, self.res) |
---|
[dbb0048] | 41 | rad_cyl,length,rad_cap,theta,phi = \ |
---|
| 42 | [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma']) |
---|
| 43 | for base in GpuCapCylinder.PD_PARS] |
---|
| 44 | |
---|
| 45 | rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True) |
---|
| 46 | rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True) |
---|
| 47 | length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True) |
---|
| 48 | theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False) |
---|
| 49 | phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False) |
---|
| 50 | |
---|
| 51 | sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0 |
---|
| 52 | size = len(theta.weight) |
---|
| 53 | sub = pars['sld_capcyl']-pars['sld_solv'] |
---|
| 54 | real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64 |
---|
| 55 | |
---|
| 56 | for i in xrange(len(rad_cyl.weight)): |
---|
| 57 | for m in xrange(len(rad_cap.weight)): |
---|
| 58 | for j in xrange(len(length.weight)): |
---|
| 59 | for k in xrange(len(theta.weight)): |
---|
| 60 | for l in xrange(len(phi.weight)): |
---|
| 61 | hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i])) |
---|
| 62 | vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*length.value[j]+2.0*4.0*atan(1.0)/3.0\ |
---|
| 63 | *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist)) |
---|
| 64 | self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, |
---|
| 65 | real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]), |
---|
| 66 | real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']), |
---|
| 67 | real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]), |
---|
| 68 | real(rad_cyl.weight[i]), real(length.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size)) |
---|
| 69 | |
---|
| 70 | vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*vol_i |
---|
| 71 | norm_vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m] |
---|
| 72 | norm += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l] |
---|
| 73 | |
---|
[a42fec0] | 74 | #if size > 1: |
---|
| 75 | # norm /= asin(1.0) |
---|
| 76 | cl.enqueue_copy(queue, self.res, self.res_b) |
---|
| 77 | sum += self.res |
---|
[dbb0048] | 78 | if vol != 0.0 and norm_vol != 0.0: |
---|
| 79 | sum *= norm_vol/vol |
---|
| 80 | |
---|
| 81 | return sum/norm + pars['background'] |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | |
---|