Changeset 8a21ba3 in sasmodels

Timestamp:

Jul 11, 2014 1:40:30 PM (10 years ago)

Author:

HMP1 <helen.park@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

09e15be

Parents:

496b252

Message:

Fixed Ellipse code, update comparison file

Files:

• Kernel-CapCyl.cpp (modified) (4 diffs)
• Kernel-Ellipse.cpp (modified) (2 diffs)
• code_capcyl.py (modified) (2 diffs)
• code_cylinder.py (modified) (1 diff)
• code_ellipse.py (modified) (4 diffs)
• compare.py (modified) (3 diffs)
• fit.py (modified) (3 diffs)

Kernel-CapCyl.cpp

@@ -1 @@
-__kernel void CapCylinderKernel(__global const float *qx, __global const float *qy, __global float *_ptvalue, __global float *vol_i,
-const float rad_cyl, const float rad_cap, const float length, const float thet, const float ph, const float sub,
-const float scale, const float phi_weight, const float theta_float, const float rad_cap_weight, const float rad_cyl_weight,
-const float length_weight, const int total, const int size, __const float Gauss76Wt, __const float Gauss76Z)
+__kernel void CapCylinderKernel(__global const real *qx, __global const real *qy, __global real *_ptvalue, __global real *vol_i,
+const real rad_cyl, const real rad_cap, const real length, const real thet, const real ph, const real sub,
+const real scale, const real phi_weight, const real theta_float, const real rad_cap_weight, const real rad_cyl_weight,
+const real length_weight, const int total, const int size, __const real Gauss76Wt, __const real Gauss76Z)
 //ph is phi, sub is sldc-slds, thet is theta
 {
@@ -8 +8 @@
     if(i < total)
     {
-        float q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
-        float pi = 4.0*atan(1.0);
-        float theta = thet*pi/180.0;
-        float phi = ph*pi/180.0;
-        float cyl_x = cos(theta)*cos(phi);
-        float cyl_y = sin(theta);
-        float cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q;
-        float alpha = acos(cos_val);
-        float yyy=0; float ans1=0; float ans2=0; float y=0; float xx=0; float ans=0; float zij=0; float be=0; float summj=0;
+        real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
+        real pi = 4.0*atan(1.0);
+        real theta = thet*pi/180.0;
+        real phi = ph*pi/180.0;
+        real cyl_x = cos(theta)*cos(phi);
+        real cyl_y = sin(theta);
+        real cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q;
+        real alpha = acos(cos_val);
+        real yyy=0; real ans1=0; real ans2=0; real y=0; real xx=0; real ans=0; real zij=0; real be=0; real summj=0;
 
-        float hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl));
+        real hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl));
         vol_i[i] = pi*rad_cyl*rad_cyl*length+2.0*pi/3.0*((rad_cap-hDist)*(rad_cap-hDist)*(2*rad_cap+hDist));
-        float vaj = -1.0*hDist/rad_cap;
+        real vaj = -1.0*hDist/rad_cap;
 
         for(int j=0;j<76;j++) //the 76 corresponds to the Gauss constants
@@ -28 +28 @@
             summj += yyy;
         }
-        float inner = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap;
-        float arg1 = q*length/2.0*cos(alpha);
-        float arg2 = q*rad_cyl*sin(alpha);
+        real inner = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap;
+        real arg1 = q*length/2.0*cos(alpha);
+        real arg2 = q*rad_cyl*sin(alpha);
         yyy = inner;
 
@@ -44 +44 @@
             yyy += pi*rad_cyl*rad_cyl*length*sin(arg1)/arg1*2.0*be;
         }
-        float answer=yyy*yyy*1.0e8*sub*sub*scale/pi*rad_cyl*rad_cyl*length+2.0*pi*(2.0*rad_cap*rad_cap*rad_cap/3.0+rad_cap*rad_cap*hDist-hDist*hDist*hDist/3.0);
+        real answer=yyy*yyy*1.0e8*sub*sub*scale/pi*rad_cyl*rad_cyl*length+2.0*pi*(2.0*rad_cap*rad_cap*rad_cap/3.0+rad_cap*rad_cap*hDist-hDist*hDist*hDist/3.0);
         answer/=sin(alpha);
 

Kernel-Ellipse.cpp

