Changeset c9e31e2 in sasmodels for sasmodels/kernelpy.py

Timestamp:

Mar 6, 2015 10:54:22 AM (9 years ago)

Author:

richardh

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:

9f6f2f8, ab87a12

Parents:

d60b433 (diff), 3c56da87 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

File:

• sasmodels/kernelpy.py (modified) (12 diffs)

sasmodels/kernelpy.py

@@ -1 +1 @@
 import numpy as np
-from numpy import pi, sin, cos, sqrt
-
-from .generate import F32, F64
+from numpy import pi, cos
+
+from .generate import F64
+
+class PyModel(object):
+    def __init__(self, info):
+        self.info = info
+
+    def __call__(self, q_input):
+        kernel = self.info['Iqxy'] if q_input.is_2D else self.info['Iq']
+        return PyKernel(kernel, self.info, q_input)
+
+    # pylint: disable=no-self-use
+    def make_input(self, q_vectors):
+        return PyInput(q_vectors, dtype=F64)
+
+    def release(self):
+        pass
 
 class PyInput(object):
@@ -27 +42 @@
         self.dtype = dtype
         self.is_2D = (len(q_vectors) == 2)
-        self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors]
+        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
         self.q_pointers = [q.ctypes.data for q in q_vectors]
 
@@ -41 +56 @@
     *info* is the module information
 
-    *input* is the DllInput q vectors at which the kernel should be
+    *q_input* is the DllInput q vectors at which the kernel should be
     evaluated.
 
@@ -52 +67 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, info, input):
+    def __init__(self, kernel, info, q_input):
         self.info = info
-        self.input = input
-        self.res = np.empty(input.nq, input.dtype)
-        dim = '2d' if input.is_2D else '1d'
+        self.q_input = q_input
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2D else '1d'
         # Loop over q unless user promises that the kernel is vectorized by
         # taggining it with vectorized=True
@@ -62 +77 @@
             if dim == '2d':
                 def vector_kernel(qx, qy, *args):
-                    return np.array([kernel(qxi,qyi,*args) for qxi,qyi in zip(qx,qy)])
+                    return np.array([kernel(qxi, qyi, *args)
+                                     for qxi, qyi in zip(qx, qy)])
             else:
                 def vector_kernel(q, *args):
-                    return np.array([kernel(qi,*args) for qi in q])
+                    return np.array([kernel(qi, *args)
+                                     for qi in q])
             self.kernel = vector_kernel
         else:
             self.kernel = kernel
-        fixed_pars = info['partype']['fixed-'+dim]
-        pd_pars = info['partype']['pd-'+dim]
+        fixed_pars = info['partype']['fixed-' + dim]
+        pd_pars = info['partype']['pd-' + dim]
         vol_pars = info['partype']['volume']
 
         # First two fixed pars are scale and background
         pars = [p[0] for p in info['parameters'][2:]]
-        offset = len(self.input.q_vectors)
-        self.args = self.input.q_vectors + [None]*len(pars)
-        self.fixed_index = np.array([pars.index(p)+offset for p in fixed_pars[2:] ])
-        self.pd_index = np.array([pars.index(p)+offset for p in pd_pars])
-        self.vol_index = np.array([pars.index(p)+offset for p in vol_pars])
-        try: self.theta_index = pars.index('theta')+offset
+        offset = len(self.q_input.q_vectors)
+        self.args = self.q_input.q_vectors + [None] * len(pars)
+        self.fixed_index = np.array([pars.index(p) + offset
+                                     for p in fixed_pars[2:]])
+        self.pd_index = np.array([pars.index(p) + offset
+                                  for p in pd_pars])
+        self.vol_index = np.array([pars.index(p) + offset
+                                   for p in vol_pars])
+        try: self.theta_index = pars.index('theta') + offset
         except ValueError: self.theta_index = -1
 
@@ -94 +114 @@
         # First two fixed
         scale, background = fixed[:2]
-        for index,value in zip(self.fixed_index, fixed[2:]):
+        for index, value in zip(self.fixed_index, fixed[2:]):
             args[index] = float(value)
