Index: DataLoader/.pydevproject
===================================================================
--- DataLoader/.pydevproject	(revision 11a0319ff1d484c300ed9c17d720b8bae9e27810)
+++ DataLoader/.pydevproject	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -5,5 +5,7 @@
 <pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
 <path>/DataLoader/src</path>
+<path>/DataLoader/extensions</path>
 </pydev_pathproperty>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.5</pydev_property>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>
 </pydev_project>
Index: DataLoader/extensions/smearer.cpp
===================================================================
--- DataLoader/extensions/smearer.cpp	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
+++ DataLoader/extensions/smearer.cpp	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -0,0 +1,164 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2009, University of Tennessee
+ */
+#include "smearer.hh"
+#include <stdio.h>
+#include <math.h>
+using namespace std;
+
+
+/**
+ * Constructor for BaseSmearer
+ *
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
+BaseSmearer :: BaseSmearer(double qmin, double qmax, int nbins) {
+	// Number of bins
+	this->nbins = nbins;
+	this->qmin = qmin;
+	this->qmax = qmax;
+	// Flag to keep track of whether we have a smearing matrix or
+	// whether we need to compute one
+	has_matrix = false;
+};
+
+/**
+ * Constructor for SlitSmearer
+ *
+ * @param width: slit width in Q units
+ * @param height: slit height in Q units
+ * @param qmin: minimum Q value
+ * @param qmax: maximum Q value
+ * @param nbins: number of Q bins
+ */
+SlitSmearer :: SlitSmearer(double width, double height, double qmin, double qmax, int nbins) :
+	BaseSmearer(qmin, qmax, nbins){
+	this->height = height;
+	this->width = width;
+};
+
+	/**
+	 * Constructor for SlitSmearer
+	 *
+	 * @param width: array slit widths for each Q point, in Q units
+	 * @param qmin: minimum Q value
+	 * @param qmax: maximum Q value
+	 * @param nbins: number of Q bins
+	 */
+QSmearer :: QSmearer(double* width, double qmin, double qmax, int nbins) :
+	BaseSmearer(qmin, qmax, nbins){
+	this->width = width;
+};
+
+/**
+ * Compute the slit smearing matrix
+ */
+void SlitSmearer :: compute_matrix(){
+
+	weights = new vector<double>(nbins*nbins,0);
+
+	// Loop over all q-values
+	for(int i=0; i<nbins; i++) {
+		double q = qmin + (double)i*(qmax-qmin)/((double)nbins-1.0);
+
+		// For each q-value, compute the weight of each other q-bin
+		// in the I(q) array
+		int npts_h = height>0 ? npts : 1;
+		int npts_w = width>0 ? npts : 1;
+
+		// If both height and width are great than zero,
+		// modify the number of points in each direction so
+		// that the total number of points is still what
+		// the user would expect (downgrade resolution)
+		if(npts_h>1 && npts_w>1){
+			npts_h = (int)ceil(sqrt((double)npts));
+			npts_w = npts_h;
+		}
+
+		double shift_h, shift_w;
+		for(int k=0; k<npts_h; k++){
+			if(npts_h==1){
+				shift_h = 0;
+			} else {
+				shift_h = height/((double)npts_h-1.0) * (double)k;
+			}
+			for(int j=0; j<npts_w; j++){
+				if(npts_w==1){
+					shift_w = 0;
+				} else {
+					shift_w = width/((double)npts_w-1.0) * (double)j;
+				}
+				double q_shifted = sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) );
+				int q_i = (int)(floor( (q_shifted-qmin) /((qmax-qmin)/((double)nbins -1.0)) ));
+
+				// Skip the entries outside our I(q) range
+				//TODO: be careful with edge effect
+				if(q_i<nbins)
+					(*weights)[i*nbins+q_i]++;
+			}
+		}
+	}
+};
+
+/**
+ * Compute the point smearing matrix
+ */
+void QSmearer :: compute_matrix(){
+	weights = new vector<double>(nbins*nbins,0);
+
+	// Loop over all q-values
+	double step = (qmax-qmin)/((double)nbins-1.0);
+	for(int i=0; i<nbins; i++) {
+		double q = qmin + (double)i*step;
+		double q_min = q - 0.5*step;
+		double q_max = q + 0.5*step;
+
+		for(int j=0; j<nbins; j++) {
+			double q_j = qmin + (double)j*step;
+
+			// Compute the fraction of the Gaussian contributing
+			// to the q bin between q_min and q_max
+			double value =  erf( (q_max-q_j)/(sqrt(2.0)*width[j]) );
+        	value -= erf( (q_min-q_j)/(sqrt(2.0)*width[j]) );
+        	(*weights)[i*nbins+j] += value;
+		}
+	}
+}
+
+/**
+ * Perform smearing by applying the smearing matrix to the input Q array
+ */
+void BaseSmearer :: smear(double *iq_in, double *iq_out, int first_bin, int last_bin){
+
+	// If we haven't computed the smearing matrix, do it now
+	if(!has_matrix) {
+		compute_matrix();
+		has_matrix = true;
+	}
+
+	// Loop over q-values and multiply apply matrix
+	for(int q_i=first_bin; q_i<=last_bin; q_i++){
+		double sum = 0.0;
+		double counts = 0.0;
+
+		for(int i=first_bin; i<=last_bin; i++){
+			sum += iq_in[i] * (*weights)[q_i*nbins+i];
+			counts += (*weights)[q_i*nbins+i];
+		}
+
+		// Normalize counts
+		iq_out[q_i] = (counts>0.0) ? sum/counts : 0;
+	}
+}
Index: DataLoader/extensions/smearer.hh
===================================================================
--- DataLoader/extensions/smearer.hh	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
+++ DataLoader/extensions/smearer.hh	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -0,0 +1,87 @@
+/**
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	If you use DANSE applications to do scientific research that leads to
+	publication, we ask that you acknowledge the use of the software with the
+	following sentence:
+
+	"This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+	copyright 2009, University of Tennessee
+ */
+#ifndef SMEARER_CLASS_H
+#define SMEARER_CLASS_H
+
+#include <vector>
+
+using namespace std;
+
+/**
+ * Base smearer class, implementing the matrix multiplication only
+ */
+class BaseSmearer {
+protected:
+	bool has_matrix;
+    // Smearing matrix
+    vector<double>* weights;
+    // Q_min (Min Q-value for I(q))
+    double qmin;
+    // Q_max (Max Q_value for I(q))
+    double qmax;
+    // Number of Q bins
+    int nbins;
+
+public:
+    // Constructor
+    BaseSmearer(double qmin, double qmax, int nbins);
+	// Smear function
+	virtual void smear(double *, double *, int, int);
+	// Compute the smearing matrix
+	virtual void compute_matrix(){};
+	// Utility function to check the number of bins
+	int get_nbins() { return nbins; }
+};
+
+
+/**
+ * Slit smearer class
+ */
+class SlitSmearer : public BaseSmearer {
+
+protected:
+    // Number of points used in the smearing computation
+    static const int npts   = 10000;
+
+public:
+    // Slit width in Q units
+    double width;
+    // Slit height in Q units
+    double height;
+
+    // Constructor
+	SlitSmearer(double width, double height, double qmin, double qmax, int nbins);
+	// Compute the smearing matrix
+	virtual void compute_matrix();
+};
+
+
+/**
+ * Point smearer class
+ */
+class QSmearer : public BaseSmearer {
+
+protected:
+    // Standard deviation in Q [A-1]
+    double* width;
+
+public:
+
+    // Constructor
+    QSmearer(double* width,double qmin, double qmax, int nbins);
+	// Compute the smearing matrix
+	virtual void compute_matrix();
+};
+
+#endif
Index: DataLoader/extensions/smearer_module.cpp
===================================================================
--- DataLoader/extensions/smearer_module.cpp	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
+++ DataLoader/extensions/smearer_module.cpp	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -0,0 +1,141 @@
+/**
+ *  Smearer module to perform point and slit smearing calculations
+
+	This software was developed by the University of Tennessee as part of the
+	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+	project funded by the US National Science Foundation.
+
+	See the license text in license.txt
+
+	copyright 2009, University of Tennessee
+ */
+#include <Python.h>
+#include <stdio.h>
+#include <smearer.hh>
+
+// Utilities
+#define INVECTOR(obj,buf,len)										\
+    do { \
+        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+#define OUTVECTOR(obj,buf,len) \
+    do { \
+        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
+        if (err < 0) return NULL; \
+        len /= sizeof(*buf); \
+    } while (0)
+
+/**
+ * Delete a SlitSmearer object
+ */
+void del_slit_smearer(void *ptr){
+	SlitSmearer* smearer = static_cast<SlitSmearer *>(ptr);
+	delete smearer;
+	return;
+}
+
+/**
+ * Create a SlitSmearer as a python object
+ */
+PyObject * new_slit_smearer(PyObject *, PyObject *args) {
+	double width;
+	double height;
+	double qmin;
+	double qmax;
+	int nbins;
+
+	if (!PyArg_ParseTuple(args, "ddddi", &width, &height, &qmin, &qmax, &nbins)) return NULL;
+	SlitSmearer* smearer = new SlitSmearer(width, height, qmin, qmax, nbins);
+	return PyCObject_FromVoidPtr(smearer, del_slit_smearer);
+}
+
+/**
+ * Delete a QSmearer object
+ */
+void del_q_smearer(void *ptr){
+	QSmearer* smearer = static_cast<QSmearer *>(ptr);
+	delete smearer;
+	return;
+}
+
+/**
+ * Create a QSmearer as a python object
+ */
+PyObject * new_q_smearer(PyObject *, PyObject *args) {
+	PyObject *width_obj;
+	Py_ssize_t n_w;
+	double* width;
+	double qmin;
+	double qmax;
+	int nbins;
+
+	if (!PyArg_ParseTuple(args, "Oddi", &width_obj, &qmin, &qmax, &nbins)) return NULL;
+	OUTVECTOR(width_obj, width, n_w);
+	QSmearer* smearer = new QSmearer(width, qmin, qmax, nbins);
+	return PyCObject_FromVoidPtr(smearer, del_q_smearer);
+}
+
+/**
+ * Smear the given input according to a given smearer object
+ */
+PyObject * smear_input(PyObject *, PyObject *args) {
+	PyObject *input_iq_obj;
+	PyObject *output_iq_obj;
+	PyObject *smear_obj;
+	Py_ssize_t n_in;
+	Py_ssize_t n_out;
+	double *input_iq;
+	double *output_iq;
+	int first_bin;
+	int last_bin;
+
+
+	if (!PyArg_ParseTuple(args, "OOOii", &smear_obj, &input_iq_obj, &output_iq_obj, &first_bin, &last_bin)) return NULL;
+	OUTVECTOR(input_iq_obj, input_iq, n_in);
+	OUTVECTOR(output_iq_obj, output_iq, n_out);
+
+	// Sanity check
+	if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+	// Set the array pointers
+	void *temp = PyCObject_AsVoidPtr(smear_obj);
+	BaseSmearer* s = static_cast<BaseSmearer *>(temp);
+
+	if(n_in!=s->get_nbins()) {
+		return Py_BuildValue("i",-2);
+	}
+
+	s->smear(input_iq, output_iq, first_bin, last_bin);
+
+	return Py_BuildValue("i",1);
+}
+
+
+/**
+ * Define module methods
+ */
+static PyMethodDef module_methods[] = {
+	{"new_slit_smearer", (PyCFunction)new_slit_smearer, METH_VARARGS,
+		  "Create a new Slit Smearer object"},
+	{"new_q_smearer", (PyCFunction)new_q_smearer, METH_VARARGS,
+		  "Create a new Q Smearer object"},
+	{"smear",(PyCFunction)smear_input, METH_VARARGS,
+		  "Smear the given input array"},
+    {NULL}
+};
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initsmearer(void)
+{
+    PyObject* m;
+
+    m = Py_InitModule3("smearer", module_methods,
+                       "Q Smearing module");
+}
Index: DataLoader/qsmearing.py
===================================================================
--- DataLoader/qsmearing.py	(revision c7ac15ecd3956c37590e94a81d320947aea7d146)
+++ DataLoader/qsmearing.py	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -6,8 +6,8 @@
 See the license text in license.txt
 
-copyright 2008, University of Tennessee
+copyright 2009, University of Tennessee
 """
 
-
+import DataLoader.extensions.smearer as smearer
 import numpy
 import math
@@ -71,6 +71,4 @@
                 _found_resolution = False
                 break
-    #print "_found_slit",_found_slit 
-    #print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0]       
     # If we found slit smearing data, return a slit smearer
     if _found_slit == True:
@@ -85,45 +83,45 @@
         self.nbins = 0
         self._weights = None
+        ## Internal flag to keep track of C++ smearer initialization
+        self._init_complete = False
+        self._smearer = None
+        
+    def __deepcopy__(self, memo={}):
+        """
+            Return a valid copy of self.
+            Avoid copying the _smearer C object and force a matrix recompute
+            when the copy is used.  
+        """
+        result = _BaseSmearer()
+        result.nbins = self.nbins
+        return result
+
         
     def _compute_matrix(self): return NotImplemented
 
-    def __call__(self, iq):
-        """
-            Return the smeared I(q) value at the given q.
-            The smeared I(q) is computed using a predetermined 
-            smearing matrix for a particular binning.
-        
-            @param q: I(q) array
-            @return: smeared I(q)
-        """
+    def __call__(self, iq_in, first_bin=0, last_bin=None):
+        """
+            Perform smearing
+        """
+        # If this is the first time we call for smearing,
+        # initialize the C++ smearer object first
+        if not self._init_complete:
+            self._initialize_smearer()
+             
+        # Get the max value for the last bin
+        if last_bin is None or last_bin>=len(iq_in):
+            last_bin = len(iq_in)-1
+        # Check that the first bin is positive
+        if first_bin<0:
+            first_bin = 0
+            
         # Sanity check
-        if len(iq) != self.nbins:
-            raise RuntimeError, "Invalid I(q) vector: inconsistent array length %s <> %s" % (len(iq), self.nbins)
-        
-        if self._weights == None:
-            self._compute_matrix()
-            
-        iq_smeared = numpy.zeros(self.nbins)
-        # Loop over q-values
-        idwb=[]
-        
-        for q_i in range(self.nbins):
-            sum = 0.0
-            counts = 0.0  
-
-            for i in range(self.nbins):
-                if iq[i]==0 or self._weights[q_i][i]==0:
-                    continue
-                else:
-                    sum += iq[i] * self._weights[q_i][i] 
-                    counts += self._weights[q_i][i]
-                    #print "i,iq[i],self._weights[q_i][i] ",i,iq[i],self._weights[q_i][i]
-            if counts == 0:
-                iq_smeared[q_i] = 0
-            else:
-                iq_smeared[q_i] = sum/counts 
-            #print "q_i,iq_smeared[q_i]",q_i,iq[i],iq_smeared[q_i]
-            #print "iq[i],iq_smeared[q_i],sum,counts,self.nbins",iq[i], iq_smeared[q_i],sum,counts,self.nbins
-        return iq_smeared    
+        if len(iq_in) != self.nbins:
+            raise RuntimeError, "Invalid I(q) vector: inconsistent array length %d != %s" % (len(iq_in), str(self.nbins))
+             
+        # Storage for smeared I(q)   
+        iq_out = numpy.zeros(self.nbins)
+        smearer.smear(self._smearer, iq_in, iq_out, first_bin, last_bin)
+        return iq_out
     
 class _SlitSmearer(_BaseSmearer):
@@ -142,4 +140,5 @@
             @param max: Q_max [A-1]
         """
+        _BaseSmearer.__init__(self)
         ## Slit width in Q units
         self.width  = width
@@ -157,47 +156,12 @@
         self._weights = None
         
-    def _compute_matrix(self):
-        """
-            Compute the smearing matrix for the current I(q) array
-        """
-        weights = numpy.zeros([self.nbins, self.nbins])
-        
-        # Loop over all q-values
-        for i in range(self.nbins):
-            q = self.min + i*(self.max-self.min)/float(self.nbins-1.0)
-            
-            # For each q-value, compute the weight of each other q-bin
-            # in the I(q) array
-            npts_h = self.npts if self.height>0 else 1 
-            npts_w = self.npts if self.width>0 else 1 
-            
-            # If both height and width are great than zero,
-            # modify the number of points in each direction so 
-            # that the total number of points is still what 
-            # the user would expect (downgrade resolution)
-            if npts_h>1 and npts_w>1:
-                npts_h = int(math.ceil(math.sqrt(self.npts)))
-                npts_w = npts_h
-                
-            for k in range(npts_h):
-                if npts_h==1:
-                    shift_h = 0
-                else:
-                    shift_h = self.height/(float(npts_h-1.0)) * k
-                for j in range(npts_w):
-                    if npts_w==1:
-                        shift_w = 0
-                    else:
-                        shift_w = self.width/(float(npts_w-1.0)) * j
-                    q_shifted = math.sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) )
-                    q_i = int(math.floor( (q_shifted-self.min)/((self.max-self.min)/(self.nbins -1.0)) ))
-                    
-                    # Skip the entries outside our I(q) range
-                    #TODO: be careful with edge effect
-                    if q_i<self.nbins:
-                        weights[i][q_i] = weights[i][q_i]+1.0
-
-        self._weights = weights
-        return self._weights
+    def _initialize_smearer(self):
+        """
+            Initialize the C++ smearer object.
+            This method HAS to be called before smearing
+        """
+        self._smearer = smearer.new_slit_smearer(self.width, self.height, self.min, self.max, self.nbins)
+        self._init_complete = True
+
 
 class SlitSmearer(_SlitSmearer):
