Index: pr_inversion/c_extensions/Cinvertor.c
===================================================================
--- pr_inversion/c_extensions/Cinvertor.c	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/c_extensions/Cinvertor.c	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
@@ -80,8 +80,23 @@
 	Py_ssize_t ndata;
 	double *data;
+	int i;
   
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
 	OUTVECTOR(data_obj,data,ndata);
-	self->params.x = data;
+	
+	free(self->params.x);
+	self->params.x = (double*) malloc(ndata*sizeof(double));
+	
+	if(self->params.x==NULL) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.set_x: problem allocating memory.");
+		return NULL;		
+	}
+	
+	for (i=0; i<ndata; i++) {
+		self->params.x[i] = data[i];
+	}
+	
+	//self->params.x = data;
 	self->params.npoints = ndata;
 	return Py_BuildValue("i", self->params.npoints);	
@@ -95,5 +110,5 @@
     
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
-	INVECTOR(data_obj, data, ndata);
+	OUTVECTOR(data_obj, data, ndata);
 	
 	// Check that the input array is large enough
@@ -120,8 +135,23 @@
 	Py_ssize_t ndata;
 	double *data;
+	int i;
   
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
 	OUTVECTOR(data_obj,data,ndata);
-	self->params.y = data;
+	
+	free(self->params.y);
+	self->params.y = (double*) malloc(ndata*sizeof(double));
+	
+	if(self->params.y==NULL) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.set_y: problem allocating memory.");
+		return NULL;		
+	}
+	
+	for (i=0; i<ndata; i++) {
+		self->params.y[i] = data[i];
+	}	
+	
+	//self->params.y = data;
 	self->params.ny = ndata;
 	return Py_BuildValue("i", self->params.ny);	
@@ -135,5 +165,5 @@
     
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
-	INVECTOR(data_obj, data, ndata);
+	OUTVECTOR(data_obj, data, ndata);
 	
 	// Check that the input array is large enough
@@ -160,8 +190,23 @@
 	Py_ssize_t ndata;
 	double *data;
+	int i;
   
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
 	OUTVECTOR(data_obj,data,ndata);
-	self->params.err = data;
+	
+	free(self->params.err);
+	self->params.err = (double*) malloc(ndata*sizeof(double));
+	
+	if(self->params.err==NULL) {
+	    PyErr_SetString(CinvertorError, 
+	    	"Cinvertor.set_err: problem allocating memory.");
+		return NULL;		
+	}
+	
+	for (i=0; i<ndata; i++) {
+		self->params.err[i] = data[i];
+	}
+	
+	//self->params.err = data;
 	self->params.nerr = ndata;
 	return Py_BuildValue("i", self->params.nerr);	
@@ -175,5 +220,5 @@
     
 	if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
-	INVECTOR(data_obj, data, ndata);
+	OUTVECTOR(data_obj, data, ndata);
 	
 	// Check that the input array is large enough
@@ -414,4 +459,12 @@
 }
 
+static PyObject * basefunc_ft(Cinvertor *self, PyObject *args) {
+	double d_max, q;
+	int n;
+	
+	if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
+	return Py_BuildValue("f", ortho_transformed(d_max, n, q));	
+	
+}
 
 static PyMethodDef Cinvertor_methods[] = {
@@ -435,4 +488,5 @@
 		   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
 		   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
+		   {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, ""},
    
    {NULL}
Index: pr_inversion/c_extensions/invertor.c
===================================================================
--- pr_inversion/c_extensions/invertor.c	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/c_extensions/invertor.c	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
@@ -8,7 +8,7 @@
  */
 void invertor_dealloc(Invertor_params *pars) {
-	//free(pars->x);
-	//free(pars->y);
-	//free(pars->err);
+	free(pars->x);
+	free(pars->y);
+	free(pars->err);
 }
 
@@ -133,3 +133,2 @@
 }
 
-
Index: pr_inversion/c_extensions/invertor.h
===================================================================
--- pr_inversion/c_extensions/invertor.h	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/c_extensions/invertor.h	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
@@ -29,4 +29,3 @@
 void pr_err(double *pars, double *err, double d_max, int n_c, 
 		double r, double *pr_value, double *pr_value_err);
-
 #endif
