Changeset 750ffa5 in sasmodels

Timestamp:

Mar 9, 2015 2:04:55 PM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

3a45c2c

Parents:

48f0194

Message:

allow test of dll using single precision

Files:

• compare.py (modified) (1 diff)
• sasmodels/bumps_model.py (modified) (1 diff)
• sasmodels/kernel_template.c (modified) (7 diffs)
• sasmodels/kernelcl.py (modified) (4 diffs)
• sasmodels/kerneldll.py (modified) (9 diffs)
• sasmodels/kernelpy.py (modified) (1 diff)

compare.py

@@ -14 +14 @@
 from sasmodels import kerneldll
 from sasmodels.convert import revert_model
+kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
 
 # List of available models

sasmodels/bumps_model.py

@@ -365 +365 @@
         if 'theory' not in self._cache:
             if self._fn is None:
-                input_value = self.model.make_input(self._fn_inputs)
-                self._fn = self.model(input_value)
+                q_input= self.model.make_input(self._fn_inputs)
+                self._fn = self.model(q_input)
 
             fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars]

sasmodels/kernel_template.c

@@ -12 +12 @@
      #include <cstdio>
      #include <cmath>
+     using namespace std;
      #if defined(_MSC_VER)
-     #define kernel extern "C" __declspec( dllexport )
+     #   define kernel extern "C" __declspec( dllexport )
+         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
      #else
-     #define kernel extern "C"
+     #   define kernel extern "C"
      #endif
-     using namespace std;
-     inline void SINCOS(double angle, double &svar, double &cvar)
-       { svar=sin(angle); cvar=cos(angle); }
+     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
 #  else
      #include <stdio.h>
-     #include <math.h>
-     #if defined(_MSC_VER)
-     #define kernel __declspec( dllexport )
-     #else
+     #include <tgmath.h> // C99 type-generic math, so sin(float) => sinf
+     // MSVC doesn't support C99, so no need for dllexport on C99 branch
      #define kernel
-     #endif
-     #define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
 #  endif
-     #if defined(_MSC_VER)
-     inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
-     #endif
 #  define global
 #  define local
 #  define constant const
+// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+// OpenCL pown(a,b) = C99 pow(a,b), b integer
 #  define powr(a,b) pow(a,b)
 #  define pown(a,b) pow(a,b)
@@ -42 +38 @@
 #    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
 #  else
-#    define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0)
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
 #  endif
 #endif