@@ -255,4 +219,5 @@
             @param max: Q_max [A-1]
         """
+        _BaseSmearer.__init__(self)
         ## Standard deviation in Q [A-1]
         self.width  = width
@@ -266,29 +231,11 @@
         self._weights = None
         
-    def _compute_matrix(self):
-        """
-            Compute the smearing matrix for the current I(q) array
-        """
-        weights = numpy.zeros([self.nbins, self.nbins])
-        
-        # Loop over all q-values
-        step = (self.max-self.min)/float(self.nbins-1.0)
-        for i in range(self.nbins):
-            q = self.min + i*step
-            q_min = q - 0.5*step
-            q_max = q + 0.5*step
-            for j in range(self.nbins):
-                q_j = self.min + j*step
-                
-                # Compute the fraction of the Gaussian contributing
-                # to the q bin between q_min and q_max
-                #value =  math.exp(-math.pow((q_max-q_j),2)/(2*math.pow(self.width[j],2) ))
-                #value +=  math.exp(-math.pow((q_max-q_j),2)/(2*math.pow(self.width[j],2) )) 
-                value =  scipy.special.erf( (q_max-q_j)/(math.sqrt(2.0)*self.width[j]) ) 
-                value -=scipy.special.erf( (q_min-q_j)/(math.sqrt(2.0)*self.width[j]) ) 
-                weights[i][j] += value
-                                
-        self._weights = weights
-        return self._weights
+    def _initialize_smearer(self):
+        """
+            Initialize the C++ smearer object.
+            This method HAS to be called before smearing
+        """
+        self._smearer = smearer.new_q_smearer(numpy.asarray(self.width), self.min, self.max, self.nbins)
+        self._init_complete = True
         
 class QSmearer(_QSmearer):
Index: DataLoader/setup.py
===================================================================
--- DataLoader/setup.py	(revision 07f8d9a9b99c1cd5ddb838400fc2d3b2438a7a4f)
+++ DataLoader/setup.py	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -2,10 +2,10 @@
      Installation script for SANS DataLoader
 """