@@ -1 @@
-__kernel void EllipsoidKernel(const float radius_a_weight, const float radius_b_weight, const float axis_theta_weight,
-const float axis_phi_weight, const float scale, const float radius_a, const float radius_b, const float sub, const float background, const float axis_theta, const float axis_phi, __global const float *qx,
-__global const float *qy, __global float *_ptvalue, const int length, const int size)
+__kernel void EllipsoidKernel(const real radius_a_weight, const real radius_b_weight, const real axis_theta_weight,
+const real axis_phi_weight, const real scale, const real radius_a, const real radius_b, const real sub, const real axis_theta, const real axis_phi, __global const real *qx,
+__global const real *qy, __global real *_ptvalue, const int length, const int size)
 {
      int i = get_global_id(0);
      if(i < length){
-         float ret = 0;
-         float q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
-         float pi = 4.0*atan(1.0);
-         float theta = axis_theta*pi/180.0;
-         float h = axis_phi*pi/180.0;
-         float cyl_x = cos(theta)*cos(h);
-         float cyl_y = sin(theta);
-         float cos_val = cyl_x*(qx[i]/q) + cyl_y*(qy[i]/q);
+         real ret = 0;
+         real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
+         real pi = 4.0*atan(1.0);
+         real theta = axis_theta*pi/180.0;
+         real cyl_x = cos(theta)*cos(axis_phi*pi/180.0);
+         real cyl_y = sin(theta);
+         real cos_val = cyl_x*(qx[i]/q) + cyl_y*(qy[i]/q);
 
-         float nu = radius_a/radius_b;
-         float arg = q*radius_b*sqrt(1.0+(cos_val*cos_val*((nu*nu)-1.0)));
+         real arg = q*radius_b*sqrt(1.0+(cos_val*cos_val*(((radius_a*radius_a/(radius_b*radius_b))-1.0))));
          if(arg == 0.0){
              ret = 1.0/3.0;
@@ -22 +20 @@
              ret = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
          }
-         ret*=ret*9.0*sub*sub;
-         ret*=(4.0/3.0*acos(-1.0)*radius_b*radius_b*radius_a)*scale*(1.0e8);
-         ret+=background;
+         ret*=ret*9.0*sub*sub*4.0/3.0*acos(-1.0)*radius_b*radius_b*radius_a*scale*(1.0e8);
 
          _ptvalue[i] = radius_a_weight*radius_b_weight*axis_theta_weight*radius_a*axis_phi_weight*ret*pow(radius_b, 2);
-
          if(size > 1){
             _ptvalue[i] *= fabs(cos(axis_theta*pi/180.0));

code_capcyl.py

@@ -66 +66 @@
         sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
         size = len(theta.weight)
-        sub = pars['sld_capcyl']-np.float32(['sld_solv'])
+        sub = pars['sld_capcyl']-pars['sld_solv']
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
 
         for i in xrange(len(rad_cyl.weight)):
@@ -75 +76 @@
 
                             self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                        self.vol_b, np.float32(rad_cyl.value[i]), np.float32(rad_cap.value[m]), np.float32(length.value[j]),
-                                        np.float32(theta.value[k]), np.float32(phi.value[l]), np.float32(sub), np.float32(pars['scale']),
-                                        np.float32(phi.weight[l]), np.float32(theta.weight[k]), np.float32(rad_cap.weight[m]),
-                                        np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size),
+                                        self.vol_b, real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]),
+                                        real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']),
+                                        real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]),
+                                        real(rad_cyl.weight[i]), real(length.weight[j]), np.uint32(self.qx.size), np.uint32(size),
                                         self.Gauss76W_b, self.Gauss76Z_b)
 

code_cylinder.py

@@ -68 +68 @@
                 for k in xrange(len(cyl_theta.weight)):
                     for l in xrange(len(cyl_phi.weight)):
-
-
                         self.prg.CylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(sub),
                                            real(radius.value[i]), real(length.value[j]), real(pars['scale']),

code_ellipse.py

@@ -6 +6 @@
 import pyopencl as cl
 from weights import GaussianDispersion
+from sasmodel import card
 
+def set_precision(src, qx, qy, dtype):
+    qx = np.ascontiguousarray(qx, dtype=dtype)
+    qy = np.ascontiguousarray(qy, dtype=dtype)
+    if np.dtype(dtype) == np.dtype('float32'):
+        header = """\
+#define real float
+"""
+    else:
+        header = """\
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
+#define real double
+"""
+    return header+src, qx, qy
 
 class GpuEllipse(object):
@@ -13 +27 @@
     }
     PD_PARS = ['radius_a', 'radius_b', 'axis_theta', 'axis_phi']
-    def __init__(self, qx, qy):
 
-        self.qx = np.asarray(qx, np.float32)
-        self.qy = np.asarray(qy, np.float32)
-        #create context, queue, and build program
-        self.ctx = cl.create_some_context()
-        self.queue = cl.CommandQueue(self.ctx)
-        self.prg = cl.Program(self.ctx, open('Kernel-Ellipse.cpp').read()).build()
+    def __init__(self, qx, qy, dtype='float32'):
+
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel-Ellipse.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
 
         #buffers
         mf = cl.mem_flags
