Changeset b3f6bc3 in sasmodels

Timestamp:

Feb 15, 2015 4:07:00 PM (10 years ago)

Author:

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:

f786ff3

Parents:

cb6ecf4

Message:

support pure python model Iq/Iqxy?

Location:

sasmodels

Files:

• dll.py (modified) (4 diffs)
• generate.py (modified) (5 diffs)
• gpu.py (modified) (3 diffs)
• models/spherepy.py (modified) (1 diff)
• pykernel.py (added)

sasmodels/dll.py

@@ -11 +11 @@
 
 from . import generate
+from .pykernel import PyInput, PyKernel
 
 from .generate import F32, F64
@@ -116 +117 @@
 
     def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, input)
+
         if self.dll is None: self._load_dll()
-
         kernel = self.Iqxy if input.is_2D else self.Iq
         return DllKernel(kernel, self.info, input)
@@ -125 +131 @@
         Make q input vectors available to the model.
 
-        This only needs to be done once for all models that operate on the
-        same input.  So for example, if you are adding two different models
-        together to compare to a data set, then only one model needs to
-        needs to call make_input, so long as the models have the same dtype.
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
         """
-        return DllInput(q_vectors)
+        return PyInput(q_vectors, dtype=F64)
 
     def release(self):
@@ -136 +141 @@
 
 
-class DllInput(object):
-    """
-    Make q data available to the gpu.
-
-    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
-    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
-    to get the best performance on OpenCL, which may involve shifting and
-    stretching the array to better match the memory architecture.  Additional
-    points will be evaluated with *q=1e-3*.
-
-    *dtype* is the data type for the q vectors. The data type should be
-    set to match that of the kernel, which is an attribute of
-    :class:`GpuProgram`.  Note that not all kernels support double
-    precision, so even if the program was created for double precision,
-    the *GpuProgram.dtype* may be single precision.
-
-    Call :meth:`release` when complete.  Even if not called directly, the
-    buffer will be released when the data object is freed.
-    """
-    def __init__(self, q_vectors):
-        self.nq = q_vectors[0].size
-        self.dtype = np.dtype('double')
-        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]
-
-    def release(self):
-        self.q_vectors = []
-
 class DllKernel(object):
     """
     Callable SAS kernel.
 
-    *kernel* is the DllKernel object to call.
+    *kernel* is the c function to call.
 
     *info* is the module information

sasmodels/generate.py

@@ -354 +354 @@
 """
 
+
+
+##########################################################
+#                                                        #
+#   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+#   !!                                              !!   #
+#   !!  KEEP THIS CODE CONSISTENT WITH PYKERNEL.PY  !!   #
+#   !!                                              !!   #
+#   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+#                                                        #
+##########################################################
+
+
 # Polydispersity loop body.
 # This computes the weight, and if it is sufficient, calls the scattering
@@ -387 +400 @@
 # to call the form_volume function from the user supplied kernel, and accumulate
 # a normalized weight.
-VOLUME_NORM="""const double vol_weight = %(weight)s;
-    vol += vol_weight*form_volume(%(pars)s);
+VOLUME_NORM="""const double vol_weight = %(vol_weight)s;
+    vol += vol_weight*form_volume(%(vol_pars)s);
     norm_vol += vol_weight;\
 """
@@ -490 +503 @@
     if vol_pars:
         subst = {
-            'weight': "*".join(p+"_w" for p in vol_pars),
-            'pars': ", ".join(vol_pars),
+            'vol_weight': "*".join(p+"_w" for p in vol_pars),
+            'vol_pars': ", ".join(vol_pars),
             }
         volume_norm = VOLUME_NORM%subst
@@ -637 +650 @@
             }
         source.append(WORK_FUNCTION%subst)