-import os
+import os, sys
 
 # Then build and install the modules
 from distutils.core import setup, Extension
 from distutils.sysconfig import get_python_lib
-
+numpy_incl_path = os.path.join(sys.prefix, "Lib", "site-packages", "numpy", "core", "include", "numpy")
 
 setup(
@@ -14,12 +14,20 @@
     description = "Python module for loading",
     author = "University of Tennessee",
-    #author_email = "",
     url = "http://danse.chem.utk.edu",
     
-    # Use the pure python modules
-    package_dir = {"DataLoader":"."},
+    package_dir = {"DataLoader":".",
+                   "DataLoader.extensions":"extensions"},
     
-    packages = ["DataLoader","DataLoader.readers"],
-    data_files=[(os.path.join(get_python_lib(),"DataLoader","readers"), ["readers/defaults.xml"])]
-    )
+    data_files=[(os.path.join(get_python_lib(),"DataLoader","readers"), ["readers/defaults.xml"])],
+    packages = ["DataLoader","DataLoader.readers","DataLoader.extensions"],
+    
+    ext_modules = [ Extension("DataLoader.extensions.smearer",
+     sources = [
+        "extensions/smearer_module.cpp",
+        "extensions/smearer.cpp",
+            ],
+         include_dirs=["extensions", numpy_incl_path]
+     
+     
+     )])
         
Index: DataLoader/test/utest_smearing.py
===================================================================
--- DataLoader/test/utest_smearing.py	(revision 4fe439454cb897f1e074e7c9de5112df51ee46b6)
+++ DataLoader/test/utest_smearing.py	(revision a3f8d58994d10df4c15d5350e0eb4ab2acfcb19d)
@@ -56,5 +56,5 @@
         answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ]
         for i in range(len(input)):
-            self.assertAlmostEqual(answer[i], output[i], 5)
+            self.assertAlmostEqual(answer[i], output[i], 2)
         
     def test_q(self):
@@ -68,5 +68,5 @@
         output = s(input)
         for i in range(len(input)):
-            self.assertEquals(input[i], output[i])
+            self.assertAlmostEquals(input[i], output[i], 5)
 
     def test_q2(self):
