Changeset 8faffcd in sasmodels

Timestamp:

Jul 10, 2014 4:06:29 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:

2de9a5e

Parents:

8a20be5

Message:

Update for Aaron

Files:

• Capcyl_Kfun.cpp (added)
• Kernel-CapCyl.cpp (modified) (2 diffs)
• Kernel-CoreShellCylinder.cpp (modified) (2 diffs)
• capcylcope.py (modified) (4 diffs)
• fit.py (modified) (3 diffs)
• fit2.py (modified) (1 diff)
• sasmodel.py (modified) (3 diffs)

Kernel-CapCyl.cpp

@@ -22 +22 @@
         float vaj = -1.0*hDist/rad_cap;
 
-        for(j=0;j<76;j++) //the value 76 was pre-set in the kernel...
+        for(j=0;j<76;j++) //the 76 corresponds to the Gauss constants
         {
-            zij = ( Gauss76Z[j]*(1.0-vaj) + vaj + 1.0 )/2.0;    //the "t" dummy
-            yyy = Gauss76Wt[j]*ConvLens_kernel(dp,q,zij,alpha);    //uses the same Kernel as the Dumbbell, here L>0
+            zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0;
+            yyy = Gauss76Wt[j]*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha);
             summj += yyy;
         }
@@ -35 +35 @@
         if(arg2 == 0) {be = 0.5;}
         else {
-            float ax=fabs(arg2);
-            if ((ax < 8.0) {
-                y=arg2*arg2;
-                ans1=arg2*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
-                ans=ans1/ans2;
-            } 
-            else {
-                y=64.0/(ax*ax);
-                xx=ax-2.356194491;
-                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
-                ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
-                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-(8.0/ax)*sin(xx)*ans2);
-                if (arg2 < 0.0) {ans *= -1;}
-                }
-                be = ans/arg2
+            be = NR_BessJ1(arg2)/arg2;
         }
 
-        if(arg1 == 0.0) {   //limiting value of sinc(0) is 1; sinc is not defined in math.h
+        if(arg1 == 0.0) {
             yyy += pi*rad_cyl*rad_cyl*length*2.0*be;
         } 

Kernel-CoreShellCylinder.cpp

@@ -22 +22 @@
 
         float si1=0; float si2=0; float be1=0; float be2=0;
-        float ax=0; float z=0; float xx=0; float y=0; float ans=0; float ans1=0; float ans2=0;
 
         float dr1 = core_sld-shell_sld;
@@ -37 +36 @@
         if (besarg1 == 0.0){be1 = 0.5;}
         else{
-            bel = NR_BessJ1(besarg1)/besarg1
+            be1 = NR_BessJ1(besarg1)/besarg1;
         }
         if (besarg2 == 0.0){be2 = 0.5;}
         else{
-            if((ax=fabs(besarg2)) < 8.0)
-            {
-                y=besarg2*besarg2;
-                ans1=besarg2*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
-                ans=ans1/ans2;
-            }
-            else
-            {
-                z=8.0/ax;
-                y=z*z;
-                xx=ax-2.356194491;
-                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
-                ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
-                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
-                if (besarg2 < 0.0) {ans *= -1;}
-            }
-            be2 = ans/besarg2;
+            be2 = NR_BessJ1(besarg2)/besarg2;
         }
         if (sinarg1 == 0.0){

capcylcope.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 GpuCapCylinder(object):
@@ -13 +28 @@
     PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
 
-    def __init__(self, qx, qy):
+    def __init__(self, qx, qy, dtype='float32'):
 
-        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)
+        ctx,_queue = card()
+        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(self.ctx, open('Kernel-CapCyl.cpp').read()).build()
+        self.prg = cl.Program(ctx, open('Kernel-CapCyl.cpp').read()).build()
 
         #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)
         self.vol_i = float(0.0)
-        self.vol_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, self.vol_i.nbytes)
+        self.vol_b = cl.Buffer(ctx, mf.WRITE_ONLY, self.vol_i.nbytes)
 
     def eval(self, pars):
 
+        _ctx,queue = card()
         rad_cyl,length,rad_cap,theta,phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -54 +69 @@
                         for l in xrange(len(phi.weight)):
 
-                            self.prg.CapCylinderKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                            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']),
@@ -60 +75 @@
                                         np.float32(rad_cyl.weight[i]), np.float32(length.weight[j]), np.uint32(self.qx.size), np.uint32(size))
 
-                            cl.enqueue_copy(self.queue, self.res, self.res_b)
-                            cl.enqueue_copy(self.queue, self.vol_i, self.vol_b)
+                            cl.enqueue_copy(queue, self.res, self.res_b)
+                            cl.enqueue_copy(queue, self.vol_i, self.vol_b)
 
                             sum += self.res

fit.py

@@ -12 +12 @@
 data = load_data('JUN03289.DAT')
 set_beam_stop(data, 0.004)
-
-
-
-
+"""
 model = SasModel(data, GpuCylinder, scale=1, radius=64.1, length=266.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=0,
                               cyl_theta=0, cyl_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,length_pd=0.1,
@@ -24 +21 @@
 model.cyl_theta.range(0,90)
 model.cyl_phi.range(0,90)
-"""
+
 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,
@@ -33 +30 @@
                  bi_thick_pd=0.1, bi_thick_pd_n=35, bi_thick_pd_nsigma=3, dtype='float')
 
-
-model = SasModel(data, GpuCoreShellCylinder, scale=1, radius=64.1, thickness=1, length=266.96, core_sld=1e-6, shell_sld=1e-6,
-                 solvent_sld=4e-6, background=0, axis_theta=0, axis_phi=0, radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
+"""
+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')
-"""
+
+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.scale.range(0,10)

fit2.py

@@ -34 +34 @@
 radial.cyl_phi.range(0,90)
 radial.scale.range(0,10)
+trans.cyl_theta = radial.cyl_theta + 90.
 
 
-fv = FreeVariables(names=['radial','trans'],
-                   cyl_theta=radial.cyl_theta,
-                   cyl_phi = radial.cyl_phi)
-problem = FitProblem([radial,trans], freevars=fv)
+problem = FitProblem([radial,trans])
 
 

sasmodel.py

@@ -62 +62 @@
 class SasModel(object):
     def __init__(self, data, model, dtype='float32', **kw):
+        self.__dict__['_parameters'] = {}
         self.index = data.mask==0
         self.iq = data.data[self.index]
@@ -79 +80 @@
         pars.update((base+'_pd_nsigma', 3) for base in model.PD_PARS)
         pars.update(kw)
-        self._parameters = dict((k, Parameter(v, name=k)) for k, v in pars.items())
+        self._parameters = dict((k, Parameter.default(v, name=k)) for k, v in pars.items())
 
     def numpoints(self):
@@ -90 +91 @@
         return self._parameters[par]
 
+    def __setattr__(self, par, val):
+        if par in self._parameters:
+            self._parameters[par] = val
+        else:
+            self.__dict__[par] = val
+
     def theory(self):
         pars = dict((k,v.value) for k,v in self._parameters.items())
-        print pars
         result = self.gpu.eval(pars)
         return result