Changeset 2de9a5e in sasmodels

Timestamp:

Jul 11, 2014 2:47:47 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:

be5d7df

Parents:

8faffcd

Message:

Update for Aaron

Files:

• Capcyl_Gauss.py (added)
• Capcyl_Kfun.cpp (modified) (3 diffs)
• Kernel-CapCyl.cpp (modified) (4 diffs)
• Kernel-CoreShellCylinder.cpp (modified) (3 diffs)
• Kernel-Lamellar.cpp (modified) (2 diffs)
• Kernel-TriaxialEllipse.cpp (added)
• code_capcyl.py (moved) (moved from capcylcope.py) (4 diffs)
• code_coreshellcyl.py (moved) (moved from coreshellcylcode.py)
• code_cylinder.py (moved) (moved from cylcode.py)
• code_ellipse.py (moved) (moved from ellipsecode.py)
• code_lamellar.py (moved) (moved from lamellarcode.py) (3 diffs)
• code_triaxialellipse.py (added)
• fit.py (modified) (2 diffs)

Capcyl_Kfun.cpp

@@ -14 +14 @@
         return(val);
 }
-const real Gauss76Z[76]={
+/*
+real Gauss76Z[76]={
          .999505948362153*(-1.0),               //0
          .997397786355355*(-1.0),
@@ -92 +93 @@
         -.999505948362153*(-1.0)                //75
 };
-const real Gauss76Wt[76]={
+real Gauss76Wt[76]={
         .00126779163408536,             //0
         .00294910295364247,
@@ -170 +171 @@
         .00126779163408536              //75 (indexed from 0)
 };
+*/

Kernel-CapCyl.cpp

@@ -1 @@
-__kernel void CapCylinderKernel(__global const float *qx, __global const float *qy, __global float *__ptvalue, __global float *vol_i
+__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 length_weight, const int total, const int size, __const float Gauss76Wt, __const float 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 q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
         float pi = 4.0*atan(1.0);
         float theta = thet*pi/180.0;
@@ -16 +16 @@
         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 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;
 
         float hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl));
-        vol_i[i] = pi*rad_cyl*rad_cyl*len_cyl+2.0*pi/3.0*((rad_cap-hDist)*(rad_cap-hDist)*(2*rad_cap+hDist));
+        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;
 
-        for(j=0;j<76;j++) //the 76 corresponds to the Gauss constants
+        for(int j=0;j<76;j++) //the 76 corresponds to the Gauss constants
         {
             zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0;
@@ -45 +45 @@
         }
         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);
-        answer/=sin(alpha)
+        answer/=sin(alpha);
 
         _ptvalue[i] = rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i[i]*answer

Kernel-CoreShellCylinder.cpp

@@ -11 +11 @@
         float pi = 4.0*atan(1.0);
         float theta = axis_theta*pi/180.0;
-        float phi = axis_phi*pi/180.0;
-        float cyl_x = cos(theta)*cos(phi);
+        float cyl_x = cos(theta)*cos(axis_phi*pi/180.0);
         float cyl_y = sin(theta);
-        float cos_val = cyl_x*qx[i]/q + cyl_y*qx[i]/q;
-        float alpha = acos(cos_val);
+        float alpha = acos(cyl_x*qx[i]/q + cyl_y*qx[i]/q);
 
         if (alpha == 0.0){
@@ -23 +21 @@
         float si1=0; float si2=0; float be1=0; float be2=0;
 
-        float dr1 = core_sld-shell_sld;
-        float dr2 = shell_sld-solvent_sld;
-        float vol1 = pi*radius*radius*(length);
         float vol2 = pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
 
@@ -54 +49 @@
             si2 = sin(sinarg2)/sinarg2;
         }
+        float tt = 2.0*vol2*(shell_sld-solvent_sld)*si2*be2 + 2.0*(pi*radius*radius*(length))*(core_sld-shell_sld)*si1*be1;
 
-        float t1 = 2.0*vol1*dr1*si1*be1;
-        float t2 = 2.0*vol2*dr2*si2*be2;
-
-        float answer = ((t1+t2)*(t1+t2))*sin(alpha)/fabs(sin(alpha));
+        float answer = (tt*tt)*sin(alpha)/fabs(sin(alpha));
         float vol=pi*(radius+thickness)*(radius+thickness)*(length+2.0*thickness);
         answer = answer/vol*1.0e8*scale;

Kernel-Lamellar.cpp

@@ -1 @@
-#ifndef real
-# define real float
-#endif
 
 __kernel void LamellarKernel(__global const real *qx, global const real *qy, __global real *ret, const real bi_thick,
- const real scale, const real sub, const real background, const int length)
+ const real scale, const real sub, const int length)
 {
     int i = get_global_id(0);
@@ -12 +9 @@
         real pi = 4.0*atan(1.0);
         real Pq = 2.0*sub*(sub/q)/q*(1.0-cos(q*bi_thick));
-        ret[i] = 2.0*pi*scale*Pq/(q*q)/bi_thick*1.0e8;
-        ret[i] += background;
+        ret[i] = 2.0*pi*scale*Pq/(q*q)/bi_thick*1.0e8 ;
     }
 }