-        res = _loops(form, form_volume, cutoff, scale, background,  args,
+        res = _loops(form, form_volume, cutoff, scale, background, args,
                      pd, self.pd_index, self.vol_index, self.theta_index)
 
@@ -102 +122 @@
 
     def release(self):
-        self.input = None
+        self.q_input = None
 
 def _loops(form, form_volume, cutoff, scale, background,
@@ -134 +154 @@
     """
 
-    ##########################################################
-    #                                                        #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #   !!                                              !!   #
-    #   !!  KEEP THIS CODE CONSISTENT WITH GENERATE.PY  !!   #
-    #   !!                                              !!   #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #                                                        #
-    ##########################################################
+    ################################################################
+    #                                                              #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                                    !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
+    #   !!                                                    !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                              #
+    ################################################################
 
     weight = np.empty(len(pd), 'd')
@@ -164 +184 @@
         stride = np.array([1])
         vol_weight_index = slice(None, None)
+        # keep lint happy
+        fast_value = [None]
+        fast_weight = [None]
 
     ret = np.zeros_like(args[0])
@@ -171 +194 @@
     for k in range(stride[-1]):
         # update polydispersity parameter values
-        fast_index = k%stride[0]
+        fast_index = k % stride[0]
         if fast_index == 0:  # bottom loop complete ... check all other loops
             if weight.size > 0:
-                for i,index, in enumerate(k%stride):
+                for i, index, in enumerate(k % stride):
                     args[pd_index[i]] = pd[i][0][index]
                     weight[i] = pd[i][1][index]
@@ -180 +203 @@
             args[pd_index[0]] = fast_value[fast_index]
             weight[0] = fast_weight[fast_index]
-        # This computes the weight, and if it is sufficient, calls the scattering
-        # function and adds it to the total.  If there is a volume normalization,
-        # it will also be added here.
+        # This computes the weight, and if it is sufficient, calls the
+        # scattering function and adds it to the total.  If there is a volume
+        # normalization, it will also be added here.
         # Note: make sure this is consistent with the code in PY_LOOP_BODY!!
         # Note: can precompute w1*w2*...*wn
@@ -188 +211 @@
         if w > cutoff:
             I = form(*args)
-            positive = (I>=0.0)
-
-            # Note: can precompute spherical correction if theta_index is not the fast index
-            # Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
-            #spherical_correction = abs(sin(pi*args[theta_index])) if theta_index>=0 else 1.0
-            #spherical_correction = abs(cos(pi*args[theta_index]))*pi/2 if theta_index>=0 else 1.0
-            spherical_correction = 1.0
-            ret += w*I*spherical_correction*positive
-            norm += w*positive
+            positive = (I >= 0.0)
+
+            # Note: can precompute spherical correction if theta_index is not
+            # the fast index. Correction factor for spherical integration
+            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
+            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
+                                    if theta_index >= 0 else 1.0)
+            #spherical_correction = 1.0
+            ret += w * I * spherical_correction * positive
+            norm += w * positive
 
             # Volume normalization.
-            # If there are "volume" polydispersity parameters, then these will be used
-            # to call the form_volume function from the user supplied kernel, and accumulate
-            # a normalized weight.
-            # Note: can precompute volume norm if the fast index is not a volume index
+            # If there are "volume" polydispersity parameters, then these
+            # will be used to call the form_volume function from the user
+            # supplied kernel, and accumulate a normalized weight.
+            # Note: can precompute volume norm if fast index is not a volume
             if form_volume:
                 vol_args = [args[index] for index in vol_index]
                 vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight*form_volume(*vol_args)*positive
-                vol_norm += vol_weight*positive
-
-    positive = (vol*vol_norm != 0.0)
+                vol += vol_weight * form_volume(*vol_args) * positive
+                vol_norm += vol_weight * positive
+
+    positive = (vol * vol_norm != 0.0)
     ret[positive] *= vol_norm[positive] / vol[positive]
-    result = scale*ret/norm+background
+    result = scale * ret / norm + background
     return result