Changeset 1726b21 in sasmodels for Models

Timestamp:

Aug 8, 2014 1:45:51 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:

ae7d639

Parents:

8cdb9f1

Message:

cylinder now MUCH faster!

Location:

Models

Files:

• code_capcyl.py (modified) (4 diffs)
• code_coreshellcyl.py (modified) (2 diffs)
• code_cylinder.py (modified) (5 diffs)
• code_cylinder_f.py (added)
• code_ellipse.py (modified) (2 diffs)
• code_triaxialellipse.py (modified) (2 diffs)

Models/code_capcyl.py

@@ -12 +12 @@
 
 class GpuCapCylinder(object):
-    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
+    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'len_cyl':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
              'theta':0, 'phi':0}
 
-    PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
+    PD_PARS = ['rad_cyl', 'len_cyl', 'rad_cap', 'theta', 'phi']
 
     def __init__(self, qx, qy, dtype='float32'):
@@ -39 +39 @@
         self.res[:] = 0
         cl.enqueue_copy(queue, self.res_b, self.res)
-        rad_cyl,length,rad_cap,theta,phi = \
+
+        rad_cyl,len_cyl,rad_cap,theta,phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
              for base in GpuCapCylinder.PD_PARS]
@@ -45 +46 @@
         rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True)
         rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True)
-        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+        len_cyl.value, len_cyl.weight = len_cyl.get_weights(pars['len_cyl'], 0, 10000, True)
         theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
         phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
@@ -56 +57 @@
         for i in xrange(len(rad_cyl.weight)):
             for m in xrange(len(rad_cap.weight)):
-                for j in xrange(len(length.weight)):
+                for j in xrange(len(len_cyl.weight)):
+
+                    hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
+                    vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*len_cyl.value[j]+2.0*4.0*atan(1.0)/3.0\
+                                    *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
+                    vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*vol_i
+                    norm_vol += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]
+
                     for k in xrange(len(theta.weight)):
                         for l in xrange(len(phi.weight)):
-                            hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
-                            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\
-                                            *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
+
                             self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]),
+                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(len_cyl.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]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
+                                        real(rad_cyl.weight[i]), real(len_cyl.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
 
-                            vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*vol_i
-                            norm_vol += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]
-                            norm += rad_cyl.weight[i]*length.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l]
+                            norm += rad_cyl.weight[i]*len_cyl.weight[j]*rad_cap.weight[m]*theta.weight[k]*phi.weight[l]
+
+
 
         #if size > 1:

Models/code_coreshellcyl.py

@@ -48 +48 @@
         for r in xrange(len(radius.weight)):
             for l in xrange(len(length.weight)):
-                for at in xrange(len(axis_theta.weight)):
-                    for p in xrange(len(axis_phi.weight)):
-                        for th in xrange(len(thickness.weight)):
+                for th in xrange(len(thickness.weight)):
+
+                    vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\
+                                   *(length.value[l]+2.0*thickness.value[th])
+                    norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
+
+                    for at in xrange(len(axis_theta.weight)):
+                        for p in xrange(len(axis_phi.weight)):
 
                             self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
@@ -59 +64 @@
                                     real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size),
                                     np.uint32(self.qx.size))
-                            vol += radius.weight[r]*length.weight[l]*thickness.weight[th]*pow(radius.value[r]+thickness.value[th],2)\
-                                   *(length.value[l]+2.0*thickness.value[th])
-                            norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
+
                             norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\
                                     *axis_phi.weight[p]

Models/code_cylinder.py

@@ -6 +6 @@
 
 from weights import GaussianDispersion
-from sasmodel import card, set_precision, set_precision_1d
+from sasmodel import card, set_precision, set_precision_1d, tic, toc
 
 
@@ -33 +33 @@
     def eval(self, pars):
 
+        tic()
         _ctx,queue = card()
         self.res[:] = 0
@@ -55 +56 @@
         for i in xrange(len(radius.weight)):
             for j in xrange(len(length.weight)):
+
+                vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
+                norm_vol += radius.weight[i]*length.weight[j]
+
                 for k in xrange(len(cyl_theta.weight)):
                     for l in xrange(len(cyl_phi.weight)):
@@ -63 +68 @@
                                            np.uint32(self.qx.size), np.uint32(size))
 
-                        vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
-                        norm_vol += radius.weight[i]*length.weight[j]
                         norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l]
+
 
        # if size > 1:
@@ -74 +78 @@
             sum *= norm_vol/vol
 
+        print toc()*1000, self.qx.shape[0]
         return sum/norm+pars['background']
 

Models/code_ellipse.py

@@ -53 +53 @@
             for j in xrange(len(radius_b.weight)):
                 #Average over theta distribution
+                vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]
+                norm_vol += radius_a.weight[i]*radius_b.weight[j]
+
                 for k in xrange(len(axis_theta.weight)):
                     #Average over phi distribution
