Index: pr_inversion/.project
===================================================================
--- pr_inversion/.project	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/.project	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,17 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<projectDescription>
+	<name>pr_inversion</name>
+	<comment></comment>
+	<projects>
+	</projects>
+	<buildSpec>
+		<buildCommand>
+			<name>org.python.pydev.PyDevBuilder</name>
+			<arguments>
+			</arguments>
+		</buildCommand>
+	</buildSpec>
+	<natures>
+		<nature>org.python.pydev.pythonNature</nature>
+	</natures>
+</projectDescription>
Index: pr_inversion/.pydevproject
===================================================================
--- pr_inversion/.pydevproject	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/.pydevproject	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,9 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<?eclipse-pydev version="1.0"?>
+
+<pydev_project>
+<pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
+<path>/pr_inversion/src</path>
+</pydev_pathproperty>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.5</pydev_property>
+</pydev_project>
Index: pr_inversion/c_extensions/Cinvertor.c
===================================================================
--- pr_inversion/c_extensions/Cinvertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/Cinvertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,424 @@
+
+#include <Python.h>
+#include "structmember.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#include "invertor.h"
+
+
+/// Error object for raised exceptions
+static PyObject * CinvertorError = NULL;
+
+#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)
+
+
+// Class definition
+typedef struct {
+    PyObject_HEAD    
+    Invertor_params params; 
+} Cinvertor;
+
+
+static void
+Cinvertor_dealloc(Cinvertor* self)
+{
+    invertor_dealloc(&(self->params));
+     
+    self->ob_type->tp_free((PyObject*)self);
+
+}
+
+static PyObject *
+Cinvertor_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
+{
+    Cinvertor *self;
+    
+    self = (Cinvertor *)type->tp_alloc(type, 0);
+   
+    return (PyObject *)self;
+}
+
+static int
+Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds)
+{
+    if (self != NULL) {    	
+    	// Create parameters
+    	invertor_init(&(self->params));
+    }
+    return 0;
+}
+
+static PyMemberDef Cinvertor_members[] = {
+    //{"params", T_OBJECT, offsetof(Cinvertor, params), 0,
+    // "Parameters"},
+    {NULL}  /* Sentinel */
+};
+
+
+/**
+ * Function to set the x data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_x(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+  
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+	self->params.x = data;
+	self->params.npoints = ndata;
+	return Py_BuildValue("i", self->params.npoints);	
+}
+
+static PyObject * get_x(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+    
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	INVECTOR(data_obj, data, ndata);
+	
+	// Check that the input array is large enough
+	if (ndata < self->params.npoints) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.get_x: input array too short for data.");
+		return NULL;		
+	}
+	
+	for(i=0; i<self->params.npoints; i++){
+		data[i] = self->params.x[i];
+	}
+	
+	return Py_BuildValue("i", self->params.npoints);	
+}
+
+/**
+ * Function to set the y data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_y(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+  
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+	self->params.y = data;
+	self->params.ny = ndata;
+	return Py_BuildValue("i", self->params.ny);	
+}
+
+static PyObject * get_y(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+    
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	INVECTOR(data_obj, data, ndata);
+	
+	// Check that the input array is large enough
+	if (ndata < self->params.ny) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.get_y: input array too short for data.");
+		return NULL;		
+	}
+	
+	for(i=0; i<self->params.ny; i++){
+		data[i] = self->params.y[i];
+	}
+	
+	return Py_BuildValue("i", self->params.npoints);	
+}
+
+/**
+ * Function to set the x data
+ * Takes an array of doubles as input
+ * Returns the number of entries found
+ */
+static PyObject * set_err(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+  
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	OUTVECTOR(data_obj,data,ndata);
+	self->params.err = data;
+	self->params.nerr = ndata;
+	return Py_BuildValue("i", self->params.nerr);	
+}
+
+static PyObject * get_err(Cinvertor *self, PyObject *args) {
+	PyObject *data_obj;
+	Py_ssize_t ndata;
+	double *data;
+    int i;
+    
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	INVECTOR(data_obj, data, ndata);
+	
+	// Check that the input array is large enough
+	if (ndata < self->params.nerr) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.get_err: input array too short for data.");
+		return NULL;		
+	}
+	
+	for(i=0; i<self->params.nerr; i++){
+		data[i] = self->params.err[i];
+	}
+	
+	return Py_BuildValue("i", self->params.npoints);	
+}
+
+/**
+ * Check the validity of the stored data
+ * Returns the number of points if it's all good, -1 otherwise
+ */
+static PyObject * is_valid(Cinvertor *self, PyObject *args) {
+	if(self->params.npoints==self->params.ny &&
+			self->params.npoints==self->params.nerr) {
+		return Py_BuildValue("i", self->params.npoints);
+	} else {
+		return Py_BuildValue("i", -1);
+	}	
+}
+
+/**
+ * Sets the maximum distance
+ */
+static PyObject * set_dmax(Cinvertor *self, PyObject *args) {
+	double d_max;
+  
+	if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL;
+	self->params.d_max = d_max;
+	return Py_BuildValue("d", self->params.d_max);	
+}
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("d", self->params.d_max);	
+}
+
+/**
+ * Gets the number of x points
+ */
+static PyObject * get_nx(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.npoints);	
+}
+
+/**
+ * Gets the number of y points
+ */
+static PyObject * get_ny(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.ny);	
+}
+
+/**
+ * Gets the number of error points
+ */
+static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
+	return Py_BuildValue("i", self->params.nerr);	
+}
+
+
+/**
+ * Function to call to evaluate the residuals
+ * @param args: input q or [q,phi]
+ * @return: function value
+ */
+static PyObject * residuals(Cinvertor *self, PyObject *args) {
+	double *pars;
+	PyObject* residuals;
+	PyObject* temp;
+	double *res;
+	int i;
+	double residual, diff;
+	// Regularization factor
+	double lmult = 0.0;
+	
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	  
+	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+	
+	OUTVECTOR(data_obj,pars,npars);
+		
+	//printf("npars: %i", npars);
+    // PyList of residuals
+	// Should create this list only once and refill it
+    residuals = PyList_New(self->params.npoints);
+    
+    for(i=0; i<self->params.npoints; i++) {
+    	//printf("\n%i: ", i);
+    	diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
+    	residual = diff*diff / (self->params.err[i]*self->params.err[i]);
+    	//printf("%i %g\n", i, residual);
+    	if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
+    	    PyErr_SetString(CinvertorError, 
+    	    	"Cinvertor.residuals: error setting residual.");
+    		return NULL;
+    	};
+    		
+    }
+    
+	return residuals;
+}
+
+/**
+ * Function to call to evaluate the scattering intensity
+ * @param args: c-parameters, and q
+ * @return: I(q)
+ */
+static PyObject * get_iq(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double q, iq_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	  
+	if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+		
+	iq_value = iq(pars, self->params.d_max, npars, q);
+	return Py_BuildValue("f", iq_value);	
+}
+
+/**
+ * Function to call to evaluate P(r)
+ * @param args: c-parameters and r
+ * @return: P(r)
+ */
+static PyObject * get_pr(Cinvertor *self, PyObject *args) {
+	double *pars;
+	double r, pr_value;
+	PyObject *data_obj;
+	Py_ssize_t npars;
+	  
+	if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
+	OUTVECTOR(data_obj,pars,npars);
+		
+	pr_value = pr(pars, self->params.d_max, npars, r);
+	return Py_BuildValue("f", pr_value);	
+}
+
+
+static PyMethodDef Cinvertor_methods[] = {
+		   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
+		   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
+		   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
+		   {"set_y", (PyCFunction)set_y, METH_VARARGS, ""},
+		   {"get_y", (PyCFunction)get_y, METH_VARARGS, ""},
+		   {"set_err", (PyCFunction)set_err, METH_VARARGS, ""},
+		   {"get_err", (PyCFunction)get_err, METH_VARARGS, ""},
+		   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
+		   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
+		   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
+		   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
+		   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, ""},
+		   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
+		   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
+		   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
+   
+   {NULL}
+};
+
+static PyTypeObject CinvertorType = {
+    PyObject_HEAD_INIT(NULL)
+    0,                         /*ob_size*/
+    "Cinvertor",             /*tp_name*/
+    sizeof(Cinvertor),             /*tp_basicsize*/
+    0,                         /*tp_itemsize*/
+    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
+    0,                         /*tp_print*/
+    0,                         /*tp_getattr*/
+    0,                         /*tp_setattr*/
+    0,                         /*tp_compare*/
+    0,                         /*tp_repr*/
+    0,                         /*tp_as_number*/
+    0,                         /*tp_as_sequence*/
+    0,                         /*tp_as_mapping*/
+    0,                         /*tp_hash */
+    0,                         /*tp_call*/
+    0,                         /*tp_str*/
+    0,                         /*tp_getattro*/
+    0,                         /*tp_setattro*/
+    0,                         /*tp_as_buffer*/
+    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
+    "Cinvertor objects",           /* tp_doc */
+    0,		               /* tp_traverse */
+    0,		               /* tp_clear */
+    0,		               /* tp_richcompare */
+    0,		               /* tp_weaklistoffset */
+    0,		               /* tp_iter */
+    0,		               /* tp_iternext */
+    Cinvertor_methods,             /* tp_methods */
+    Cinvertor_members,             /* tp_members */
+    0,                         /* tp_getset */
+    0,                         /* tp_base */
+    0,                         /* tp_dict */
+    0,                         /* tp_descr_get */
+    0,                         /* tp_descr_set */
+    0,                         /* tp_dictoffset */
+    (initproc)Cinvertor_init,      /* tp_init */
+    0,                         /* tp_alloc */
+    Cinvertor_new,                 /* tp_new */
+};
+
+
+static PyMethodDef module_methods[] = {
+    {NULL} 
+};
+
+/**
+ * Function used to add the model class to a module
+ * @param module: module to add the class to
+ */ 
+void addCinvertor(PyObject *module) {
+	PyObject *d;
+	
+    if (PyType_Ready(&CinvertorType) < 0)
+        return;
+
+    Py_INCREF(&CinvertorType);
+    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
+    
+    d = PyModule_GetDict(module);
+    CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL);
+    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
+}
+
+
+#ifndef PyMODINIT_FUNC	/* declarations for DLL import/export */
+#define PyMODINIT_FUNC void
+#endif
+PyMODINIT_FUNC
+initpr_inversion(void) 
+{
+    PyObject* m;
+
+    m = Py_InitModule3("pr_inversion", module_methods,
+                       "C extension module for inversion to P(r).");
+                       
+    addCinvertor(m);
+}
Index: pr_inversion/c_extensions/invertor.c
===================================================================
--- pr_inversion/c_extensions/invertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/invertor.c	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,86 @@
+#include <math.h>
+#include "invertor.h"
+
+double pi = 3.1416;
+
+/**
+ * Deallocate memory
+ */
+void invertor_dealloc(Invertor_params *pars) {
+	//free(pars->x);
+	//free(pars->y);
+	//free(pars->err);
+}
+
+void invertor_init(Invertor_params *pars) {
+	pars->d_max = 180;
+}
+
+
+/**
+ * P(r) of a sphere, for test purposes
+ * 
+ * @param R: radius of the sphere
+ * @param r: distance, in the same units as the radius
+ * @return: P(r)
+ */
+double pr_sphere(double R, double r) {
+    if (r <= 2.0*R) {
+        return 12.0* pow(0.5*r/R, 2.0) * pow(1.0-0.5*r/R, 2.0) * ( 2.0 + 0.5*r/R );
+    } else {
+        return 0.0;
+    }
+}
+
+/**
+ * Orthogonal functions:
+ * B(r) = 2r sin(pi*nr/d)
+ * 
+ */
+double ortho(double d_max, int n, double r) {
+	return 2.0*r*sin(pi*n*r/d_max);
+	
+}
+
+/**
+ * Fourier transform of the nth orthogonal function
+ * 
+ */
+double ortho_transformed(double d_max, int n, double q) {
+	return 8.0*pow(pi, 2.0)/q * d_max * n * pow(-1.0, n+1)
+		*sin(q*d_max) / ( pow(pi*n, 2.0) - pow(q*d_max, 2.0) );
+}
+
+/**
+ * First derivative in of the orthogonal function dB(r)/dr
+ * 
+ */
+double ortho_derived(double d_max, int n, double r) {
+    return 2.0*sin(pi*n*r/d_max) + 2.0*r*cos(pi*n*r/d_max);
+}
+
+/**
+ * Scattering intensity calculated from the expansion.
+ */
+double iq(double *pars, double d_max, int n_c, double q) {
+    double sum = 0.0;
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho_transformed(d_max, i+1, q);
+    }
+    return sum;
+}
+
+/**
+ * P(r) calculated from the expansion.
+ */
+double pr(double *pars, double d_max, int n_c, double r) {
+    double sum = 0.0; 
+	int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * ortho(d_max, i+1, r);
+    }
+    return sum;
+}
+
+
Index: pr_inversion/c_extensions/invertor.h
===================================================================
--- pr_inversion/c_extensions/invertor.h	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/c_extensions/invertor.h	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,27 @@
+#if !defined(invertor_h)
+#define invertor_h
+
+
+typedef struct {
+    double d_max;
+    double *x;
+    double *y;
+    double *err;
+    int npoints;    
+    int ny;    
+    int nerr;    
+} Invertor_params;
+
+void invertor_dealloc(Invertor_params *pars);
+
+void invertor_init(Invertor_params *pars);
+
+
+double pr_sphere(double R, double r);
+double ortho(double d_max, int n, double r);
+double ortho_transformed(double d_max, int n, double q);
+double ortho_derived(double d_max, int n, double r);
+double iq(double *c, double d_max, int n_c, double q);
+double pr(double *c, double d_max, int n_c, double r);
+
+#endif
Index: pr_inversion/invertor.py
===================================================================
--- pr_inversion/invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,72 @@
+from sans.pr.core.pr_inversion import Cinvertor
+import numpy
+
+class Invertor(Cinvertor):
+    
+    def __init__(self):
+        Cinvertor.__init__(self)
+        
+    def __setattr__(self, name, value):
+        """
+            Set the value of an attribute.
+            Access the parent class methods for
+            x, y, err and d_max.
+        """
+        if   name=='x':
+            return self.set_x(value)
+        elif name=='y':
+            return self.set_y(value)
+        elif name=='err':
+            return self.set_err(value)
+        elif name=='d_max':
+            return self.set_dmax(value)
+            
+        return Cinvertor.__setattr__(self, name, value)
+    
+    def __getattr__(self, name):
+        """
+           Return the value of an attribute
+           For the moment x, y, err and d_max are write-only
+           TODO: change that!
+        """
+        import numpy
+        if   name=='x':
+            out = numpy.ones(self.get_nx())
+            self.get_x(out)
+            return out
+        elif name=='y':
+            out = numpy.ones(self.get_ny())
+            self.get_y(out)
+            return out
+        elif name=='err':
+            out = numpy.ones(self.get_nerr())
+            self.get_err(out)
+            return out
+        elif name=='d_max':
+            return self.get_dmax()
+        elif name in self.__dict__:
+            return self.__dict__[name]
+        return None
+    
+    def invert(self, nfunc=5):
+        """
+            Perform inversion to P(r)
+        """
+        from scipy import optimize
+        
+        # First, check that the current data is valid
+        if self.is_valid()<=0:
+            raise RuntimeError, "Invertor.invert: Data array are of different length"
+        
+        p = numpy.ones(nfunc)
+        out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
+        
+        return out, cov_x
+    
+if __name__ == "__main__":
+    o = Invertor()
+
+    
+    
+    
+    
Index: pr_inversion/pr_test.py
===================================================================
--- pr_inversion/pr_test.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/pr_test.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,252 @@
+import math
+import pylab, numpy
+
+class PrInverse:
+    
+    D = 70.0
+    #alpha = 0.0001
+    alpha = 0.0
+    
+    data_x   = []
+    data_y   = []
+    data_err = []
+    
+    r = []
+    pr = []
+    pr_err = []
+    
+    xmult = []
+    
+    c = None
+    
+    def __init__(self, D=60.0):
+        ## Max distance
+        self.D = D
+        self.xmult = pylab.arange(0.001, self.D, self.D/21.0)
+
+    def l_multiplier(self, c):
+        return self.alpha * self._lmult2(c)
+        
+    def _lmult1(self, c):
+        sum = 0
+        for i in range(len(c)-1):
+            sum += (c[i+1]-c[i])**2
+        return sum
+    
+    def _lmult2(self, c):
+        sum = 0
+        for i in range(len(self.xmult)):
+            dp = 0
+            for j in range(len(c)):
+                dp += c[j] * self._ortho_derived(j+1, self.xmult[i])
+            sum += (dp**2)
+        return sum * self.D/len(self.xmult)
+    
+    def _ortho_transformed(self, n, q):
+        """
+            Fourier transform of the nth orthogonal function
+        """
+        return 8.0*math.pi**2/q * self.D * n * (-1)**(n+1) * \
+            math.sin(q*self.D) / ( (math.pi*n)**2 - (q*self.D)**2)
+            
+    def _ortho_derived(self, n, r):
+        
+        return 2.0*math.sin(math.pi*n*r/self.D) + 2.0*r*math.cos(math.pi*n*r/self.D)  
+    
+    def iq(self, c, q):
+        sum = 0.0
+        N = len(c)
+        #str =  ''
+        for i in range(N):
+            # n goes from 1 to N in the formulas
+            sum += c[i] * self._ortho_transformed(i+1, q)
+            #str += "c%i=%g  " % (i, c[i]);
+        #print str
+        #import time
+        #time.sleep(10)
+        return sum
+    
+    def residuals(self, params):
+        """
+            Calculates the vector of residuals for each point 
+            in y for a given set of input parameters.
+            @param params: list of parameter values
+            @return: vector of residuals
+        """
+        residuals = []
+        M = len(self.data_x)
+        #lmult = self.l_multiplier(params)
+        lmult = 0
+        for j in range(M):
+            if self.data_x[j]>0.0:
+                #print j
+                value = (self.data_y[j] - self.iq(params, self.data_x[j]))**2 / self.data_err[j]**2
+                #print self.data_x[j], value, lmult/float(M)
+                value += lmult/float(M)
+                residuals.append( value )
+            
+        return residuals
+    
+    def convert_res(self, res):
+        v = numpy.asarray(res)
+        if not v.flags.contiguous:
+            print "not contiguous" 
+            v = numpy.array(res)
+        
+        print v
+        res = self.c.residuals(v)
+        print res
+        #return self.c.residuals(v)
+        return res
+    
+    def new_fit(self):
+        from sans.pr.invertor import Invertor
+        from scipy import optimize
+        
+        c = Invertor()
+        c.set_dmax(self.D)
+        c.set_x(self.data_x)
+        c.set_y(self.data_y)
+        c.set_err(self.data_err)
+        
+        out, cov_x = c.invert()
+        
+        return out, cov_x
+
+    def new_fit_0(self):
+        from sans.pr.core.pr_inversion import Cinvertor
+        from scipy import optimize
+        
+        self.c = Cinvertor()
+        self.c.set_dmax(self.D)
+        self.c.set_x(self.data_x)
+        self.c.set_y(self.data_y)
+        self.c.set_err(self.data_err)
+        print "Valid", self.c.is_valid()
+        M=10
+        p = numpy.ones(M)
+        
+        p[0] = 10.0
+        
+        out, cov_x, info, mesg, success = optimize.leastsq(self.c.residuals, p, full_output=1, warning=True)
+        #out, cov_x, info, mesg, success = optimize.leastsq(self.convert_res, p, full_output=1, warning=True)
+        
+        
+        return out, cov_x
+    
+    def fit(self):
+        return self.new_fit()
+        
+        import numpy
+        from scipy import optimize
+
+
+        M = len(self.data_x)-3
+        M=3
+        print "Number of functions: ", M
+        p = numpy.zeros(M)
+        p[0] = 1000.0
+        
+        out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
+        
+        
+        return out, cov_x
+
+    def pr_theory(self, r, R):
+        """
+           
+        """
+        if r<=2*R:
+            return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R )
+        else:
+            return 0.0
+
+
+    def residuals_pr(self, params):
+        import pylab
+        
+        residuals = []
+        M = len(self.r)
+        
+        
+        for j in range(M):
+            
+            if self.r[j]>0:
+                err2 = (self.pr_err[j])**2
+                if j>0:
+                    err2 += 0.25*(self.pr[j-1]-self.pr[j-1])**2
+                value = (self.pr[j]-self.pr_calc(params, self.r[j]))**2/err2
+                residuals.append( value )
+            
+        return residuals        
+
+    def pr_calc(self, c, r):
+        sum = 0
+        N = len(c)
+        for i in range(N):
+            sum += c[i]*2.0*r*math.sin(math.pi*(i+1)*r/self.D)
+            
+        return sum
+
+    def fill_sphere(self, radius=60):
+        import pylab
+        import numpy
+        self.r = pylab.arange(0.0, radius, radius/50.0)
+        M = len(self.r)
+        self.pr = numpy.zeros(M)
+        self.pr_err = numpy.zeros(M)
+        
+        for j in range(M):
+            self.pr[j] = self.pr_theory(self.r[j], radius)
+            self.pr_err[j] = math.sqrt(self.pr[j])
+        
+
+    def fit_pr(self):
+        import numpy
+        from scipy import optimize
+        
+        
+        M = len(self.r)-3
+        M=2
+        print "Number of functions: ", M
+        p = numpy.ones(M)
+        #p = numpy.ones(15)
+        
+        out, cov_x, info, mesg, success = optimize.leastsq(self.residuals_pr, p, full_output=1, warning=True)
+      
+        return out, cov_x
+        
+        
+
+    def load(self, path = "sphere_test_data.txt"):
+        # Read the data from the data file
+        self.data_x = numpy.zeros(0)
+        self.data_y = numpy.zeros(0)
+        self.data_err = numpy.zeros(0)
+        if not path == None:
+            input_f = open(path,'r')
+            buff = input_f.read()
+            lines = buff.split('\n')
+            for line in lines:
+                try:
+                    toks = line.split()
+                    x = float(toks[0])
+                    y = float(toks[1])
+                    self.data_x = numpy.append(self.data_x, x)
+                    self.data_y = numpy.append(self.data_y, y)
+                    self.data_err = numpy.append(self.data_err, math.sqrt(y)+1000.0)
+                except:
+                    print "Error reading line: ", line
+                   
+        print "Lines read:", len(self.data_x)
+        
+
+if __name__ == "__main__": 
+    from sans.pr.core.pr_inversion import Cinvertor
+    pr = PrInverse()
+    pr.load("sphere_80.txt")
+    out, cov = pr.fit()
+    
+    print "output:"
+    print out
+    
Index: pr_inversion/setup.py
===================================================================
--- pr_inversion/setup.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/setup.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,35 @@
+import sys, os
+
+# Then build and install the modules
+from distutils.core import setup, Extension
+
+#numpy_incl = "Lib\site-packages\numpy\core\include\numpy"
+numpy_incl_path = os.path.join(sys.prefix, "Lib", "site-packages", "numpy", "core", "include", "numpy")
+print "NUMPY", numpy_incl_path
+
+# Build the module name
+srcdir  = "c_extensions"
+
+setup(
+    name="pr_inversion",
+    version = "0.1",
+    description = "Python module inversion of the scattering intensity to P(r)",
+    author = "University of Tennessee",
+    url = "danse.chem.utk.edu",
+    
+    # Use the pure python modules
+    package_dir = {"sans.pr.core":"c_extensions",
+                   "sans.pr": "."},
+    
+    packages = ["sans.pr","sans.pr.core"],
+    
+    ext_modules = [ Extension("sans.pr.core.pr_inversion",
+     sources = [
+        srcdir+"/Cinvertor.c",
+        srcdir+"/invertor.c",
+            ],
+         include_dirs=[numpy_incl_path]
+     
+     
+     )])
+        
Index: pr_inversion/sphere_80.txt
===================================================================
--- pr_inversion/sphere_80.txt	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/sphere_80.txt	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,100 @@
+<X>   <Y>
+0.0010101  2.14186e+008
+0.0020202  2.13348e+008
+0.0030303  2.11958e+008
+0.0040404  2.10025e+008
+0.00505051  2.07561e+008
+0.00606061  2.04584e+008
+0.00707071  2.01112e+008
+0.00808081  1.97169e+008
+0.00909091  1.92781e+008
+0.010101  1.87975e+008
+0.0111111  1.82783e+008
+0.0121212  1.77237e+008
+0.0131313  1.71373e+008
+0.0141414  1.65226e+008
+0.0151515  1.58833e+008
+0.0161616  1.52233e+008
+0.0171717  1.45465e+008
+0.0181818  1.38567e+008
+0.0191919  1.31578e+008
+0.020202  1.24536e+008
+0.0212121  1.17479e+008
+0.0222222  1.10442e+008
+0.0232323  1.03461e+008
+0.0242424  9.65692e+007
+0.0252525  8.97973e+007
+0.0262626  8.31748e+007
+0.0272727  7.67288e+007
+0.0282828  7.04839e+007
+0.0292929  6.4462e+007
+0.030303  5.86828e+007
+0.0313131  5.31627e+007
+0.0323232  4.79158e+007
+0.0333333  4.29533e+007
+0.0343434  3.82836e+007
+0.0353535  3.39126e+007
+0.0363636  2.98434e+007
+0.0373737  2.60768e+007
+0.0383838  2.2611e+007
+0.0393939  1.9442e+007
+0.040404  1.65639e+007
+0.0414141  1.39687e+007
+0.0424242  1.16467e+007
+0.0434343  9.58662e+006
+0.0444444  7.77602e+006
+0.0454545  6.20118e+006
+0.0464646  4.84751e+006
+0.0474747  3.69971e+006
+0.0484848  2.74197e+006
+0.0494949  1.95815e+006
+0.0505051  1.33203e+006
+0.0515152  847402
+0.0525253  488319
+0.0535354  239201
+0.0545455  84994.7
+0.0555556  11292.7
+0.0565657  4445.63
+0.0575758  51653.5
+0.0585859  141041
+0.059596  261717
+0.0606061  403816
+0.0616162  558522
+0.0626263  718084
+0.0636364  875805
+0.0646465  1.02603e+006
+0.0656566  1.1641e+006
+0.0666667  1.28634e+006
+0.0676768  1.38998e+006
+0.0686869  1.47309e+006
+0.069697  1.53454e+006
+0.0707071  1.57391e+006
+0.0717172  1.59141e+006
+0.0727273  1.5878e+006
+0.0737374  1.5643e+006
+0.0747475  1.52256e+006
+0.0757576  1.46449e+006
+0.0767677  1.39227e+006
+0.0777778  1.30823e+006
+0.0787879  1.21477e+006
+0.079798  1.11433e+006
+0.0808081  1.00933e+006
+0.0818182  902079
+0.0828283  794770
+0.0838384  689421
+0.0848485  587856
+0.0858586  491670
+0.0868687  402219
+0.0878788  320609
+0.0888889  247687
+0.089899  184050
+0.0909091  130047
+0.0919192  85790.3
+0.0929293  51176.4
+0.0939394  25900.5
+0.0949495  9481.12
+0.0959596  1283.93
+0.0969697  547.497
+0.0979798  6409.98
+0.0989899  17935.9
+0.1  34142.4
Index: pr_inversion/sphere_test_data.txt
===================================================================
--- pr_inversion/sphere_test_data.txt	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/sphere_test_data.txt	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,51 @@
+<q>   <I_1D>
+0.01  1.30173e+008
+0.02  9.59565e+007
+0.03  5.57754e+007
+0.04  2.37208e+007
+0.05  6.0257e+006
+0.06  332209
+0.07  335963
+0.08  1.03749e+006
+0.09  816241
+0.1  234630
+0.11  230.644
+0.12  100168
+0.13  182129
+0.14  111803
+0.15  16326.8
+0.16  6931.71
+0.17  45283.6
+0.18  50971.8
+0.19  19812.5
+0.2  146.532
+0.21  9676.19
+0.22  21749.2
+0.23  15602.5
+0.24  2700.12
+0.25  1044.09
+0.26  8110.76
+0.27  10054.2
+0.28  4251.65
+0.29  38.5548
+0.3  2295.45
+0.31  5499.87
+0.32  4158.67
+0.33  752.025
+0.34  308.434
+0.35  2475.65
+0.36  3174.91
+0.37  1380.1
+0.38  11.2373
+0.39  815.932
+0.4  1988.73
+0.41  1531.36
+0.42  276.75
+0.43  127.297
+0.44  1001.82
+0.45  1297.13
+0.46  566.229
+0.47  3.50245
+0.48  365.098
+0.49  887.09
+0.5  684.787
Index: pr_inversion/test/test_module.py
===================================================================
--- pr_inversion/test/test_module.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/test/test_module.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,27 @@
+from sans.pr.core.pr_inversion import Cinvertor
+import numpy
+
+c = Cinvertor()
+print c.__class__.__name__
+
+
+x = numpy.ones(3)
+print x.__class__.__name__
+print c.set_x(x)
+print c.set_y(x)
+print c.set_err(x)
+print c.is_valid()
+print c.set_dmax(60.0)
+p = numpy.ones(1)
+print c.residuals(p)
+
+
+
+v = numpy.zeros(0)
+v = numpy.append(v, 1.0)
+v = numpy.append(v, 2.0)
+print v[0], v[1]
+print len(v)
+
+
+print numpy.ndarray.
Index: pr_inversion/test/utest_invertor.py
===================================================================
--- pr_inversion/test/utest_invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
+++ pr_inversion/test/utest_invertor.py	(revision 9e8dc223ee841261f4cf577f54920873d455660d)
@@ -0,0 +1,116 @@
+"""
+    Unit tests for Invertor class
+"""
+# Disable "missing docstring" complaint
+# pylint: disable-msg=C0111
+# Disable "too many methods" complaint 
+# pylint: disable-msg=R0904 
+
+
+import unittest, math, numpy
+from sans.pr.invertor import Invertor
+        
+class TestBasicComponent(unittest.TestCase):
+    
+    def setUp(self):
+        self.invertor = Invertor()
+        self.invertor.d_max = 100.0
+        
+        # Test array
+        self.ntest = 5
+        self.x_in = numpy.ones(self.ntest)
+        for i in range(self.ntest):
+            self.x_in[i] = 1.0*i
+
+
+    def testset_dmax(self):
+        """
+            Set and read d_max
+        """
+        value = 15.0
+        self.invertor.d_max = value
+        self.assertEqual(self.invertor.d_max, value)
+        
+    def testset_x_1(self):
+        """
+            Setting and reading the x array the hard way
+        """
+        # Set x
+        self.invertor.x = self.x_in
+        
+        # Read it back
+        npts = self.invertor.get_nx()
+        x_out = numpy.ones(npts)
+        
+        self.invertor.get_x(x_out)
+
+        for i in range(self.ntest):
+            self.assertEqual(self.x_in[i], x_out[i])
+            
+    def testset_x_2(self):
+        """
+            Setting and reading the x array the easy way
+        """
+        # Set x
+        self.invertor.x = self.x_in
+        
+        # Read it back
+        x_out = self.invertor.x
+        
+        for i in range(self.ntest):
+            self.assertEqual(self.x_in[i], x_out[i])
+       
+    def testset_y(self):
+        """
+            Setting and reading the y array the easy way
+        """
+        # Set y
+        self.invertor.y = self.x_in
+        
+        # Read it back
+        y_out = self.invertor.y
+        
+        for i in range(self.ntest):
+            self.assertEqual(self.x_in[i], y_out[i])
+       
+    def testset_err(self):
+        """
+            Setting and reading the err array the easy way
+        """
+        # Set err
+        self.invertor.err = self.x_in
+        
+        # Read it back
+        err_out = self.invertor.err
+        
+        for i in range(self.ntest):
+            self.assertEqual(self.x_in[i], err_out[i])
+       
+    def test_iq(self):
+        """
+            Test iq calculation
+        """
+        q = 0.11
+        v1 = 8.0*math.pi**2/q * self.invertor.d_max *math.sin(q*self.invertor.d_max)
+        v1 /= ( math.pi**2 - (q*self.invertor.d_max)**2.0 )
+        
+        pars = numpy.ones(1)
+        self.assertAlmostEqual(self.invertor.iq(pars, q), v1, 2)
+        
+    def test_pr(self):
+        """
+            Test pr calculation
+        """
+        r = 10.0
+        v1 = 2.0*r*math.sin(math.pi*r/self.invertor.d_max)
+        pars = numpy.ones(1)
+        self.assertAlmostEqual(self.invertor.pr(pars, r), v1, 2)
+        
+    def test_getsetters(self):
+        self.invertor.new_data = 1.0
+        self.assertEqual(self.invertor.new_data, 1.0)
+        
+        self.assertEqual(self.invertor.test_no_data, None)
+        
+if __name__ == '__main__':
+    unittest.main()