Index: pr_inversion/invertor.py
===================================================================
--- pr_inversion/invertor.py	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/invertor.py	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
@@ -5,5 +5,7 @@
     
     ## Chisqr of the last computation
-    chisqr = 0
+    chi2  = 0
+    ## Time elapsed for last computation
+    elapsed = 0
     
     def __init__(self):
@@ -58,4 +60,20 @@
         return None
     
+    def clone(self):
+        """
+            Return a clone of this instance
+        """
+        invertor = Invertor()
+        invertor.chi2    = self.chi2 
+        invertor.elapsed = self.elapsed 
+        invertor.alpha   = self.alpha
+        invertor.d_max   = self.d_max
+        
+        invertor.x = self.x
+        invertor.y = self.y
+        invertor.err = self.err
+        
+        return invertor
+    
     def invert(self, nfunc=5):
         """
@@ -63,4 +81,5 @@
         """
         from scipy import optimize
+        import time
         
         # First, check that the current data is valid
@@ -69,4 +88,5 @@
         
         p = numpy.ones(nfunc)
+        t_0 = time.time()
         out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
         
@@ -78,4 +98,7 @@
         
         self.chi2 = chisqr
+
+        # Store computation time
+        self.elapsed = time.time() - t_0
         
         return out, cov_x
@@ -92,4 +115,5 @@
         
         p = numpy.ones(nfunc)
+        t_0 = time.time()
         out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True)
         
@@ -102,4 +126,7 @@
         self.chisqr = chisqr
         
+        # Store computation time
+        self.elapsed = time.time() - t_0
+
         return out, cov_x
     
@@ -117,4 +144,105 @@
 
         return self.get_pr_err(c, c_err, r)
+       
+    def lstsq(self, nfunc=5):
+        import math
+        from scipy.linalg.basic import lstsq
+        
+        # a -- An M x N matrix.
+        # b -- An M x nrhs matrix or M vector.
+        npts = len(self.x)
+        nq   = 20
+        sqrt_alpha = math.sqrt(self.alpha)
+        
+        a = numpy.zeros([npts+nq, nfunc])
+        b = numpy.zeros(npts+nq)
+        err = numpy.zeros(nfunc)
+        
+        for j in range(nfunc):
+            for i in range(npts):
+                a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
+            for i_q in range(nq):
+                r = self.d_max/nq*i_q
+                #a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*math.fabs(math.sin(math.pi*(j+1)*r/self.d_max) + math.pi*(j+1)*r/self.d_max * math.cos(math.pi*(j+1)*r/self.d_max))     
+                a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*(2.0*math.pi*(j+1)/self.d_max*math.cos(math.pi*(j+1)*r/self.d_max) + math.pi**2*(j+1)**2*r/self.d_max**2 * math.sin(math.pi*(j+1)*r/self.d_max))     
+        
+        for i in range(npts):
+            b[i] = self.y[i]/self.err[i]
+            
+        c, chi2, rank, n = lstsq(a, b)
+        self.chi2 = chi2
+                
+        at = numpy.transpose(a)
+        inv_cov = numpy.zeros([nfunc,nfunc])
+        for i in range(nfunc):
+            for j in range(nfunc):
+                inv_cov[i][j] = 0.0
+                for k in range(npts):
+                    inv_cov[i][j] = at[i][k]*a[k][j]
+                    
+        # Compute the reg term size for the output
+        sum_sig = 0.0
+        sum_reg = 0.0
+        for j in range(nfunc):
+            for i in range(npts):
+                sum_sig += (a[i][j])**2
+            for i in range(nq):
+                sum_reg += (a[i_q+npts][j])**2
+                    
+        new_alpha = sum_sig/(sum_reg/self.alpha)
+        print "Suggested alpha =", 0.1*new_alpha
+        
+        try:
+            err = math.fabs(chi2/(npts-nfunc))* inv_cov
+        except:
+            print "Error estimating uncertainties"
+            
+        
+        return c, err
+        
+    def svd(self, nfunc=5):
+        import math, time
+        # Ac - b = 0
+        
+        A = numpy.zeros([nfunc, nfunc])
+        y = numpy.zeros(nfunc)
+        
+        t_0 = time.time()
+        for i in range(nfunc):
+            # A
+            for j in range(nfunc):
+                A[i][j] = 0.0
+                for k in range(len(self.x)):
+                    err = self.err[k]
+                    A[i][j] += 1.0/err/err*self.basefunc_ft(self.d_max, j+1, self.x[k]) \
+                            *self.basefunc_ft(self.d_max, i+1, self.x[k]);
+                #print A[i][j]
+                #A[i][j] -= self.alpha*(math.cos(math.pi*(i+j)) - math.cos(math.pi*(i-j)));
+                if i==j:
+                    A[i][j] += -1.0*self.alpha
+                elif i-j==1 or i-j==-1:
+                    A[i][j] += 1.0*self.alpha
+                #print "   ",A[i][j]
+            # y
+            y[i] = 0.0
+            for k in range(len(self.x)):
+                y[i] = self.y[k]/self.err[k]/self.err[k]*self.basefunc_ft(self.d_max, i+1, self.x[k])
+            
+        print time.time()-t_0, 'secs'
+        
+        # use numpy.pinv(A)
+        #inv_A = numpy.linalg.inv(A)
+        #c = y*inv_A
+        print y
+        c = numpy.linalg.solve(A, y)
+        
+        
+        print c
+        
+        err = numpy.zeros(len(c))
+        return c, err
+        
+        
+        
         
     
Index: pr_inversion/test/utest_invertor.py
===================================================================
--- pr_inversion/test/utest_invertor.py	(revision eca05c8ee390007a16bb46b22163b8750ee5ab98)
+++ pr_inversion/test/utest_invertor.py	(revision 2d06beb109031c888260d5a73af0ae5e57628b01)
@@ -174,4 +174,66 @@
             print "chi2 =", chi2/51.0
             raise