code_capcyl.py

@@ -7 +7 @@
 from weights import GaussianDispersion
 from sasmodel import card
+from Capcyl_Gauss import Gauss76Wt, Gauss76Z
 
 def set_precision(src, qx, qy, dtype):
@@ -34 +35 @@
         trala = open('NR_BessJ1.cpp').read()+"\n"+open('Capcyl_Kfun.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read()
         src, qx, qy = set_precision(trala, qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
 
-        self.prg = cl.Program(ctx, open('Kernel-CapCyl.cpp').read()).build()
 
         #buffers
@@ -41 +43 @@
         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)
+        G, Z = Gauss76Wt, Gauss76Z
+        self.Gauss76W_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=G)
+        self.Gauss76Z_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=Z)
         self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, qx.nbytes)
         self.res = np.empty_like(self.qx)
@@ -73 +78 @@
                                         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))
+                                        np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size),
+                                        self.Gauss76W_b, self.Gauss76Z_b)
 
                             cl.enqueue_copy(queue, self.res, self.res_b)

code_lamellar.py

@@ -11 +11 @@
     if dtype == 'double':
         header = """\
+#define real float
+"""
+    else:
+        header = """\
 #pragma OPENCL EXTENSION cl_khr_fp64: enable
 #define real double
 """
-        return header+src,qx,qy
-    else:
-        return src,qx,qy
+    return header+src, qx, qy
 
 
@@ -23 +25 @@
         'scale':1, 'bi_thick':1, 'sld_bi':1e-6, 'sld_sol':0, 'background':0,
     }
-
-    def __init__(self, qx, qy, dtype='float'):
+    PD_PARS = {'bi_thick'}
+    def __init__(self, qx, qy, dtype='float32'):
 
         #create context, queue, and build program
@@ -51 +53 @@
         for i in xrange(len(bi_thick.weight)):
             self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(bi_thick.value[i]),
-                                    real(pars['scale']), real(sub), real(pars['background']), np.uint32(self.qx.size))
+                                    real(pars['scale']), real(sub), np.uint32(self.qx.size))
             cl.enqueue_copy(self.queue, self.res, self.res_b)
 

fit.py

@@ -3 +3 @@
 
 from bumps.names import *
-from cylcode import GpuCylinder
-from lamellarcode import GpuLamellar
-from ellipsecode import GpuEllipse
-from coreshellcylcode import GpuCoreShellCylinder
+from code_cylinder import GpuCylinder
+from code_lamellar import GpuLamellar
+from code_ellipse import GpuEllipse
+from code_coreshellcyl import GpuCoreShellCylinder
+from code_capcyl import GpuCapCylinder
+from code_triaxialellipse import GpuTriEllipse
 from sasmodel import SasModel, load_data, set_beam_stop
 
@@ -25 +27 @@
                  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')
+                 axis_phi_pd_n=6, axis_phi_pd_nsigma=3, dtype='float32')
 
 model = SasModel(data, GpuLamellar, scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0,
-                 bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float')
+                 bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float32')
 
-"""
+
 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,
                  length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, thickness_pd=0.1, thickness_pd_n=2, thickness_pd_nsigma=3,
                  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')
+                 axis_phi_pd_nsigma=3, dtype='float32')
 
-model = SasModel(data, GpuCapCylinder, scale, rad_cyl, rad_cap, length, sld_capcyl, sld_solv, background, theta, phi,
-                 rad_cyl_pd, rad_cyl_pd_n, rad_cyl_nsigma, rad_cap_pd, rad_cap_pd_n, rad_cap_pd_nsigma, length_)
+model = SasModel(data, GpuCapCylinder, scale=1, rad_cyl=20, rad_cap=40, length=400, sld_capcyl=1e-6, sld_solv=6.3e-6,
+                 background=0, theta=0, phi=0, rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_nsigma=3, rad_cap_pd=.1, rad_cap_pd_n=1,
+                 rad_cap_pd_nsigma=3, length_pd=.1, length_pd_n=10, length_pd_nsigma=3, theta_pd=.1, theta_pd_n=4,
+                 theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, dtype='float32')
+"""
+model = SasModel(data, GpuTriEllipse, scale=1, axisA=35, axisB=100, axisC=400, sldEll=1e-6, sldSolv=6.3e-6,
+                 background=0, theta=57, phi=57, psi=0, theta_pd=.1, theta_pd_n=4,
+                 theta_pd_nsigma=3, phi_pd=.1, phi_pd_n=4, phi_pd_nsigma=3, psi_pd=.1, psi_pd_n=4, psi_pd_nsigma=3,
+                 axisA_pd=.1, axisA_pd_n=4, axisA_pd_nsigma=3, axisB_pd=.1, axisB_pd_n=4, axisB_pd_nsigma=3,
+                 axisC_pd=.1, axisC_pd_n=4, axisC_pd_nsigma=3, dtype='float')
 
-model.scale.range(0,10)
+model.scale.range(0,1)
 
 problem = FitProblem(model)