@@ -126 +122 @@
     const double weight = IQ_WEIGHT_PRODUCT;
     if (weight > cutoff) {
-      const double I = Iq(qi, IQ_PARAMETERS);
-      if (I>=0.0) { // scattering cannot be negative
-        ret += weight*I;
+      const double scattering = Iq(qi, IQ_PARAMETERS);
+      if (scattering >= 0.0) { // scattering cannot be negative
+        ret += weight*scattering;
         norm += weight;
       #ifdef VOLUME_PARAMETERS
@@ -136 +132 @@
       #endif
       }
-    //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,I); }
+    //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }
     IQ_CLOSE_LOOPS
@@ -199 +195 @@
     if (weight > cutoff) {
 
-      const double I = Iqxy(qxi, qyi, IQXY_PARAMETERS);
-      if (I>=0.0) { // scattering cannot be negative
+      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
+      if (scattering >= 0.0) { // scattering cannot be negative
         // TODO: use correct angle for spherical correction
         // Definition of theta and phi are probably reversed relative to the
@@ -211 +207 @@
         const double spherical_correction
           = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0);
-        ret += spherical_correction * weight * I;
+        ret += spherical_correction * weight * scattering;
       #else
-        ret += weight * I;
+        ret += weight * scattering;
       #endif
         norm += weight;
@@ -222 +218 @@
         norm_vol += vol_weight;
       }
-      //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,I); }
+      //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }
     IQXY_CLOSE_LOOPS

sasmodels/kernelcl.py

@@ -68 +68 @@
     #open(info['name']+'.cl','w').write(source)
     #source = open(info['name']+'.cl','r').read()
-    return GpuModel(source, info, dtype)
+    return GpuModel(source, info, np.dtype(dtype))
 
 ENV = None
@@ -92 +92 @@
     Return the size of an execution batch for *kernel* running on *queue*.
     """
-    return kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
-                                      queue.device)
+    return kernel.get_work_group_info(
+            cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
+            queue.device)
 
 def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
@@ -219 +220 @@
         self.info = info
         self.source = source
-        self.dtype = dtype
+        self.dtype = np.dtype(dtype)
         self.program = None # delay program creation
 
@@ -230 +231 @@
         self.__dict__ = state.copy()
 
-    def __call__(self, input_value):
-        if self.dtype != input_value.dtype:
-            raise TypeError("data and kernel have different types")
+    def __call__(self, q_input):
+        if self.dtype != q_input.dtype:
+            raise TypeError("data is %s kernel is %s" % (q_input.dtype, self.dtype))
         if self.program is None:
             compiler = environment().compile_program
             self.program = compiler(self.info['name'], self.source, self.dtype)
-        kernel_name = generate.kernel_name(self.info, input_value.is_2D)
+        kernel_name = generate.kernel_name(self.info, q_input.is_2D)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, input_value)
+        return GpuKernel(kernel, self.info, q_input)
 
     def release(self):

sasmodels/kerneldll.py

@@ -1 +1 @@
 """
 C types wrapper for sasview models.
+
+The global attribute *ALLOW_SINGLE_PRECISION_DLLS* should be set to *True* if
+you wish to allow single precision floating point evaluation for the compiled
+models, otherwise it defaults to *False*.
 """
+
 import sys
 import os
 import tempfile
 import ctypes as ct
-from ctypes import c_void_p, c_int, c_double
+from ctypes import c_void_p, c_int, c_double, c_float
 
 import numpy as np
@@ -36 +41 @@
 DLL_PATH = tempfile.gettempdir()
 
-
-def dll_path(info):
+ALLOW_SINGLE_PRECISION_DLLS = False
+
+
+def dll_path(info, dtype="double"):
     """
     Path to the compiled model defined by *info*.
@@ -43 +50 @@
     from os.path import join as joinpath, split as splitpath, splitext
     basename = splitext(splitpath(info['filename'])[1])[0]
+    if np.dtype(dtype) == generate.F32:
+        basename += "32"
     return joinpath(DLL_PATH, basename+'.so')
 
 
-def load_model(kernel_module, dtype=None):
+def load_model(kernel_module, dtype="double"):
     """
     Load the compiled model defined by *kernel_module*.
@@ -57 +66 @@
     be defined without using too many resources.
     """
+    if not ALLOW_SINGLE_PRECISION_DLLS: dtype = "double"   # Force 64-bit dll
+    dtype = np.dtype(dtype)
+
     source, info = generate.make(kernel_module)
     if callable(info.get('Iq',None)):
         return PyModel(info)
+
+    if dtype == generate.F32: # 32-bit dll
+        source = generate.use_single(source)
+        tempfile_prefix = 'sas_'+info['name']+'32_'
+    else:
+        tempfile_prefix = 'sas_'+info['name']+'_'
+
     source_files = generate.sources(info) + [info['filename']]
+    dllpath = dll_path(info, dtype)
     newest = max(os.path.getmtime(f) for f in source_files)
-    dllpath = dll_path(info)
     if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest:
         # Replace with a proper temp file
-        fid, filename = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name'])
+        fid, filename = tempfile.mkstemp(suffix=".c",prefix=tempfile_prefix)
         os.fdopen(fid,"w").write(source)
         command = COMPILE%{"source":filename, "output":dllpath}
@@ -76 +95 @@
             #os.unlink(filename); print "saving compiled file in %r"%filename
             pass
-    return DllModel(dllpath, info)
+    return DllModel(dllpath, info, dtype=dtype)
 
 
@@ -96 +115 @@
     Call :meth:`release` when done with the kernel.
     """
-    def __init__(self, dllpath, info):
+    def __init__(self, dllpath, info, dtype=generate.F32):
         self.info = info
         self.dllpath = dllpath
         self.dll = None
+        self.dtype = np.dtype(dtype)
 
     def _load_dll(self):
@@ -110 +130 @@
         self.dll = ct.CDLL(self.dllpath)
 
-        pd_args_1d = [c_void_p, c_double] + [c_int]*Npd1d if Npd1d else []
-        pd_args_2d= [c_void_p, c_double] + [c_int]*Npd2d if Npd2d else []
+        fp = c_float if self.dtype == generate.F32 else c_double
+        pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
+        pd_args_2d= [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
         self.Iq = self.dll[generate.kernel_name(self.info, False)]
-        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [c_double]*Nfixed1d
+        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d
 
         self.Iqxy = self.dll[generate.kernel_name(self.info, True)]
-        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [c_double]*Nfixed2d
+        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d
 
     def __getstate__(self):
@@ -125 +146 @@
 
     def __call__(self, q_input):
+        if self.dtype != q_input.dtype:
+            raise TypeError("data is %s kernel is %s" % (q_input.dtype, self.dtype))
         if self.dll is None: self._load_dll()
         kernel = self.Iqxy if q_input.is_2D else self.Iq
@@ -138 +161 @@
         ctypes and some may be pure python.
         """
-        return PyInput(q_vectors, dtype=F64)
+        return PyInput(q_vectors, dtype=self.dtype)
 
     def release(self):

sasmodels/kernelpy.py

@@ -43 +43 @@
         self.is_2D = (len(q_vectors) == 2)
         self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
-        self.q_pointers = [q.ctypes.data for q in q_vectors]
+        self.q_pointers = [q.ctypes.data for q in self.q_vectors]
 
     def release(self):