+        
+    def test_lstsq(self):
+        """
+            Test an inversion for which we know the answer
+        """
+        x, y, err = load("sphere_80.txt")
+
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = .0007
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        #out, cov = self.invertor.invert(10)
+        
+        out, cov = self.invertor.lstsq(10)
+         
+        
+        # This is a very specific case
+        # We should make sure it always passes
+        try:
+            self.assertTrue(self.invertor.chi2/len(x)<200.00)
+        except:
+            print "Chi2(I(q)) =", self.invertor.chi2/len(x)
+            raise
+        
+        # Check the computed P(r) with the theory
+        # for shpere of radius 80
+        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0)
+        y = numpy.zeros(len(x))
+        dy = numpy.zeros(len(x))
+        y_true = numpy.zeros(len(x))
+
+        sum = 0.0
+        sum_true = 0.0
+        for i in range(len(x)):
+            #y[i] = self.invertor.pr(out, x[i])
+            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i])
+            sum += y[i]
+            if x[i]<80.0:
+                y_true[i] = pr_theory(x[i], 80.0)
+            else:
+                y_true[i] = 0
+            sum_true += y_true[i]
+            
+        y = y/sum*self.invertor.d_max/len(x)
+        dy = dy/sum*self.invertor.d_max/len(x)
+        y_true = y_true/sum_true*self.invertor.d_max/len(x)
+        
+        chi2 = 0.0
+        for i in range(len(x)):
+            res = (y[i]-y_true[i])/dy[i]
+            chi2 += res*res
+            
+        try:
+            self.assertTrue(chi2/51.0<50.0)
+        except:
+            print "chi2(P(r)) =", chi2/51.0
+            raise
             
     def test_q_zero(self):
@@ -239,7 +301,37 @@
             print "Chi2 =", self.invertor.chi2
             raise
+        
+    def no_test_time(self):
+        x, y, err = load("sphere_80.txt")
+
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = 1e-7
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+    
+        # time scales like nfunc**2
+        # on a Lenovo Intel Core 2 CPU T7400 @ 2.16GHz, 
+        # I get time/(nfunc)**2 = 0.022 sec
                             
-        
-
+        out, cov = self.invertor.invert(15)
+        t16 = self.invertor.elapsed
+        
+        out, cov = self.invertor.invert(30)
+        t30 = self.invertor.elapsed
+        
+        t30s = t30/30.0**2
+        self.assertTrue( (t30s-t16/16.0**2)/t30s <1.2 )
+        
+    def test_clone(self):
+        self.invertor.x = self.x_in
+        clone = self.invertor.clone()
+        
+        for i in range(len(self.x_in)):
+            self.assertEqual(self.x_in[i], clone.x[i])
+        
 def pr_theory(r, R):
     """
