Index: Models/code_capcyl.py
===================================================================
--- Models/code_capcyl.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_capcyl.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from math import asin, sqrt, fabs, atan
+import pyopencl as cl
+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
+
+class GpuCapCylinder(object):
+    PARS = {'scale':1, 'rad_cyl':1, 'rad_cap':1, 'length':1, 'sld_capcyl':1e-6, 'sld_solv':0, 'background':0,
+             'theta':0, 'phi':0}
+
+    PD_PARS = ['rad_cyl', 'length', 'rad_cap', 'theta', 'phi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        trala = open('NR_BessJ1.cpp').read()+"\n"+open('Capcyl_Kfun.cpp').read()+"\n"+open('Kernel-CapCyl.cpp').read()
+        src, qx, qy = set_precision(trala, qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+
+        #buffers
+        mf = cl.mem_flags
+        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 = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        rad_cyl,length,rad_cap,theta,phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuCapCylinder.PD_PARS]
+
+        rad_cyl.value, rad_cyl.weight = rad_cyl.get_weights(pars['rad_cyl'], 0, 10000, True)
+        rad_cap.value, rad_cap.weight = rad_cap.get_weights(pars['rad_cap'], 0, 10000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
+        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
+
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(theta.weight)
+        sub = pars['sld_capcyl']-pars['sld_solv']
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+
+        for i in xrange(len(rad_cyl.weight)):
+            for m in xrange(len(rad_cap.weight)):
+                for j in xrange(len(length.weight)):
+                    for k in xrange(len(theta.weight)):
+                        for l in xrange(len(phi.weight)):
+                            hDist = -1.0*sqrt(fabs(rad_cap.value[m]*rad_cap.value[m]-rad_cyl.value[i]*rad_cyl.value[i]))
+                            vol_i = 4.0*atan(1.0)*rad_cyl.value[i]*rad_cyl.value[i]*length.value[j]+2.0*4.0*atan(1.0)/3.0\
+                                            *((rad_cap.value[m]-hDist)*(rad_cap.value[m]-hDist)*(2*rad_cap.value[m]+hDist))
+                            self.prg.CapCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                                        real(vol_i), real(hDist), real(rad_cyl.value[i]), real(rad_cap.value[m]), real(length.value[j]),
+                                        real(theta.value[k]), real(phi.value[l]), real(sub), real(pars['scale']),
+                                        real(phi.weight[l]), real(theta.weight[k]), real(rad_cap.weight[m]),
+                                        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 vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: Models/code_coreshellcyl.py
===================================================================
--- Models/code_coreshellcyl.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_coreshellcyl.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import math
+import pyopencl as cl
+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
+
+class GpuCoreShellCylinder(object):
+    PARS = {'scale':1, 'radius':1, 'thickness':1, 'length':1, 'core_sld':1e-6, 'shell_sld':-1e-6, 'solvent_sld':0,
+            'background':0, 'axis_theta':0, 'axis_phi':0}
+    PD_PARS = ['radius', 'length', 'thickness', 'axis_phi', 'axis_theta']
+
+    def __init__(self, qx, qy, dtype='float32'):
+        #create context, queue, and build program
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-CoreShellCylinder.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+
+        #buffers
+        mf = cl.mem_flags
+        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 = np.empty_like(qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        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]
+
+        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+        thickness.value, thickness.weight = thickness.get_weights(pars['thickness'], 0, 10000, True)
+        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
+        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
+
+        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)
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        for r in xrange(len(radius.weight)):
+            for l in xrange(len(length.weight)):
+                for at in xrange(len(axis_theta.weight)):
+                    for p in xrange(len(axis_phi.weight)):
+                        for th in xrange(len(thickness.weight)):
+
+                            self.prg.CoreShellCylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                                    real(axis_theta.value[at]), real(axis_phi.value[p]), real(thickness.value[th]),
+                                    real(length.value[l]), real(radius.value[r]), real(pars['scale']),
+                                    real(radius.weight[r]), real(length.weight[l]), real(thickness.weight[th]),
+                                    real(axis_theta.weight[at]), real(axis_phi.weight[p]), real(pars['core_sld']),
+                                    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])
+                            norm_vol += radius.weight[r]*length.weight[l]*thickness.weight[th]
+                            norm += radius.weight[r]*length.weight[l]*thickness.weight[th]*axis_theta.weight[at]\
+                                    *axis_phi.weight[p]
+
+        #if size>1:
+         #   norm /= math.asin(1.0)
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
Index: Models/code_cylinder.py
===================================================================
--- Models/code_cylinder.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_cylinder.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pyopencl as cl
+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
+
+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):
+    PARS = {
+        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
+        'cyl_theta':0,'cyl_phi':0,
+    }
+    PD_PARS = ['radius', 'length', 'cyl_theta', 'cyl_phi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('NR_BessJ1.cpp').read()+"\n"+open('Kernel-Cylinder.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        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 = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        radius, length, cyl_theta, cyl_phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuCylinder.PD_PARS]
+
+        #Get the weights for each
+        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)
+
+        #Perform the computation, with all weight points
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(cyl_theta.weight)
+        sub = pars['sldCyl'] - pars['sldSolv']
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        #Loop over radius, length, theta, phi weight points
+        for i in xrange(len(radius.weight)):
+            for j in xrange(len(length.weight)):
+                for k in xrange(len(cyl_theta.weight)):
+                    for l in xrange(len(cyl_phi.weight)):
+                        self.prg.CylinderKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(sub),
+                                           real(radius.value[i]), real(length.value[j]), real(pars['scale']),
+                                           real(radius.weight[i]), real(length.weight[j]), real(cyl_theta.weight[k]),
+                                           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]
+                        norm += radius.weight[i]*length.weight[j]*cyl_theta.weight[k]*cyl_phi.weight[l]
+
+       # if size > 1:
+        #    norm /= math.asin(1.0)
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm+pars['background']
+
+class OneDGpuCylinder(object):
+    PARS = {
+        'scale':1,'radius':1,'length':1,'sldCyl':1e-6,'sldSolv':0,'background':0,
+        'bolim':0, 'uplim':90
+    }
+    PD_PARS = ['radius', 'length']
+
+    def __init__(self, q, dtype='float32'):
+
+        #create context, queue, and build program
+        ctx,_queue = card()
+        trala = open('NR_BessJ1.cpp').read()+"\n"+open('OneDCyl_Kfun.cpp').read()+"\n"+open('Kernel-OneDCylinder.cpp').read()
+        src, self.q = set_precision_1d(trala, q, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+
+        #buffers
+        mf = cl.mem_flags
+        self.q_b = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.q)
+        self.res_b = cl.Buffer(ctx, mf.WRITE_ONLY, q.nbytes)
+        self.res = np.empty_like(self.q)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        radius, length = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in OneDGpuCylinder.PD_PARS]
+
+        #Get the weights for each
+        radius.value, radius.weight = radius.get_weights(pars['radius'], 0, 10000, True)
+        length.value, length.weight = length.get_weights(pars['length'], 0, 10000, True)
+
+        #Perform the computation, with all weight points
+        sum, norm, vol = 0.0, 0.0, 0.0,
+        sub = pars['sldCyl'] - pars['sldSolv']
+
+        real = np.float32 if self.q.dtype == np.dtype('float32') else np.float64
+        #Loop over radius, length, theta, phi weight points
+        for r in xrange(len(radius.weight)):
+            for l in xrange(len(length.weight)):
+                        self.prg.OneDCylKernel(queue, self.q.shape, None, self.q_b, self.res_b, real(sub),
+                                           real(length.value[l]), real(radius.value[r]), real(pars['scale']),
+                                           np.uint32(self.q.size), real(pars['uplim']), real(pars['bolim']))
+                        cl.enqueue_copy(queue, self.res, self.res_b)
+                        sum += radius.weight[r]*length.weight[l]*self.res*pow(radius.value[r],2)*length.value[l]
+                        vol += radius.weight[r]*length.weight[l] *pow(radius.value[r],2)*length.value[l]
+                        norm += radius.weight[r]*length.weight[l]
+
+        if vol != 0.0 and norm != 0.0:
+            sum *= norm/vol
+
+        return sum/norm + pars['background']
+
+       
+
+        
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: Models/code_ellipse.py
===================================================================
--- Models/code_ellipse.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_ellipse.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import math
+import pyopencl as cl
+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
+
+class GpuEllipse(object):
+    PARS = {
+    'scale':1, 'radius_a':1, 'radius_b':1, 'sldEll':1e-6, 'sldSolv':0, 'background':0, 'axis_theta':0, 'axis_phi':0,
+    }
+    PD_PARS = ['radius_a', 'radius_b', 'axis_theta', 'axis_phi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel-Ellipse.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        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 = np.empty_like(self.qx)
+
+    def eval(self, pars):
+    #b_n = radius_b # want, a_n = radius_a # want, etc
+        _ctx,queue = card()
+        radius_a, radius_b, axis_theta, axis_phi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuEllipse.PD_PARS]
+
+        radius_a.value, radius_a.weight = radius_a.get_weights(pars['radius_a'], 0, 10000, True)
+        radius_b.value, radius_b.weight = radius_b.get_weights(pars['radius_b'], 0, 10000, True)
+        axis_theta.value, axis_theta.weight = axis_theta.get_weights(pars['axis_theta'], -90, 180, False)
+        axis_phi.value, axis_phi.weight = axis_phi.get_weights(pars['axis_phi'], -90, 180, False)
+
+        #Perform the computation, with all weight points
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(axis_theta.weight)
+        sub = pars['sldEll'] - pars['sldSolv']
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+
+        #Loop over radius weight points
+        for i in xrange(len(radius_a.weight)):
+            #Loop over length weight points
+            for j in xrange(len(radius_b.weight)):
+                #Average over theta distribution
+                for k in xrange(len(axis_theta.weight)):
+                    #Average over phi distribution
+                    for l in xrange(len(axis_phi.weight)):
+                        #call the kernel
+                        self.prg.EllipsoidKernel(queue, self.qx.shape, None, real(radius_a.weight[i]),
+                                        real(radius_b.weight[j]), real(axis_theta.weight[k]),
+                                        real(axis_phi.weight[l]), real(pars['scale']), real(radius_a.value[i]),
+                                        real(radius_b.value[j]), real(sub), real(axis_theta.value[k]),
+                                        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]
+                        norm += radius_a.weight[i]*radius_b.weight[j]*axis_theta.weight[k]*axis_phi.weight[l]
+        # Averaging in theta needs an extra normalization
+        # factor to account for the sin(theta) term in the
+        # integration (see documentation).
+
+    #    if size > 1:
+     #       norm /= math.asin(1.0)
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm+pars['background']
+
+
+def demo():
+    from time import time
+    import matplotlib.pyplot as plt
+
+    #create qx and qy evenly spaces
+    qx = np.linspace(-.02, .02, 128)
+    qy = np.linspace(-.02, .02, 128)
+    qx, qy = np.meshgrid(qx, qy)
+
+    #saved shape of qx
+    r_shape = qx.shape
+    #reshape for calculation; resize as float32
+    qx = qx.flatten()
+    qy = qy.flatten()
+
+    #int main
+    pars = EllipsoidParameters(.027, 60, 180, .297e-6, 5.773e-06, 4.9, 0, 90)
+
+    t = time()
+    result = GpuEllipse(qx, qy)
+    result.x = result.ellipsoid_fit(qx, qy, pars, b_n=35, t_n=35, a_n=1, p_n=1, sigma=3, b_w=.1, t_w=.1, a_w=.1, p_w=.1)
+    result.x = np.reshape(result.x, r_shape)
+    tt = time()
+    print("Time taken: %f" % (tt - t))
+
+    plt.pcolormesh(result.x)
+    plt.show()
+
+
+if __name__ == "__main__":
+    demo()
+
+
+
+
Index: Models/code_lamellar.py
===================================================================
--- Models/code_lamellar.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_lamellar.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pyopencl as cl
+from 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
+
+
+class GpuLamellar(object):
+    PARS = {
+        'scale':1, 'bi_thick':1, 'sld_bi':1e-6, 'sld_sol':0, 'background':0,
+    }
+    PD_PARS = {'bi_thick'}
+    def __init__(self, qx, qy, dtype='float32'):
+
+        #create context, queue, and build program
+        self.ctx = cl.create_some_context()
+        self.queue = cl.CommandQueue(self.ctx)
+        src,qx,qy = set_precision(open('Kernel-Lamellar.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(self.ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        self.qx_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qx)
+        self.qy_b = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.qy)
+        self.res_b = cl.Buffer(self.ctx, mf.WRITE_ONLY, qx.nbytes)
+        self.res = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        bi_thick = GaussianDispersion(int(pars['bi_thick_pd_n']), pars['bi_thick_pd'], pars['bi_thick_pd_nsigma'])
+        bi_thick.value, bi_thick.weight = bi_thick.get_weights(pars['bi_thick'], 0, 10000, True)
+
+        sum, norm = 0.0, 0.0
+        sub = pars['sld_bi'] - pars['sld_sol']
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        for i in xrange(len(bi_thick.weight)):
+            self.prg.LamellarKernel(self.queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b, real(bi_thick.value[i]),
+                                    real(pars['scale']), real(sub), np.uint32(self.qx.size))
+            cl.enqueue_copy(self.queue, self.res, self.res_b)
+
+            sum += bi_thick.weight[i]*self.res
+            norm += bi_thick.weight[i]
+
+        return sum/norm + pars['background']
+
+def demo():
+
+    from time import time
+    import matplotlib.pyplot as plt
+
+    #create qx and qy evenly spaces
+    qx = np.linspace(-.01, .01, 128)
+    qy = np.linspace(-.01, .01, 128)
+    qx, qy = np.meshgrid(qx, qy)
+
+    #saved shape of qx
+    r_shape = qx.shape
+
+    #reshape for calculation; resize as float32
+    qx = qx.flatten()
+    qy = qy.flatten()
+
+    pars = LamellarParameters(scale=1, bi_thick=100, sld_bi=.291e-6, sld_sol=5.77e-6, background=0)
+    t = time()
+    result = GpuLamellar(qx, qy)
+    result.x = result.lamellar_fit(pars, b_n=35, b_w=.1, sigma=3)
+    result.x = np.reshape(result.x, r_shape)
+    tt = time()
+
+    print("Time taken: %f" % (tt - t))
+
+    f = open("r.txt", "w")
+    for x in xrange(len(r_shape)):
+        f.write(str(result.x[x]))
+        f.write("\n")
+
+    plt.pcolormesh(np.log10(result.x))
+    plt.show()
+
+
+if __name__ == "__main__":
+    demo()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: Models/code_triaxialellipse.py
===================================================================
--- Models/code_triaxialellipse.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/code_triaxialellipse.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+from math import asin
+import pyopencl as cl
+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
+
+class GpuTriEllipse:
+    PARS = {'scale':1, 'axisA':35, 'axisB':100, 'axisC':400, 'sldEll':1e-6, 'sldSolv':6.3e-6, 'background':0,
+            'theta':0, 'phi':0, 'psi':0}
+
+    PD_PARS = ['axisA', 'axisB', 'axisC', 'theta', 'phi', 'psi']
+
+    def __init__(self, qx, qy, dtype='float32'):
+        ctx,_queue = card()
+        src, qx, qy = set_precision(open('Kernel-TriaxialEllipse.cpp').read(), qx, qy, dtype=dtype)
+        self.prg = cl.Program(ctx, src).build()
+        self.qx, self.qy = qx, qy
+
+        #buffers
+        mf = cl.mem_flags
+        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 = np.empty_like(self.qx)
+
+    def eval(self, pars):
+
+        _ctx,queue = card()
+        axisA, axisB, axisC, theta, phi, psi = \
+            [GaussianDispersion(int(pars[base+'_pd_n']), pars[base+'_pd'], pars[base+'_pd_nsigma'])
+             for base in GpuTriEllipse.PD_PARS]
+
+        axisA.value, axisA.weight = axisA.get_weights(pars['axisA'], 0, 10000, True)
+        axisB.value, axisB.weight = axisB.get_weights(pars['axisB'], 0, 10000, True)
+        axisC.value, axisC.weight = axisC.get_weights(pars['axisC'], 0, 10000, True)
+        theta.value, theta.weight = theta.get_weights(pars['theta'], -90, 180, False)
+        phi.value, phi.weight = phi.get_weights(pars['phi'], -90, 180, False)
+        psi.value, psi.weight = psi.get_weights(pars['psi'], -90, 180, False)
+
+        sum, norm, norm_vol, vol = 0.0, 0.0, 0.0, 0.0
+        size = len(theta.weight)
+        sub = pars['sldEll'] - pars['sldSolv']
+
+        real = np.float32 if self.qx.dtype == np.dtype('float32') else np.float64
+        for a in xrange(len(axisA.weight)):
+            for b in xrange(len(axisB.weight)):
+                for c in xrange(len(axisC.weight)):
+                    for t in xrange(len(theta.weight)):
+                        for i in xrange(len(phi.weight)):
+                            for s in xrange(len(psi.weight)):
+                                self.prg.TriaxialEllipseKernel(queue, self.qx.shape, None, self.qx_b, self.qy_b, self.res_b,
+                                            real(sub), real(pars['scale']), real(axisA.value[a]), real(axisB.value[b]),
+                                            real(axisC.value[c]), real(phi.value[i]), real(theta.value[t]), real(psi.value[s]),
+                                            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]
+                                norm_vol += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]
+                                norm += axisA.weight[a]*axisB.weight[b]*axisC.weight[c]*theta.weight[t]*phi.weight[i]*psi.weight[s]
+
+      #  if size > 1:
+       #     norm /= asin(1.0)
+
+        if vol != 0.0 and norm_vol != 0.0:
+            sum *= norm_vol/vol
+
+        return sum/norm + pars['background']
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: Models/weights.py
===================================================================
--- Models/weights.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
+++ Models/weights.py	(revision dbb0048466f8b528c1953a0eb5a496b6c1c25e32)
@@ -0,0 +1,26 @@
+import numpy as np
+
+class GaussianDispersion(object):
+    def __init__(self, npts=35, width=0, nsigmas=3): #number want, percent deviation, #standard deviations from mean
+        self.type = 'gaussian'
+        self.npts = npts
+        self.width = width
+        self.nsigmas = nsigmas
+
+    def get_pars(self):
+        return self.__dict__
+
+    def get_weights(self, center, min, max, relative):
+        """ *center* is the center of the distribution
+        *min*,*max* are the min, max allowed values
+        *relative* is True if the width is relative to the center instead of absolute
+        For polydispersity use relative.  For orientation parameters use absolute."""
+        npts, width, nsigmas = self.npts, self.width, self.nsigmas
+        sigma = width * center if relative else width
+        if sigma == 0:
+            return np.array([center, 1.], 'd')
+        x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts)
+        x = x[(x >= min) & (x <= max)]
+        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
+        return x, px
+