-        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
-        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
-        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
+        self.qx_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
+        self.qy_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
+        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes)
         self.res = np.empty_like(self.qx)
 
     def eval(self, pars):
     #b_n = radius_b # want, a_n = radius_a # want, etc
+        _ctx,queue = card()
         radius_a, radius_b, axis_theta, axis_phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -43 +57 @@
         sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
         size = len(axis_theta.weight)
-        sub =  pars['sldEll'] - pars['sldSolv']
+        sub = pars['sldEll'] - pars['sldSolv']
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
 
         #Loop over radius weight points
@@ -54 +69 @@
                     for l in xrange(len(axis_phi.weight)):
                         #call the kernel
-                        self.prg.EllipsoidKernel(self.queue, self.qx.shape, None, np.float32(radius_a.weight[i]),
-                                        np.float32(radius_b.weight[j]), np.float32(axis_theta.weight[k]),
-                                        np.float32(axis_phi.weight[l]), np.float32(pars['scale']), np.float32(radius_a.value[i]),
-                                        np.float32(radius_b.value[j]), np.float32(sub),np.float32(pars['background']),
-                                        np.float32(axis_theta.value[k]), np.float32(axis_phi.value[l]), self.qx_b, self.qy_b,
-                                        self.res_b, np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
+                        self.prg.EllipsoidKernel(queue, self.qx.shape, None, real(radius_a.weight[i]),
+                                        real(radius_b.weight[j]), real(axis_theta.weight[k]),
+                                        real(axis_phi.weight[l]), real(pars['scale']), real(radius_a.value[i]),
+                                        real(radius_b.value[j]), real(sub), real(axis_theta.value[k]),
+                                        real(axis_phi.value[l]), self.qx_b, self.qy_b, self.res_b,
+                                        np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
                         #copy result back from buffer
-                        cl.enqueue_copy(self.queue, self.res, self.res_b)
+                        cl.enqueue_copy(queue, self.res, self.res_b)
                         sum += self.res
                         vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]

compare.py

@@ -38 +38 @@
     return theory
 
-def demo(N=1):
+def cyl(N=1):
     import sys
     import matplotlib.pyplot as plt
-    import numpy as np
 
     if len(sys.argv) > 1:
@@ -61 +60 @@
     cpu_time = toc()*1000./N
 
-    from cylcode import GpuCylinder
+    from code_cylinder import GpuCylinder
     model = SasModel(data, GpuCylinder, dtype='f', **pars)
+    tic()
+    for i in range(N):
+        gpu = model.theory()
+    gpu_time = toc()*1000./N
+
+    relerr = (gpu - cpu)/cpu
+    print "max(|(ocl-omp)/ocl|)", max(abs(relerr))
+    print "omp t=%.1f ms"%cpu_time
+    print "ocl t=%.1f ms"%gpu_time
+
+    plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time)
+    plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time)
+    plt.subplot(133); plot_data(data, 1e8*relerr); plt.title("relerr x 10^8"); plt.colorbar()
+    plt.show()
+
+def ellipse(N=1):
+    import sys
+    import matplotlib.pyplot as plt
+
+    if len(sys.argv) > 1:
+        N = int(sys.argv[1])
+    data = load_data('JUN03289.DAT')
+    set_beam_stop(data, 0.004)
+
+    pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9,
+                axis_theta=0, axis_phi=90, radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, radius_b_pd=0.1, radius_b_pd_n=10,
+                radius_b_pd_nsigma=3, axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, axis_phi_pd=0.1,
+                axis_phi_pd_n=6, axis_phi_pd_nsigma=3,)
+
+    model = sasview_model('ellipsoid', **pars)
+    tic()
+    for i in range(N):
+        cpu = sasview_eval(model, data)
+    cpu_time = toc()*1000./N
+
+    from code_ellipse import GpuEllipse
+    model = SasModel(data, GpuEllipse, dtype='f', **pars)
     tic()
     for i in range(N):
@@ -80 +116 @@
 
 if __name__ == "__main__":
-    demo()
+    ellipse()
 
 

fit.py

@@ -25 +25 @@
 """
 
-"""
+
 model = SasModel(data, GpuEllipse, scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9,
                  axis_theta=0, axis_phi=90, radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, radius_b_pd=0.1, radius_b_pd_n=10,
                  radius_b_pd_nsigma=3, axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, axis_phi_pd=0.1,
                  axis_phi_pd_n=6, axis_phi_pd_nsigma=3, dtype='float')
-"""
+
 
 """
@@ -37 +37 @@
 """
 
-
+"""
 model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=4e-6,
                  solvent_sld=1e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
@@ -43 +43 @@
                  axis_theta_pd=0.1, axis_theta_pd_n=2, axis_theta_pd_nsigma=3, axis_phi_pd=0.1, axis_phi_pd_n=2,
                  axis_phi_pd_nsigma=3, dtype='float')
-
+"""