@@ -64 +67 @@
                                         np.uint32(self.qx.size), np.uint32(len(axis_theta.weight)))
 
-                        vol += radius_a.weight[i]*radius_b.weight[j]*pow(radius_b.value[j], 2)*radius_a.value[i]
-                        norm_vol += radius_a.weight[i]*radius_b.weight[j]
                         norm += radius_a.weight[i]*radius_b.weight[j]*axis_theta.weight[k]*axis_phi.weight[l]
+
         # Averaging in theta needs an extra normalization
         # factor to account for the sin(theta) term in the

Models/code_triaxialellipse.py

@@ -9 +9 @@
 
 class GpuTriEllipse:
-    PARS = {'scale':1, 'axisA':35, 'axisB':100, 'axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
-            'theta':0, 'phi':0, 'psi':0}
+    PARS = {'scale':1, 'semi_axisA':35, 'semi_axisB':100, 'semi_axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
+            'axis_theta':0, 'axis_phi':0, 'axis_psi':0}
 
-    PD_PARS = ['axisA', 'axisB', 'axisC', 'theta', 'phi', 'psi']
+    PD_PARS = ['semi_axisA', 'semi_axisB', 'semi_axisC', 'axis_theta', 'axis_phi', 'axis_psi']
 
     def __init__(self, qx, qy, dtype='float32'):
@@ -32 +32 @@
         self.res[:] = 0
         cl.enqueue_copy(queue, self.res_b, self.res)
-        axisA, axisB, axisC, theta, phi, psi = \
+        semi_axisA, semi_axisB, semi_axisC, axis_theta, axis_phi, axis_psi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
              for base in GpuTriEllipse.PD_PARS]
 
-        axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 10000, True)
-        axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 10000, True)
-        axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 10000, True)
-        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
-        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
-        psi.value, psi.weight = psi.get_weights(pars['psi'], -90, 180, False)
+        semi_axisA.value, semi_axisA.weight = semi_axisA.get_weights(pars['semi_axisA'], 0, 10000, True)
+        semi_axisB.value, semi_axisB.weight = semi_axisB.get_weights(pars['semi_axisB'], 0, 10000, True)
+        semi_axisC.value, semi_axisC.weight = semi_axisC.get_weights(pars['semi_axisC'], 0, 10000, True)
+        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
+        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
+        axis_psi.value, axis_psi.weight = axis_psi.get_weights(pars['axis_psi'], -90, 180, False)
 
         sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-        size = len(theta.weight)
+        size = len(axis_theta.weight)
         sub = pars['sldEll'] - pars['sldSolv']
 
         real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
-        for a in xrange(len(axisA.weight)):
-            for b in xrange(len(axisB.weight)):
-                for c in xrange(len(axisC.weight)):
-                    for t in xrange(len(theta.weight)):
-                        for i in xrange(len(phi.weight)):
-                            for s in xrange(len(psi.weight)):
+        for a in xrange(len(semi_axisA.weight)):
+            for b in xrange(len(semi_axisB.weight)):
+                for c in xrange(len(semi_axisC.weight)):
+
+                    vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*semi_axisA.value[a]*semi_axisB.value[b]*semi_axisC.value[c]
+                    norm_vol += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]
+
+                    for t in xrange(len(axis_theta.weight)):
+                        for i in xrange(len(axis_phi.weight)):
+                            for s in xrange(len(axis_psi.weight)):
                                 self.prg.TriaxialEllipseKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
-                                            real(sub), real(pars['scale']), real(axisA.value[a]), real(axisB.value[b]),
-                                            real(axisC.value[c]), real(phi.value[i]), real(theta.value[t]), real(psi.value[s]),
-                                            real(axisA.weight[a]), real(axisB.weight[b]), real(axisC.weight[c]), real(psi.weight[s]),
-                                            real(phi.weight[i]), real(theta.weight[t]), np.uint32(self.qx.size), np.uint32(size))
+                                            real(sub), real(pars['scale']), real(semi_axisA.value[a]), real(semi_axisB.value[b]),
+                                            real(semi_axisC.value[c]), real(axis_phi.value[i]), real(axis_theta.value[t]), real(axis_psi.value[s]),
+                                            real(semi_axisA.weight[a]), real(semi_axisB.weight[b]), real(semi_axisC.weight[c]), real(axis_psi.weight[s]),
+                                            real(axis_phi.weight[i]), real(axis_theta.weight[t]), np.uint32(self.qx.size), np.uint32(size))
 
-                                vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*axisA.value[a]*axisB.value[b]*axisC.value[c]
-                                norm_vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]
-                                norm += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*theta.weight[t]*phi.weight[i]*psi.weight[s]
+                                norm += semi_axisA.weight[a]*semi_axisB.weight[b]*semi_axisC.weight[c]*axis_theta.weight[t]*axis_phi.weight[i]*axis_psi.weight[s]
+
 
       #  if size > 1: