Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 35647abd3c79d3ca6880d3b9913181dffce3a86f)
+++ sasmodels/kernelpy.py	(revision 119bd3da8a328c15484a2c3040e2d4d741bfb89c)
@@ -128,6 +128,6 @@
 
     def __call__(self, fixed, pd, cutoff=1e-5):
-        print("fixed",fixed)
-        print("pd", pd)
+        #print("fixed",fixed)
+        #print("pd", pd)
         args = self.args[:]  # grab a copy of the args
         form, form_volume = self.kernel, self.info['form_volume']
@@ -187,8 +187,4 @@
     ################################################################
 
-    #TODO: Wojtek's notes
-    #TODO: The goal is to restructure polydispersity loop
-    #so it allows fitting arbitrary polydispersity parameters
-    #and it accepts cases like coupled parameters
     weight = np.empty(len(pd), 'd')
     if weight.size > 0:
@@ -264,264 +260,2 @@
     result = scale * ret / norm + background
     return result
-=======
-"""
-Python driver for python kernels
-
-Calls the kernel with a vector of $q$ values for a single parameter set.
-Polydispersity is supported by looping over different parameter sets and
-summing the results.  The interface to :class:`PyModel` matches those for
-:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`.
-"""
-import numpy as np
-from numpy import pi, cos
-
-from .generate import F64
-
-class PyModel(object):
-    """
-    Wrapper for pure python models.
-    """
-    def __init__(self, model_info):
-        self.info = model_info
-
-    def __call__(self, q_vectors):
-        q_input = PyInput(q_vectors, dtype=F64)
-        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']
-        return PyKernel(kernel, self.info, q_input)
-
-    def release(self):
-        """
-        Free resources associated with the model.
-        """
-        pass
-
-class PyInput(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, dtype):
-        self.nq = q_vectors[0].size
-        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_pointers = [q.ctypes.data for q in self.q_vectors]
-
-    def release(self):
-        """
-        Free resources associated with the model inputs.
-        """
-        self.q_vectors = []
-
-class PyKernel(object):
-    """
-    Callable SAS kernel.
-
-    *kernel* is the DllKernel object to call.
-
-    *model_info* is the module information
-
-    *q_input* is the DllInput q vectors at which the kernel should be
-    evaluated.
-
-    The resulting call method takes the *pars*, a list of values for
-    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
-    vectors for the polydisperse parameters.  *cutoff* determines the
-    integration limits: any points with combined weight less than *cutoff*
-    will not be calculated.
-
-    Call :meth:`release` when done with the kernel instance.
-    """
-    def __init__(self, kernel, model_info, q_input):
-        self.info = model_info
-        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
-        if not getattr(kernel, 'vectorized', False):
-            if dim == '2d':
-                def vector_kernel(qx, qy, *args):
-                    """
-                    Vectorized 2D kernel.
-                    """
-                    return np.array([kernel(qxi, qyi, *args)
-                                     for qxi, qyi in zip(qx, qy)])
-            else:
-                def vector_kernel(q, *args):
-                    """
-                    Vectorized 1D kernel.
-                    """
-                    return np.array([kernel(qi, *args)
-                                     for qi in q])
-            self.kernel = vector_kernel
-        else:
-            self.kernel = kernel
-        fixed_pars = model_info['partype']['fixed-' + dim]
-        pd_pars = model_info['partype']['pd-' + dim]
-        vol_pars = model_info['partype']['volume']
-
-        # First two fixed pars are scale and background
-        pars = [p.name for p in model_info['parameters'][2:]]
-        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
-
-        # Caller needs fixed_pars and pd_pars
-        self.fixed_pars = fixed_pars
-        self.pd_pars = pd_pars
-
-    def __call__(self, fixed, pd, cutoff=1e-5):
-        #print("fixed",fixed)
-        #print("pd", pd)
-        args = self.args[:]  # grab a copy of the args
-        form, form_volume = self.kernel, self.info['form_volume']
-        # First two fixed
-        scale, background = 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,
-                     pd, self.pd_index, self.vol_index, self.theta_index)
-
-        return res
-
-    def release(self):
-        """
-        Free resources associated with the kernel.
-        """
-        self.q_input = None
-
-def _loops(form, form_volume, cutoff, scale, background,
-           args, pd, pd_index, vol_index, theta_index):
-    """
-    *form* is the name of the form function, which should be vectorized over
-    q, but otherwise have an interface like the opencl kernels, with the
-    q parameters first followed by the individual parameters in the order
-    given in model.parameters (see :mod:`sasmodels.generate`).
-
-    *form_volume* calculates the volume of the shape.  *vol_index* gives
-    the list of volume parameters
-
-    *cutoff* ignores the corners of the dispersion hypercube
-
-    *scale*, *background* multiplies the resulting form and adds an offset
-
-    *args* is the prepopulated set of arguments to the form function, starting
-    with the q vectors, and including slots for all the parameters.  The
-    values for the parameters will be substituted with values from the
-    dispersion functions.
-
-    *pd* is the list of dispersion parameters
-
-    *pd_index* are the indices of the dispersion parameters in *args*
-
-    *vol_index* are the indices of the volume parameters in *args*
-
-    *theta_index* is the index of the theta parameter for the sasview
-    spherical correction, or -1 if there is no angular dispersion
-    """
-
-    ################################################################
-    #                                                              #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #   !!                                                    !!   #
-    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
-    #   !!                                                    !!   #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #                                                              #
-    ################################################################
-
-    weight = np.empty(len(pd), 'd')
-    if weight.size > 0:
-        # weight vector, to be populated by polydispersity loops
-
-        # identify which pd parameters are volume parameters
-        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
-
-        # Sort parameters in decreasing order of pd length
-        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
-        order = np.argsort(Npd)[::-1]
-        stride = np.cumprod(Npd[order])
-        pd = [pd[index] for index in order]
-        pd_index = pd_index[order]
-        vol_weight_index = vol_weight_index[order]
-
-        fast_value = pd[0][0]
-        fast_weight = pd[0][1]
-    else:
-        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])
-    norm = np.zeros_like(ret)
-    vol = np.zeros_like(ret)
-    vol_norm = np.zeros_like(ret)
-    for k in range(stride[-1]):
-        # update polydispersity parameter values
-        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):
-                    args[pd_index[i]] = pd[i][0][index]
-                    weight[i] = pd[i][1][index]
-        else:
-            args[pd_index[0]] = fast_value[fast_index]
-            weight[0] = fast_weight[fast_index]
-
-        # Computes the weight, and if it is not sufficient then ignore this
-        # parameter set.
-        # Note: could precompute w1*...*wn so we only need to multiply by w0
-        w = np.prod(weight)
-        if w > cutoff:
-            # 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
-
-            # Call the scattering function and adds it to the total.
-            I = form(*args)
-            if np.isnan(I).any(): continue
-            ret += w * I * spherical_correction
-            norm += w
-
-            # 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 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)
-                vol_norm += vol_weight
-
-    positive = (vol * vol_norm != 0.0)
-    ret[positive] *= vol_norm[positive] / vol[positive]
-    result = scale * ret / norm + background
-    return result
Index: sasmodels/src
===================================================================
--- sasmodels/src	(revision 9c490bcf3c2afca690a8982acef8d9ab4213aff7)
+++ sasmodels/src	(revision 9c490bcf3c2afca690a8982acef8d9ab4213aff7)
@@ -0,0 +1,339 @@
+// GENERATED CODE --- DO NOT EDIT ---
+// Code is produced by sasmodels.gen from sasmodels/models/MODEL.c
+
+#ifdef __OPENCL_VERSION__
+# define USE_OPENCL
+#endif
+
+#define USE_KAHAN_SUMMATION 0
+
+// If opencl is not available, then we are compiling a C function
+// Note: if using a C++ compiler, then define kernel as extern "C"
+#ifndef USE_OPENCL
+#  ifdef __cplusplus
+      #include <cstdio>
+      #include <cmath>
+      using namespace std;
+      #if defined(_MSC_VER)
+         #include <limits>
+         #include <float.h>
+         #define kernel extern "C" __declspec( dllexport )
+         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+         inline double fmin(double x, double y) { return x>y ? y : x; }
+         inline double fmax(double x, double y) { return x<y ? y : x; }
+         inline double isnan(double x) { return _isnan(x); }
+         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
+         static double cephes_expm1(double x) {
+            // Adapted from the cephes math library.
+            // Copyright 1984 - 1992 by Stephen L. Moshier
+            if (x != x || x == 0.0) {
+               return x; // NaN and +/- 0
+            } else if (x < -0.5 || x > 0.5) {
+               return exp(x) - 1.0;
+            } else {
+               const double xsq = x*x;
+               const double p = (((
+                  +1.2617719307481059087798E-4)*xsq
+                  +3.0299440770744196129956E-2)*xsq
+                  +9.9999999999999999991025E-1);
+               const double q = ((((
+                  +3.0019850513866445504159E-6)*xsq
+                  +2.5244834034968410419224E-3)*xsq
+                  +2.2726554820815502876593E-1)*xsq
+                  +2.0000000000000000000897E0);
+               double r = x * p;
+               r =  r / (q - r);
+               return r+r;
+             }
+         }
+         #define expm1 cephes_expm1
+     #else
+         #define kernel extern "C"
+     #endif
+     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
+#  else
+     #include <stdio.h>
+     #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
+     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  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)
+#else
+#  if defined(USE_SINCOS)
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+#endif
+
+// Standard mathematical constants:
+//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4,
+//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2)
+// OpenCL defines M_constant_F for float constants, and nothing if double
+// is not enabled on the card, which is why these constants may be missing
+#ifndef M_PI
+#  define M_PI 3.141592653589793
+#endif
+#ifndef M_PI_2
+#  define M_PI_2 1.570796326794897
+#endif
+#ifndef M_PI_4
+#  define M_PI_4 0.7853981633974483
+#endif
+#ifndef M_E
+#  define M_E 2.718281828459045091
+#endif
+
+// Non-standard function library
+// pi/180, used for converting between degrees and radians
+// 4/3 pi for computing sphere volumes
+// square and cube for computing squares and cubes
+#ifndef M_PI_180
+#  define M_PI_180 0.017453292519943295
+#endif
+#ifndef M_4PI_3
+#  define M_4PI_3 4.18879020478639
+#endif
+//inline double square(double x) { return pow(x,2.0); }
+//inline double square(double x) { return pown(x,2); }
+inline double square(double x) { return x*x; }
+inline double cube(double x) { return x*x*x; }
+inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+
+
+#define VOLUME_PARAMETER_DECLARATIONS void
+#define IQ_KERNEL_NAME bessel_Iq
+#define IQ_PARAMETERS ignored
+#define IQ_FIXED_PARAMETER_DECLARATIONS const double scale, \
+    const double background, \
+    const double ignored
+#define IQ_PARAMETER_DECLARATIONS double ignored
+#define IQXY_KERNEL_NAME bessel_Iqxy
+#define IQXY_PARAMETERS ignored
+#define IQXY_FIXED_PARAMETER_DECLARATIONS const double scale, \
+    const double background, \
+    const double ignored
+#define IQXY_PARAMETER_DECLARATIONS double ignored
+
+/*
+The wrapper for gamma function from OpenCL and standard libraries
+The OpenCL gamma function fails misserably on values lower than 1.0
+while works fine on larger values.
+We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
+to function for values lower than 1.0. Namely Gamma(t) = 1/t * Gamma(t + 1)
+*/
+
+
+inline double sas_gamma( double x) { return tgamma(x+1)/x; }
+
+double form_volume(VOLUME_PARAMETER_DECLARATIONS);
+double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
+    
+    
+}
+
+
+double Iq(double q, IQ_PARAMETER_DECLARATIONS);
+double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
+    
+    return sas_gamma(q);
+    
+}
+
+
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
+    
+    // never called since no orientation or magnetic parameters.
+    //return -1.0;
+    
+}
+
+
+/*
+    ##########################################################
+    #                                                        #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #   !!                                              !!   #
+    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
+    #   !!                                              !!   #
+    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
+    #                                                        #
+    ##########################################################
+*/
+
+#ifdef IQ_KERNEL_NAME
+kernel void IQ_KERNEL_NAME(
+    global const double *q,
+    global double *result,
+    const int Nq,
+#ifdef IQ_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQ_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQ_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQ_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQ_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qi = q[i];
+#ifdef IQ_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+  #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+  #endif
+    IQ_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQ_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+      const double scattering = Iq(qi, IQ_PARAMETERS);
+      // allow kernels to exclude invalid regions by returning NaN
+      if (!isnan(scattering)) {
+        ret += weight*scattering;
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+        norm_vol += vol_weight;
+      #endif
+      }
+    //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
+    }
+    IQ_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
+
+
+#ifdef IQXY_KERNEL_NAME
+kernel void IQXY_KERNEL_NAME(
+    global const double *qx,
+    global const double *qy,
+    global double *result,
+    const int Nq,
+#ifdef IQXY_OPEN_LOOPS
+  #ifdef USE_OPENCL
+    global double *loops_g,
+  #endif
+    local double *loops,
+    const double cutoff,
+    IQXY_DISPERSION_LENGTH_DECLARATIONS,
+#endif
+    IQXY_FIXED_PARAMETER_DECLARATIONS
+    )
+{
+#ifdef USE_OPENCL
+  #ifdef IQXY_OPEN_LOOPS
+  // copy loops info to local memory
+  event_t e = async_work_group_copy(loops, loops_g, (IQXY_DISPERSION_LENGTH_SUM)*2, 0);
+  wait_group_events(1, &e);
+  #endif
+
+  int i = get_global_id(0);
+  if (i < Nq)
+#else
+  #pragma omp parallel for
+  for (int i=0; i < Nq; i++)
+#endif
+  {
+    const double qxi = qx[i];
+    const double qyi = qy[i];
+    #if USE_KAHAN_SUMMATION
+    double accumulated_error = 0.0;
+    #endif
+#ifdef IQXY_OPEN_LOOPS
+    double ret=0.0, norm=0.0;
+    #ifdef VOLUME_PARAMETERS
+    double vol=0.0, norm_vol=0.0;
+    #endif
+    IQXY_OPEN_LOOPS
+    //for (int radius_i=0; radius_i < Nradius; radius_i++) {
+    //  const double radius = loops[2*(radius_i)];
+    //  const double radius_w = loops[2*(radius_i)+1];
+
+    const double weight = IQXY_WEIGHT_PRODUCT;
+    if (weight > cutoff) {
+
+      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
+      if (!isnan(scattering)) { // if scattering is bad, exclude it from sum
+      //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
+        // equation which gave rise to this correction, leading to an
+        // attenuation of the pattern as theta moves through pi/2.  Either
+        // reverse the meanings of phi and theta in the forms, or use phi
+        // rather than theta in this correction.  Current code uses cos(theta)
+        // so that values match those of sasview.
+      #if defined(IQXY_HAS_THETA) // && 0
+        const double spherical_correction
+          = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0);
+        const double next = spherical_correction * weight * scattering;
+      #else
+        const double next = weight * scattering;
+      #endif
+      #if USE_KAHAN_SUMMATION
+        const double y = next - accumulated_error;
+        const double t = ret + y;
+        accumulated_error = (t - ret) - y;
+        ret = t;
+      #else
+        ret += next;
+      #endif
+        norm += weight;
+      #ifdef VOLUME_PARAMETERS
+        const double vol_weight = VOLUME_WEIGHT_PRODUCT;
+        vol += vol_weight*form_volume(VOLUME_PARAMETERS);
+        norm_vol += vol_weight;
+      #endif
+      }
+      //else { printf("exclude qx,qy,I:%g,%g,%g\n",qi,scattering); }
+    }
+    IQXY_CLOSE_LOOPS
+  #ifdef VOLUME_PARAMETERS
+    if (vol*norm_vol != 0.0) {
+      ret *= norm_vol/vol;
+    }
+  #endif
+    result[i] = scale*ret/norm+background;
+#else
+    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background;
+#endif
+  }
+}
+#endif
+
