Changeset a42fec0 in sasmodels for Models

Timestamp:

Aug 4, 2014 5:20:07 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:

8cdb9f1

Parents:

099e053

Message:

Speed-up of 3X, compare.py working

Location:

Models

Files:

• code_capcyl.py (modified) (3 diffs)
• code_coreshellcyl.py (modified) (5 diffs)
• code_cylinder.py (modified) (6 diffs)
• code_ellipse.py (modified) (4 diffs)
• code_lamellar.py (modified) (1 diff)
• code_triaxialellipse.py (modified) (4 diffs)
• weights.py (modified) (1 diff)

Models/code_capcyl.py

@@ -3 +3 @@
 
 import numpy as np
-from math import asin, sqrt, fabs, atan
+from math import sqrt, fabs, atan
 import pyopencl as cl
 
 from weights import GaussianDispersion
-from sasmodel import card
+from sasmodel import card, set_precision
 
 
-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):
@@ -50 +37 @@
 
         _ctx,queue = card()
+        self.res[:] = 0
+        cl.enqueue_copy(queue, self.res_b, self.res)
         rad_cyl,length,rad_cap,theta,phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -79 +68 @@
                                         real(rad_cyl.weight[i]), real(length.weight[j]), real(theta.weight[k]), np.uint32(self.qx.size), np.uint32(size))
 
-                            cl.enqueue_copy(queue, self.res, self.res_b)
-
-                            sum += self.res
                             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]
 
-        if size > 1:
-            norm /= asin(1.0)
-
+        #if size > 1:
+         #   norm /= asin(1.0)
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        sum += self.res
         if vol != 0.0 and norm_vol != 0.0:
             sum *= norm_vol/vol

Models/code_coreshellcyl.py

@@ -6 +6 @@
 
 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
+from sasmodel import card, set_precision
 
 class GpuCoreShellCylinder(object):
@@ -46 +31 @@
 
         _ctx,queue = card()
+        self.res[:] = 0
+        cl.enqueue_copy(queue, self.res_b, self.res)
         radius, length, thickness, axis_phi, axis_theta = [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
                                      for base in GpuCoreShellCylinder.PD_PARS]
@@ -56 +43 @@
 
         sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
-
-        print radius.value
-        print thickness.weight
-        print axis_phi.weight
-        print axis_theta.weight
-        print length.value
-
         size = len(axis_theta.weight)
 
@@ -79 +59 @@
                                     real(pars['shell_sld']), real(pars['solvent_sld']),np.uint32(size),
                                     np.uint32(self.qx.size))
-                            cl.enqueue_copy(queue, self.res, self.res_b)
-
-                            sum += self.res
                             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])
@@ -90 +67 @@
         #if size>1:
          #   norm /= math.asin(1.0)
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        sum = self.res
         if vol != 0.0 and norm_vol != 0.0:
             sum *= norm_vol/vol

Models/code_cylinder.py

@@ -6 +6 @@
 
 from weights import GaussianDispersion
-from sasmodel import card
+from sasmodel import card, set_precision, set_precision_1d
 
-
-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
-
-def set_precision_1d(src, q, dtype):
-    q = np.ascontiguousarray(q, 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, q
 
 class GpuCylinder(object):
@@ -55 +28 @@
         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_b = cl.Buffer(ctx, cl.mem_flags.READ_WRITE, self.qx.nbytes)
         self.res = np.empty_like(self.qx)
 
@@ -61 +34 @@
 
         _ctx,queue = card()
+        self.res[:] = 0
+        cl.enqueue_copy(queue, self.res_b, self.res)
         radius, length, cyl_theta, cyl_phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -68 +43 @@
         radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
         length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
-        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -90, 180, False)
-        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -90, 180, False)
+        cyl_theta.value, cyl_theta.weight = cyl_theta.get_weights(pars['cyl_theta'], -np.inf, np.inf, False)
+        cyl_phi.value, cyl_phi.weight = cyl_phi.get_weights(pars['cyl_phi'], -np.inf, np.inf, False)
 
         #Perform the computation, with all weight points
@@ -87 +62 @@
                                            real(cyl_phi.weight[l]), real(cyl_theta.value[k]), real(cyl_phi.value[l]),
                                            np.uint32(self.qx.size), np.uint32(size))
-                        cl.enqueue_copy(queue, self.res, self.res_b)
-                        sum += self.res
+
                         vol += radius.weight[i]*length.weight[j]*pow(radius.value[i], 2)*length.value[j]
                         norm_vol += radius.weight[i]*length.weight[j]
@@ -95 +69 @@
        # if size > 1:
         #    norm /= math.asin(1.0)
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        sum = self.res
         if vol != 0.0 and norm_vol != 0.0:
             sum *= norm_vol/vol

Models/code_ellipse.py

@@ -6 +6 @@
 
 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
+from sasmodel import card, set_precision
 
 class GpuEllipse(object):
@@ -46 +31 @@
     #b_n = radius_b # want, a_n = radius_a # want, etc
         _ctx,queue = card()
+        self.res[:] = 0
+        cl.enqueue_copy(queue, self.res_b, self.res)
         radius_a, radius_b, axis_theta, axis_phi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -76 +63 @@
                                         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(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]
                         norm_vol += radius_a.weight[i]*radius_b.weight[j]
@@ -88 +73 @@
     #    if size > 1:
      #       norm /= math.asin(1.0)
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        sum += self.res
         if vol != 0.0 and norm_vol != 0.0:
             sum *= norm_vol/vol

Models/code_lamellar.py

@@ -5 +5 @@
 import pyopencl as cl
 from Models.weights import GaussianDispersion
-
-def set_precision(src, qx, qy, dtype):
-    qx = np.ascontiguousarray(qx, dtype=dtype)
-    qy = np.ascontiguousarray(qy, dtype=dtype)
-    if dtype == 'double':
-        header = """\
-#define real float
-"""
-    else:
-        header = """\
-#pragma OPENCL EXTENSION cl_khr_fp64: enable
-#define real double
-"""
-    return header+src, qx, qy
+from sasmodel import set_precision
 
 

Models/code_triaxialellipse.py

@@ -6 +6 @@
 
 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
+from sasmodel import card, set_precision
 
 class GpuTriEllipse:
@@ -45 +30 @@
 
         _ctx,queue = card()
+        self.res[:] = 0
+        cl.enqueue_copy(queue, self.res_b, self.res)
         axisA, axisB, axisC, theta, phi, psi = \
             [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
@@ -72 +59 @@
                                             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))
-                                cl.enqueue_copy(queue, self.res, self.res_b)
-                                sum += self.res
 
                                 vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*axisA.value[a]*axisB.value[b]*axisC.value[c]
@@ -81 +66 @@
       #  if size > 1:
        #     norm /= asin(1.0)
-
+        cl.enqueue_copy(queue, self.res, self.res_b)
+        sum = self.res
         if vol != 0.0 and norm_vol != 0.0:
             sum *= norm_vol/vol

Models/weights.py

@@ -19 +19 @@
         sigma = width * center if relative else width
         if sigma == 0:
-            return np.array([center, 1.], 'd')
+            return np.array([center],'d'), np.array([1.], 'd')
         x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts)
         x = x[(x >= min) & (x <= max)]