-    kernel_Iq = make_kernel(info, is_2D=False)
-    kernel_Iqxy = make_kernel(info, is_2D=True)
+    kernel_Iq = make_kernel(info, is_2D=False) if not callable(info['Iq']) else ""
+    kernel_Iqxy = make_kernel(info, is_2D=True) if not callable(info['Iqxy']) else ""
     kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy])
     return kernel
@@ -718 +731 @@
 def demo_time():
     import datetime
+    from .models import cylinder
+    toc = lambda: (datetime.datetime.now()-tic).total_seconds()
     tic = datetime.datetime.now()
-    toc = lambda: (datetime.datetime.now()-tic).total_seconds()
-    path = os.path.dirname("__file__")
-    doc, c = make_model(os.path.join(path, "models", "cylinder.c"))
+    source, info = make(cylinder)
     print "time:",toc()
 
 def demo():
-    from os.path import join as joinpath, dirname
-    c, info, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c"))
-    #print doc
-    #print c
+    from .models import cylinder
+    source, info = make(cylinder)
+    #print doc(cylinder)
+    print source
 
 if __name__ == "__main__":

sasmodels/gpu.py

@@ -24 +24 @@
 """
 import numpy as np
+
 import pyopencl as cl
 from pyopencl import mem_flags as mf
 
 from . import generate
+from .pykernel import PyInput, PyKernel
 
 F64_DEFS = """\
@@ -185 +187 @@
 
     def __call__(self, input):
+        # Support pure python kernel call
+        if input.is_2D and callable(self.info['Iqxy']):
+            return PyKernel(self.info['Iqxy'], self.info, input)
+        elif not input.is_2D and callable(self.info['Iq']):
+            return PyKernel(self.info['Iq'], self.info, input)
+
         if self.dtype != input.dtype:
             raise TypeError("data and kernel have different types")
@@ -202 +210 @@
         Make q input vectors available to the model.
 
-        This only needs to be done once for all models that operate on the
-        same input.  So for example, if you are adding two different models
-        together to compare to a data set, then only one model needs to
-        needs to call make_input, so long as the models have the same dtype.
+        Note that each model needs its own q vector even if the case of
+        mixture models because some models may be OpenCL, some may be
+        ctypes and some may be pure python.
         """
-        # Note: the weird interface, where make_input doesn't care which
-        # model calls it, allows us to ask the model to define the data
-        # and the caller does not need to know if it is opencl or ctypes.
-        # The returned data object is opaque.
-        return GpuInput(q_vectors, dtype=self.dtype)
+        # Support pure python kernel call
+        if len(q_vectors) == 1 and callable(self.info['Iq']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        elif callable(self.info['Iqxy']):
+            return PyInput(q_vectors, dtype=self.dtype)
+        else:
+            return GpuInput(q_vectors, dtype=self.dtype)
 
 # TODO: check that we don't need a destructor for buffers which go out of scope

sasmodels/models/spherepy.py

@@ -83 +83 @@
 
 
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
 def form_volume(radius):
-    return 1.333333333333333*pi*radius*radius*radius
+    return 1.333333333333333*pi*radius**3
 
 def Iq(q, sld, solvent_sld, radius):
+    #print "q",q
+    #print "sld,r",sld,solvent_sld,radius
     qr = q*radius
     sn, cn = sin(qr), cos(qr)
-    bes = 3 * (sn-qr*cn)/qr**3
+    # FOR VECTORIZED VERSION, UNCOMMENT THE NEXT TWO LINES
+    bes = 3 * (sn-qr*cn)/qr**3 # may be 0/0 but we fix that next line
     bes[qr==0] = 1
+    # FOR NON VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+    #bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
     fq = bes * (sld - solvent_sld) * form_volume(radius)
-    return 1.0e-4*fq*fq
+    return 1.0e-4*fq**2
+# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
+Iq.vectorized = True
 
 def Iqxy(qx, qy, sld, solvent_sld, radius):
     return Iq(sqrt(qx**2 + qy**2), sld, solvent_sld, radius)
+Iqxy.vectorized = True
 
 def ER(radius):