Changeset 9a23253e in sasview

Timestamp:

Jun 30, 2008 2:18:00 PM (16 years ago)

Author:

Mathieu Doucet <doucetm@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

a4bd2ac

Parents:

7bb61da

Message:

Started to add slit smearing. Added Rg and I(q=0) as outputs.

Location:

pr_inversion

Files:

• c_extensions/Cinvertor.c (modified) (43 diffs)
• c_extensions/invertor.c (modified) (15 diffs)
• c_extensions/invertor.h (modified) (3 diffs)
• invertor.py (modified) (13 diffs)
• test/test_output.txt (modified) (1 diff)
• test/utest_invertor.py (modified) (2 diffs)

pr_inversion/c_extensions/Cinvertor.c

@@ -3 +3 @@
  * Cinvertor is the base class for the Invertor class
  * and provides the underlying computations.
- * 
+ *
  */
 #include <Python.h>
@@ -38 +38 @@
  * Cinvertor is the base class for the Invertor class
  * and provides the underlying computations.
- * 
+ *
  */
 typedef struct {
-    PyObject_HEAD    
+    PyObject_HEAD
     /// Internal data structure
-    Invertor_params params; 
+    Invertor_params params;
 } Cinvertor;
 
@@ -51 +51 @@
 {
     invertor_dealloc(&(self->params));
-     
+
     self->ob_type->tp_free((PyObject*)self);
 
@@ -60 +60 @@
 {
     Cinvertor *self;
-    
+
     self = (Cinvertor *)type->tp_alloc(type, 0);
-   
+
     return (PyObject *)self;
 }
@@ -69 +69 @@
 Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds)
 {
-    if (self != NULL) {         
+    if (self != NULL) {
         // Create parameters
         invertor_init(&(self->params));
@@ -82 +82 @@
 };
 
-const char set_x_doc[] = 
+const char set_x_doc[] =
         "Function to set the x data\n"
         "Takes an array of doubles as input.\n"
@@ -97 +97 @@
         double *data;
         int i;
-  
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,data,ndata);
-        
+
         free(self->params.x);
         self->params.x = (double*) malloc(ndata*sizeof(double));
-        
+
         if(self->params.x==NULL) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.set_x: problem allocating memory.");
-                return NULL;            
-        }
-        
+                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);        
-}
-
-const char get_x_doc[] = 
+        return Py_BuildValue("i", self->params.npoints);
+}
+
+const char get_x_doc[] =
         "Function to get the x data\n"
         "Takes an array of doubles as input.\n"
@@ -129 +129 @@
         double *data;
     int i;
-    
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj, data, ndata);
-        
+
         // Check that the input array is large enough
         if (ndata < self->params.npoints) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.get_x: input array too short for data.");
-                return NULL;            
-        }
-        
+                return NULL;
+        }
+
         for(i=0; i<self->params.npoints; i++){
                 data[i] = self->params.x[i];
         }
-        
-        return Py_BuildValue("i", self->params.npoints);        
-}
-
-const char set_y_doc[] = 
+
+        return Py_BuildValue("i", self->params.npoints);
+}
+
+const char set_y_doc[] =
         "Function to set the y data\n"
         "Takes an array of doubles as input.\n"
@@ -162 +162 @@
         double *data;
         int i;
-  
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,data,ndata);
-        
+
         free(self->params.y);
         self->params.y = (double*) malloc(ndata*sizeof(double));
-        
+
         if(self->params.y==NULL) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.set_y: problem allocating memory.");
-                return NULL;            
-        }
-        
+                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);     
-}
-
-const char get_y_doc[] = 
+        return Py_BuildValue("i", self->params.ny);
+}
+
+const char get_y_doc[] =
         "Function to get the y data\n"
         "Takes an array of doubles as input.\n"
@@ -194 +194 @@
         double *data;
     int i;
-    
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj, data, ndata);
-        
+
         // Check that the input array is large enough
         if (ndata < self->params.ny) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.get_y: input array too short for data.");
-                return NULL;            
-        }
-        
+                return NULL;
+        }
+
         for(i=0; i<self->params.ny; i++){
                 data[i] = self->params.y[i];
         }
-        
-        return Py_BuildValue("i", self->params.npoints);        
-}
-
-const char set_err_doc[] = 
+
+        return Py_BuildValue("i", self->params.npoints);
+}
+
+const char set_err_doc[] =
         "Function to set the err data\n"
         "Takes an array of doubles as input.\n"
@@ -227 +227 @@
         double *data;
         int i;
-  
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,data,ndata);
-        
+
         free(self->params.err);
         self->params.err = (double*) malloc(ndata*sizeof(double));
-        
+
         if(self->params.err==NULL) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.set_err: problem allocating memory.");
-                return NULL;            
-        }
-        
+                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);   
-}
-
-const char get_err_doc[] = 
+        return Py_BuildValue("i", self->params.nerr);
+}
+
+const char get_err_doc[] =
         "Function to get the err data\n"
         "Takes an array of doubles as input.\n"
@@ -259 +259 @@
         double *data;
     int i;
-    
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj, data, ndata);
-        
+
         // Check that the input array is large enough
         if (ndata < self->params.nerr) {
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.get_err: input array too short for data.");
-                return NULL;            
-        }
-        
+                return NULL;
+        }
+
         for(i=0; i<self->params.nerr; i++){
                 data[i] = self->params.err[i];
         }
-        
-        return Py_BuildValue("i", self->params.npoints);        
-}
-
-const char is_valid_doc[] = 
+
+        return Py_BuildValue("i", self->params.npoints);
+}
+
+const char is_valid_doc[] =
         "Check the validity of the stored data\n"
         " @return: Returns the number of points if it's all good, -1 otherwise";
@@ -291 +291 @@
         } else {
                 return Py_BuildValue("i", -1);
-        }       
-}
-
-const char set_dmax_doc[] = 
+        }
+}
+
+const char set_has_bck_doc[] =
+        "Sets background flag\n";
+
+/**
+ * Sets the maximum distance
+ */
+static PyObject * set_has_bck(Cinvertor *self, PyObject *args) {
+        int has_bck;
+
+        if (!PyArg_ParseTuple(args, "i", &has_bck)) return NULL;
+        self->params.has_bck = has_bck;
+        return Py_BuildValue("i", self->params.has_bck);
+}
+
+const char get_has_bck_doc[] =
+        "Gets background flag\n";
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_has_bck(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("i", self->params.has_bck);
+}
+
+const char set_dmax_doc[] =
         "Sets the maximum distance\n";
 
@@ -302 +326 @@
 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);  
-}
-
-const char get_dmax_doc[] = 
+        return Py_BuildValue("d", self->params.d_max);
+}
+
+const char get_dmax_doc[] =
         "Gets the maximum distance\n";
 
@@ -315 +339 @@
  */
 static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("d", self->params.d_max);  
-}
-
-const char set_qmin_doc[] = 
+        return Py_BuildValue("d", self->params.d_max);
+}
+
+const char set_slit_height_doc[] =
+        "Sets the slit height in units of q [A-1]\n";
+
+/**
+ * Sets the slit height
+ */
+static PyObject * set_slit_height(Cinvertor *self, PyObject *args) {
+        double slit_height;
+
+        if (!PyArg_ParseTuple(args, "d", &slit_height)) return NULL;
+        self->params.slit_height = slit_height;
+        return Py_BuildValue("d", self->params.slit_height);
+}
+
+const char get_slit_height_doc[] =
+        "Gets the slit height\n";
+
+/**
+ * Gets the slit height
+ */
+static PyObject * get_slit_height(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("d", self->params.slit_height);
+}
+
+const char set_slit_width_doc[] =
+        "Sets the slit width in units of q [A-1]\n";
+
+/**
+ * Sets the slit width
+ */
+static PyObject * set_slit_width(Cinvertor *self, PyObject *args) {
+        double slit_width;
+
+        if (!PyArg_ParseTuple(args, "d", &slit_width)) return NULL;
+        self->params.slit_width = slit_width;
+        return Py_BuildValue("d", self->params.slit_width);
+}
+
+const char get_slit_width_doc[] =
+        "Gets the slit width\n";
+
+/**
+ * Gets the slit width
+ */
+static PyObject * get_slit_width(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("d", self->params.slit_width);
+}
+
+
+const char set_qmin_doc[] =
         "Sets the minimum q\n";
 
@@ -326 +399 @@
 static PyObject * set_qmin(Cinvertor *self, PyObject *args) {
         double q_min;
-  
+
         if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL;
         self->params.q_min = q_min;
-        return Py_BuildValue("d", self->params.q_min);  
-}
-
-const char get_qmin_doc[] = 
+        return Py_BuildValue("d", self->params.q_min);
+}
+
+const char get_qmin_doc[] =
         "Gets the minimum q\n";
 
@@ -339 +412 @@
  */
 static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("d", self->params.q_min);  
-}
-
-const char set_qmax_doc[] = 
+        return Py_BuildValue("d", self->params.q_min);
+}
+
+const char set_qmax_doc[] =
         "Sets the maximum q\n";
 
@@ -350 +423 @@
 static PyObject * set_qmax(Cinvertor *self, PyObject *args) {
         double q_max;
-  
+
         if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL;
         self->params.q_max = q_max;
-        return Py_BuildValue("d", self->params.q_max);  
-}
-
-const char get_qmax_doc[] = 
+        return Py_BuildValue("d", self->params.q_max);
+}
+
+const char get_qmax_doc[] =
         "Gets the maximum q\n";
 
@@ -363 +436 @@
  */
 static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("d", self->params.q_max);  
-}
-
-const char set_alpha_doc[] = 
+        return Py_BuildValue("d", self->params.q_max);
+}
+
+const char set_alpha_doc[] =
         "Sets the alpha parameter\n";
 
 static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
         double alpha;
-  
+
         if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
         self->params.alpha = alpha;
-        return Py_BuildValue("d", self->params.alpha);  
-}
-
-const char get_alpha_doc[] = 
+        return Py_BuildValue("d", self->params.alpha);
+}
+
+const char get_alpha_doc[] =
         "Gets the alpha parameter\n";
 
@@ -384 +457 @@
  */
 static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("d", self->params.alpha);  
-}
-
-const char get_nx_doc[] = 
+        return Py_BuildValue("d", self->params.alpha);
+}
+
+const char get_nx_doc[] =
         "Gets the number of x points\n";
 
@@ -394 +467 @@
  */
 static PyObject * get_nx(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("i", self->params.npoints);        
-}
-
-const char get_ny_doc[] = 
+        return Py_BuildValue("i", self->params.npoints);
+}
+
+const char get_ny_doc[] =
         "Gets the number of y points\n";
 
@@ -404 +477 @@
  */
 static PyObject * get_ny(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("i", self->params.ny);     
-}
-
-const char get_nerr_doc[] = 
+        return Py_BuildValue("i", self->params.ny);
+}
+
+const char get_nerr_doc[] =
         "Gets the number of err points\n";
 
@@ -414 +487 @@
  */
 static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
-        return Py_BuildValue("i", self->params.nerr);   
-}
-
-
-const char residuals_doc[] = 
+        return Py_BuildValue("i", self->params.nerr);
+}
+
+
+const char residuals_doc[] =
         "Function to call to evaluate the residuals\n"
         "for P(r) inversion\n"
@@ -441 +514 @@
         // Number of slices in regularization term estimate
         int nslice = 25;
-        
+
         PyObject *data_obj;
         Py_ssize_t npars;
-          
-        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
-        
+
+        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+
         OUTVECTOR(data_obj,pars,npars);
-                
+
     // PyList of residuals
         // Should create this list only once and refill it
@@ -454 +527 @@
 
     regterm = reg_term(pars, self->params.d_max, npars, nslice);
-    
+
     for(i=0; i<self->params.npoints; 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]);
         tmp = residual;
-        
+
         // regularization term
         residual += self->params.alpha * regterm;
-        
+
         if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.residuals: error setting residual.");
                 return NULL;
         };
-                
+
     }
-    
+
         return residuals;
 }
 
-const char pr_residuals_doc[] = 
+const char pr_residuals_doc[] =
         "Function to call to evaluate the residuals\n"
         "for P(r) minimization (for testing purposes)\n"
@@ -498 +571 @@
         // Number of slices in regularization term estimate
         int nslice = 25;
-        
+
         PyObject *data_obj;
         Py_ssize_t npars;
-          
-        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
-        
+
+        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+
         OUTVECTOR(data_obj,pars,npars);
-                
+
         // Should create this list only once and refill it
     residuals = PyList_New(self->params.npoints);
@@ -511 +584 @@
     regterm = reg_term(pars, self->params.d_max, npars, nslice);
 
-    
+
     for(i=0; i<self->params.npoints; i++) {
         diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]);
         residual = diff*diff / (self->params.err[i]*self->params.err[i]);
         tmp = residual;
-        
+
         // regularization term
         residual += self->params.alpha * regterm;
-        
+
         if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
-            PyErr_SetString(CinvertorError, 
+            PyErr_SetString(CinvertorError,
                 "Cinvertor.residuals: error setting residual.");
                 return NULL;
         };
-                
+
     }
-    
+
         return residuals;
 }
 
-const char get_iq_doc[] = 
+const char get_iq_doc[] =
         "Function to call to evaluate the scattering intensity\n"
         " @param args: c-parameters, and q\n"
@@ -546 +619 @@
         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);    
-}
-
-const char get_pr_doc[] = 
+        return Py_BuildValue("f", iq_value);
+}
+
+const char get_pr_doc[] =
         "Function to call to evaluate P(r)\n"
         " @param args: c-parameters and r\n"
@@ -569 +642 @@
         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);    
-}
-
-const char get_pr_err_doc[] = 
+        return Py_BuildValue("f", pr_value);
+}
+
+const char get_pr_err_doc[] =
         "Function to call to evaluate P(r) with errors\n"
         " @param args: c-parameters and r\n"
@@ -596 +669 @@
         PyObject *err_obj;
         Py_ssize_t npars2;
-        int i; 
-        
+        int i;
+
         if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
-        OUTVECTOR(data_obj,pars,npars); 
+        OUTVECTOR(data_obj,pars,npars);
         OUTVECTOR(err_obj,pars_err,npars2);
 
         pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value);
-        return Py_BuildValue("ff", pr_value, pr_err_value);     
-}
-
-const char basefunc_ft_doc[] = 
+        return Py_BuildValue("ff", pr_value, pr_err_value);
+}
+
+const char basefunc_ft_doc[] =
         "Returns the value of the nth Fourier transofrmed base function\n"
         " @param args: c-parameters, n and q\n"
@@ -614 +687 @@
         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));      
-        
-}
-
-const char oscillations_doc[] = 
+        return Py_BuildValue("f", ortho_transformed(d_max, n, q));
+
+}
+
+const char oscillations_doc[] =
         "Returns the value of the oscillation figure of merit for\n"
         "the given set of coefficients. For a sphere, the oscillation\n"
@@ -632 +705 @@
         Py_ssize_t npars;
         double oscill, norm;
-        
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,pars,npars);
-        
+
         oscill = reg_term(pars, self->params.d_max, npars, 100);
         norm   = int_p2(pars, self->params.d_max, npars, 100);
-        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );    
-        
-}
-
-const char get_peaks_doc[] = 
+        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );
+
+}
+
+const char get_peaks_doc[] =
         "Returns the number of peaks in the output P(r) distrubution\n"
         "for the given set of coefficients.\n"
@@ -653 +726 @@
         Py_ssize_t npars;
         int count;
-        
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,pars,npars);
-        
+
         count = npeaks(pars, self->params.d_max, npars, 100);
 
-        return Py_BuildValue("i", count );      
-        
-}
-
-const char get_positive_doc[] = 
+        return Py_BuildValue("i", count );
+
+}
+
+const char get_positive_doc[] =
         "Returns the fraction of P(r) that is positive over\n"
         "the full range of r for the given set of coefficients.\n"
@@ -674 +747 @@
         Py_ssize_t npars;
         double fraction;
-         
+
         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
         OUTVECTOR(data_obj,pars,npars);
-        
+
         fraction = positive_integral(pars, self->params.d_max, npars, 100);
 
-        return Py_BuildValue("f", fraction );   
-        
-}
-
-const char get_pos_err_doc[] = 
+        return Py_BuildValue("f", fraction );
+
+}
+
+const char get_pos_err_doc[] =
         "Returns the fraction of P(r) that is 1 standard deviation\n"
         "above zero over the full range of r for the given set of coefficients.\n"
@@ -698 +771 @@
         Py_ssize_t npars2;
         double fraction;
-        
+
         if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL;
-        OUTVECTOR(data_obj,pars,npars); 
+        OUTVECTOR(data_obj,pars,npars);
         OUTVECTOR(err_obj,pars_err,npars2);
-        
+
         fraction = positive_errors(pars, pars_err, self->params.d_max, npars, 51);
 
-        return Py_BuildValue("f", fraction );   
-        
-}
-
+        return Py_BuildValue("f", fraction );
+
+}
+
+const char get_rg_doc[] =
+        "Returns the value of the radius of gyration Rg.\n"
+        " @param args: c-parameters\n"
+        " @return: Rg";
+
+static PyObject * get_rg(Cinvertor *self, PyObject *args) {
+        double *pars;
+        double *pars_err;
+        PyObject *data_obj;
+        Py_ssize_t npars;
+        Py_ssize_t npars2;
+        double value;
+
+        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+        OUTVECTOR(data_obj,pars,npars);
+
+        value = rg(pars, self->params.d_max, npars, 101);
+
+        return Py_BuildValue("f", value );
+
+}
+
+const char get_iq0_doc[] =
+        "Returns the value of I(q=0).\n"
+        " @param args: c-parameters\n"
+        " @return: I(q=0)";
+
+static PyObject * get_iq0(Cinvertor *self, PyObject *args) {
+        double *pars;
+        double *pars_err;
+        PyObject *data_obj;
+        Py_ssize_t npars;
+        Py_ssize_t npars2;
+        double value;
+
+        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+        OUTVECTOR(data_obj,pars,npars);
+
+        value = 4.0*acos(-1.0)*int_pr(pars, self->params.d_max, npars, 101);
+
+        return Py_BuildValue("f", value );
+
+}
+
+/**
+ * Check whether a q-value is within acceptabel limits
+ * Return 1 if accepted, 0 if rejected.
+ */
+int accept_q(Cinvertor *self, double q) {
+    if (self->params.q_min>0 && q<self->params.q_min) return 0;
+    if (self->params.q_max>0 && q>self->params.q_max) return 0;
+    return 1;
+}
+
+const char get_matrix_doc[] =
+        "Returns A matrix and b vector for least square problem.\n"
+        " @param nfunc: number of base functions\n"
+        " @param nr: number of r-points used when evaluating reg term.\n"
+        " @param a: A array to fill\n"
+        " @param b: b vector to fill\n"
+        " @return: 0";
+
+static PyObject * get_matrix(Cinvertor *self, PyObject *args) {
+        double *a;
+        double *b;
+        PyObject *a_obj;
+        PyObject *b_obj;
+        Py_ssize_t n_a;
+        Py_ssize_t n_b;
+        // Number of bins for regularization term evaluation
+        int nr, nfunc;
+        int i, j, i_r;
+        double r, sqrt_alpha, pi;
+        double tmp;
+        int offset;
+
+        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &b_obj)) return NULL;
+        OUTVECTOR(a_obj,a,n_a);
+        OUTVECTOR(b_obj,b,n_b);
+
+        assert(n_b>=nfunc);
+        assert(n_a>=nfunc*(nr*self->params.npoints));
+
+        sqrt_alpha = sqrt(self->params.alpha);
+        pi = acos(-1.0);
+        offset = (self->params.has_bck==1) ? 0 : 1;
+
+    for (j=0; j<nfunc; j++) {
+        for (i=0; i<self->params.npoints; i++) {
+            if (accept_q(self, self->params.x[i])){
+                if (self->params.has_bck==1 && j==0) {
+                    a[i*nfunc+j] = 1.0/self->params.err[i];
+                } else {
+                        a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i];
+                }
+            }
+        }
+        for (i_r=0; i_r<nr; i_r++){
+            if (self->params.has_bck==1 && j==0) {
+                a[(i_r+self->params.npoints)*nfunc+j] = 0.0;
+            } else {
+                    r = self->params.d_max/nr*i_r;
+                    tmp = pi*(j+offset)/self->params.d_max;
+                    a[(i_r+self->params.npoints)*nfunc+j] = sqrt_alpha * 1.0/nr*self->params.d_max*2.0*
+                    (2.0*pi*(j+offset)/self->params.d_max*cos(pi*(j+offset)*r/self->params.d_max) +
+                    tmp*tmp*r * sin(pi*(j+offset)*r/self->params.d_max));
+                }
+        }
+    }
+
+    for (i=0; i<self->params.npoints; i++) {
+        if (accept_q(self, self->params.x[i])){
+            b[i] = self->params.y[i]/self->params.err[i];
+         }
+        }
+
+        return Py_BuildValue("i", 0);
+
+}
+
+const char get_invcov_matrix_doc[] =
+        " Compute the inverse covariance matrix, defined as inv_cov = a_transposed x a.\n"
+        " @param nfunc: number of base functions\n"
+        " @param nr: number of r-points used when evaluating reg term.\n"
+        " @param a: A array to fill\n"
+        " @param inv_cov: inverse covariance array to be filled\n"
+        " @return: 0";
+
+static PyObject * get_invcov_matrix(Cinvertor *self, PyObject *args) {
+        double *a;
+        PyObject *a_obj;
+        Py_ssize_t n_a;
+        double *inv_cov;
+        PyObject *cov_obj;
+        Py_ssize_t n_cov;
+        int nr, nfunc;
+        int i, j, k;
+
+        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &cov_obj)) return NULL;
+        OUTVECTOR(a_obj,a,n_a);
+        OUTVECTOR(cov_obj,inv_cov,n_cov);
+
+        assert(n_cov>=nfunc*nfunc);
+        assert(n_a>=nfunc*(nr+self->params.npoints));
+
+        for (i=0; i<nfunc; i++) {
+                for (j=0; j<nfunc; j++) {
+                        inv_cov[i*nfunc+j] = 0.0;
+                        for (k=0; k<nr+self->params.npoints; k++) {
+                                inv_cov[i*nfunc+j] += a[k*nfunc+i]*a[k*nfunc+j];
+                        }
+                }
+        }
+        return Py_BuildValue("i", 0);
+}
+
+const char get_reg_size_doc[] =
+        " Compute the covariance matrix, defined as inv_cov = a_transposed x a.\n"
+        " @param nfunc: number of base functions\n"
+        " @param nr: number of r-points used when evaluating reg term.\n"
+        " @param a: A array to fill\n"
+        " @param inv_cov: inverse covariance array to be filled\n"
+        " @return: 0";
+
+static PyObject * get_reg_size(Cinvertor *self, PyObject *args) {
+        double *a;
+        PyObject *a_obj;
+        Py_ssize_t n_a;
+        int nr, nfunc;
+        int i, j;
+        double sum_sig, sum_reg;
+
+        if (!PyArg_ParseTuple(args, "iiO", &nfunc, &nr, &a_obj)) return NULL;
+        OUTVECTOR(a_obj,a,n_a);
+
+        assert(n_a>=nfunc*(nr+self->params.npoints));
+
+    sum_sig = 0.0;
+    sum_reg = 0.0;
+    for (j=0; j<nfunc; j++){
+        for (i=0; i<self->params.npoints; i++){
+            if (accept_q(self, self->params.x[i])==1)
+                sum_sig += (a[i*nfunc+j])*(a[i*nfunc+j]);
+        }
+        for (i=0; i<nr; i++){
+            sum_reg += (a[(i+self->params.npoints)*nfunc+j])*(a[(i+self->params.npoints)*nfunc+j]);
+        }
+    }
+    return Py_BuildValue("ff", sum_sig, sum_reg);
+}
 
 const char eeeget_qmin_doc[] = "\
@@ -714 +977 @@
 \n\
 This is the second line.";
-const char eeeset_qmin_doc[] = 
+const char eeeset_qmin_doc[] =
         "This is a multiline doc string.\n"
         "\n"
@@ -736 +999 @@
                    {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, set_alpha_doc},
                    {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, get_alpha_doc},
+                   {"set_slit_width", (PyCFunction)set_slit_width, METH_VARARGS, set_slit_width_doc},
+                   {"get_slit_width", (PyCFunction)get_slit_width, METH_VARARGS, get_slit_width_doc},
+                   {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc},
+                   {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc},
+                   {"set_has_bck", (PyCFunction)set_has_bck, METH_VARARGS, set_has_bck_doc},
+                   {"get_has_bck", (PyCFunction)get_has_bck, METH_VARARGS, get_has_bck_doc},
                    {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc},
                    {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc},
@@ -748 +1017 @@
                    {"get_positive", (PyCFunction)get_positive, METH_VARARGS, get_positive_doc},
                    {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, get_pos_err_doc},
-   
+                   {"rg", (PyCFunction)get_rg, METH_VARARGS, get_rg_doc},
+                   {"iq0", (PyCFunction)get_iq0, METH_VARARGS, get_iq0_doc},
+                   {"_get_matrix", (PyCFunction)get_matrix, METH_VARARGS, get_matrix_doc},
+                   {"_get_invcov_matrix", (PyCFunction)get_invcov_matrix, METH_VARARGS, get_invcov_matrix_doc},
+                   {"_get_reg_size", (PyCFunction)get_reg_size, METH_VARARGS, get_reg_size_doc},
+
    {NULL}
 };
@@ -796 +1070 @@
 
 static PyMethodDef module_methods[] = {
-    {NULL} 
+    {NULL}
 };
 
@@ -802 +1076 @@
  * 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;
@@ -811 +1085 @@
     Py_INCREF(&CinvertorType);
     PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
-    
+
     d = PyModule_GetDict(module);
     CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL);
@@ -822 +1096 @@
 #endif
 PyMODINIT_FUNC
-initpr_inversion(void) 
+initpr_inversion(void)
 {
     PyObject* m;
@@ -828 +1102 @@
     m = Py_InitModule3("pr_inversion", module_methods,
                        "C extension module for inversion to P(r).");
-                       
+
     addCinvertor(m);
 }

pr_inversion/c_extensions/invertor.c

@@ -20 +20 @@
         pars->q_min = -1.0;
         pars->q_max = -1.0;
+        pars->has_bck = 0;
 }
 
@@ -25 +26 @@
 /**
  * P(r) of a sphere, for test purposes
- * 
+ *
  * @param R: radius of the sphere
  * @param r: distance, in the same units as the radius
@@ -41 +42 @@
  * Orthogonal functions:
  * B(r) = 2r sin(pi*nr/d)
- * 
+ *
  */
 double ortho(double d_max, int n, double r) {
@@ -49 +50 @@
 /**
  * Fourier transform of the nth orthogonal function
- * 
+ *
  */
 double ortho_transformed(double d_max, int n, double q) {
@@ -57 +58 @@
 
 /**
+ * Slit-smeared Fourier transform of the nth orthogonal function.
+ * Smearing follows Lake, Acta Cryst. (1967) 23, 191.
+ */
+double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) {
+        double sum, value, y, z;
+        int i, j;
+        double fnpts;
+        sum = 0.0;
+        fnpts = (float)npts-1.0;
+
+        for(i=0; i<npts; i++) {
+                y = -width/2.0+width/fnpts*(float)i;
+                for(j=0; j<npts; j++) {
+                        z = height/fnpts*(float)j;
+                        sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z));
+                }
+        }
+
+        return sum/npts/npts/height/width;
+}
+
+/**
  * First derivative in of the orthogonal function dB(r)/dr
- * 
+ *
  */
 double ortho_derived(double d_max, int n, double r) {
@@ -80 +103 @@
  */
 double pr(double *pars, double d_max, int n_c, double r) {
-    double sum = 0.0;  
+    double sum = 0.0;
         int i;
     for (i=0; i<n_c; i++) {
@@ -91 +114 @@
  * P(r) calculated from the expansion, with errors
  */
-void pr_err(double *pars, double *err, double d_max, int n_c, 
+void pr_err(double *pars, double *err, double d_max, int n_c,
                 double r, double *pr_value, double *pr_value_err) {
-    double sum = 0.0; 
+    double sum = 0.0;
     double sum_err = 0.0;
     double func_value;
@@ -109 +132 @@
         *pr_value_err = sum;
     }
-} 
+}
 
 /**
@@ -115 +138 @@
  */
 double dprdr(double *pars, double d_max, int n_c, double r) {
-    double sum = 0.0; 
-        int i; 
+    double sum = 0.0;
+        int i;
     for (i=0; i<n_c; i++) {
         sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max));
@@ -127 +150 @@
  */
 double reg_term(double *pars, double d_max, int n_c, int nslice) {
-    double sum = 0.0; 
+    double sum = 0.0;
     double r;
     double deriv;
@@ -143 +166 @@
  */
 double int_p2(double *pars, double d_max, int n_c, int nslice) {
-    double sum = 0.0; 
-    double r; 
+    double sum = 0.0;
+    double r;
     double value;
         int i;
@@ -156 +179 @@
 
 /**
+ * Integral of P(r)
+ */
+double int_pr(double *pars, double d_max, int n_c, int nslice) {
+    double sum = 0.0;
+    double r;
+    double value;
+        int i;
+    for (i=0; i<nslice; i++) {
+        r = d_max/(1.0*nslice)*i;
+        value = pr(pars, d_max, n_c, r);
+        sum += value;
+    }
+    return sum/(1.0*nslice)*d_max;
+}
+
+/**
  * Get the number of P(r) peaks.
  */
 int npeaks(double *pars, double d_max, int n_c, int nslice) {
-    double r; 
+    double r;
     double value;
         int i;
@@ -187 +226 @@
  */
 double positive_integral(double *pars, double d_max, int n_c, int nslice) {
-    double r; 
+    double r;
     double value;
         int i;
         double sum_pos = 0.0;
         double sum = 0.0;
-        
+
     for (i=0; i<nslice; i++) {
         r = d_max/(1.0*nslice)*i;
@@ -207 +246 @@
  */
 double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) {
-    double r; 
-    double value;
-        int i; 
+    double r;
+    double value;
+        int i;
         double sum_pos = 0.0;
         double sum = 0.0;
         double pr_val;
         double pr_val_err;
-        
+
     for (i=0; i<nslice; i++) {
         r = d_max/(1.0*nslice)*i;
@@ -220 +259 @@
         if (pr_val>pr_val_err) sum_pos += pr_val;
         sum += fabs(pr_val);
-        
+
 
     }
     return sum_pos/sum;
 }
+
+/**
+ * R_g radius of gyration calculation
+ *
+ * R_g**2 = integral[r**2 * p(r) dr] /  (2.0 * integral[p(r) dr])
+ */
+double rg(double *pars, double d_max, int n_c, int nslice) {
+    double sum_r2 = 0.0;
+    double sum    = 0.0;
+    double r;
+    double value;
+        int i;
+    for (i=0; i<nslice; i++) {
+        r = d_max/(1.0*nslice)*i;
+        value = pr(pars, d_max, n_c, r);
+        sum += value;
+        sum_r2 += r*r*value;
+    }
+    return sqrt(sum_r2/(2.0*sum));
+}
+

pr_inversion/c_extensions/invertor.h

@@ -15 +15 @@
     double *err;
     /// Number of q points
-    int npoints;    
+    int npoints;
     /// Number of I(q) points
-    int ny;    
+    int ny;
     /// Number of dI(q) points
-    int nerr;  
+    int nerr;
     /// Alpha value
     double alpha;
@@ -26 +26 @@
     /// Maximum q to include in inversion
     double q_max;
-} Invertor_params; 
+    /// Flag for whether or not to evalute a constant background while inverting
+    int has_bck;
+    /// Slit height in units of q [A-1]
+    double slit_height;
+    /// Slit width in units of q [A-1]
+    double slit_width;
+} Invertor_params;
 
 void invertor_dealloc(Invertor_params *pars);
 
 void invertor_init(Invertor_params *pars);
-
 
 double pr_sphere(double R, double r);
@@ -42 +47 @@
 double reg_term(double *pars, double d_max, int n_c, int nslice);
 double int_p2(double *pars, double d_max, int n_c, int nslice);
-void pr_err(double *pars, double *err, double d_max, int n_c, 
+void pr_err(double *pars, double *err, double d_max, int n_c,
                 double r, double *pr_value, double *pr_value_err);
 int npeaks(double *pars, double d_max, int n_c, int nslice);
 double positive_integral(double *pars, double d_max, int n_c, int nslice);
 double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice);
+double rg(double *pars, double d_max, int n_c, int nslice);
+double int_pr(double *pars, double d_max, int n_c, int nslice);
+double ortho_transformed_smeared(double d_max, int n, double heigth, double width, double q, int npts);
 
 #endif

pr_inversion/invertor.py

@@ -63 +63 @@
         - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero
     """
+    #TODO: Allow for slit smearing. Smear each base function once before filling
+    # the A matrix.
+    
     ## Chisqr of the last computation
     chi2  = 0
@@ -75 +78 @@
     ## Last errors on output values
     cov = None
+    ## Flag to allow I(q) data with constant background
+    #has_bck = False
+    ## Background value
+    background = 0
+    
     
     def __init__(self):
@@ -106 +114 @@
         elif name=='alpha':
             return self.set_alpha(value)
+        elif name=='slit_height':
+            return self.set_slit_height(value)
+        elif name=='slit_width':
+            return self.set_slit_width(value)
+        elif name=='has_bck':
+            if value==True:
+                return self.set_has_bck(1)
+            elif value==False:
+                return self.set_has_bck(0)
+            else:
+                raise ValueError, "Invertor: has_bck can only be True or False"
             
         return Cinvertor.__setattr__(self, name, value)
@@ -142 +161 @@
         elif name=='alpha':
             return self.get_alpha()
+        elif name=='slit_height':
+            return self.get_slit_height()
+        elif name=='slit_width':
+            return self.get_slit_width()
+        elif name=='has_bck':
+            value = self.get_has_bck()
+            if value==1:
+                return True
+            else:
+                return False
         elif name in self.__dict__:
             return self.__dict__[name]
@@ -161 +190 @@
         invertor.y = self.y
         invertor.err = self.err
+        invertor.has_bck = self.has_bck
         
         return invertor
@@ -194 +224 @@
             @return: c_out, c_cov - the coefficients with covariance matrix 
         """
-        #TODO: call the pyhton implementation for now. In the future, translate this to C.
+        # Reset the background value before proceeding
+        self.background = 0.0
         return self.lstsq(nfunc, nr=nr)
+    
+    def iq(self, out, q):
+        """
+            Function to call to evaluate the scattering intensity
+            @param args: c-parameters, and q
+            @return: I(q)
+        """
+        return Cinvertor.iq(self, out, q)+self.background
     
     def invert_optimize(self, nfunc=10, nr=20):
@@ -325 +364 @@
 
         """
-        import math
+        #TODO: Allow for background by starting at n=0 (since the base function
+        # is zero for n=0).
+        import math, time
         from scipy.linalg.basic import lstsq
         
@@ -339 +380 @@
         if sqrt_alpha<0.0:
             nq = 0
-        
+
+        # If we need to fit the background, add a term
+        if self.has_bck==True:
+            nfunc_0 = nfunc
+            nfunc += 1
+
         a = numpy.zeros([npts+nq, nfunc])
         b = numpy.zeros(npts+nq)
         err = numpy.zeros([nfunc, nfunc])
         
-        for j in range(nfunc):
-            for i in range(npts):
-                if self._accept_q(self.x[i]):
-                    a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
-                    
-            #TODO: refactor this: i_q should really be i_r
-            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):
-            if self._accept_q(self.x[i]):
-                b[i] = self.y[i]/self.err[i]
-            
+        # Construct the a matrix and b vector that represent the problem
+        self._get_matrix(nfunc, nq, a, b)
+             
+        # Perform the inversion (least square fit)
         c, chi2, rank, n = lstsq(a, b)
         # Sanity check
@@ -367 +402 @@
         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+nr):
-                    #if self._accept_q(self.x[i]):
-                    inv_cov[i][j] += at[i][k]*a[k][j]
+        # Get the covariance matrix, defined as inv_cov = a_transposed * a
+        self._get_invcov_matrix(nfunc, nr, a, inv_cov)
                     
         # 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):
-                if self._accept_q(self.x[i]):
-                    sum_sig += (a[i][j])**2
-            for i in range(nq):
-                sum_reg += (a[i+npts][j])**2
+        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a)
                     
         if math.fabs(self.alpha)>0:
@@ -403 +426 @@
             
         # Keep a copy of the last output
-        self.out = c
-        self.cov = err
-        
-        return c, err
+        if self.has_bck==False:
+            self.background = 0
+            self.out = c
+            self.cov = err
+        else:
+            self.background = c[0]
+            
+            err_0 = numpy.zeros([nfunc, nfunc])
+            c_0 = numpy.zeros(nfunc)
+            
+            for i in range(nfunc_0):
+                c_0[i] = c[i+1]
+                for j in range(nfunc_0):
+                    err_0[i][j] = err[i+1][j+1]
+                    
+            self.out = c_0
+            self.cov = err_0
+            
+        return self.out, self.cov
+        
+    def lstsq_bck(self, nfunc=5, nr=20):
+        #TODO: Allow for background by starting at n=0 (since the base function
+        # is zero for n=0).
+        import math
+        from scipy.linalg.basic import lstsq
+        
+        if self.is_valid()<0:
+            raise RuntimeError, "Invertor: invalid data; incompatible data lengths."
+        
+        self.nfunc = nfunc
+        # a -- An M x N matrix.
+        # b -- An M x nrhs matrix or M vector.
+        npts = len(self.x)
+        nq   = nr
+        sqrt_alpha = math.sqrt(math.fabs(self.alpha))
+        if sqrt_alpha<0.0:
+            nq = 0
+        
+        err_0 = numpy.zeros([nfunc, nfunc])
+        c_0 = numpy.zeros(nfunc)
+        nfunc_0 = nfunc
+        nfunc += 1
+        
+        a = numpy.zeros([npts+nq, nfunc])
+        b = numpy.zeros(npts+nq)
+        err = numpy.zeros([nfunc, nfunc])
+        
+        # Construct the a matrix and b vector that represent the problem
+        self._get_matrix(nfunc, nq, a, b)
+            
+        c, chi2, rank, n = lstsq(a, b)
+        # Sanity check
+        try:
+            float(chi2)
+        except:
+            chi2 = -1.0
+        self.chi2 = chi2
+
+        inv_cov = numpy.zeros([nfunc,nfunc])
+        # Get the covariance matrix, defined as inv_cov = a_transposed * a
+        self._get_invcov_matrix(nfunc, nr, a, inv_cov)
+                    
+        # Compute the reg term size for the output
+        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a)
+        
+        if math.fabs(self.alpha)>0:
+            new_alpha = sum_sig/(sum_reg/self.alpha)
+        else:
+            new_alpha = 0.0
+        self.suggested_alpha = new_alpha
+        
+        try:
+            cov = numpy.linalg.pinv(inv_cov)
+            err = math.fabs(chi2/float(npts-nfunc)) * cov
+        except:
+            # We were not able to estimate the errors,
+            # returns an empty covariance matrix
+            print "lstsq:", sys.exc_value
+            print chi2
+            pass
+            
+        # Keep a copy of the last output
+        
+        print "BACKGROUND =", c[0]
+        self.background = c[0]
+        
+        for i in range(nfunc_0):
+            c_0[i] = c[i+1]
+            for j in range(nfunc_0):
+                err_0[i][j] = err[i+1][j+1]
+                
+        self.out = c_0
+        self.cov = err_0
+        
+        return c_0, err_0
         
     def estimate_alpha(self, nfunc):
@@ -433 +547 @@
                  
             # Perform inversion to find the largest alpha
-            out, cov = pr.lstsq(nfunc)
+            out, cov = pr.invert(nfunc)
             elapsed = time.time()-starttime
             initial_alpha = pr.alpha
@@ -440 +554 @@
             # Try the inversion with the estimated alpha
             pr.alpha = pr.suggested_alpha
-            out, cov = pr.lstsq(nfunc)
+            out, cov = pr.invert(nfunc)
     
             npeaks = pr.get_peaks(out)
@@ -458 +572 @@
                 for i in range(10):
                     pr.alpha = (0.33)**(i+1)*alpha
-                    out, cov = pr.lstsq(nfunc)
+                    out, cov = pr.invert(nfunc)
                     
                     peaks = pr.get_peaks(out)

pr_inversion/test/test_output.txt

@@ -5 +5 @@
 #elapsed=0
 #alpha_estimate=4.58991e-006
-#C_0=217.876478953+-0.0394946442786
-#C_1=0.30164617344+-0.407412922728
-#C_2=7.15694574584+-2.73726481828
-#C_3=-0.79690812244+-5.7183647083
-#C_4=0.948333847403+-16.2086432642
-#C_5=-0.751612279429+-28.2108893047
-#C_6=-0.0796650783501+-52.0568306543
-#C_7=-0.518109559543+-88.6068891309
-#C_8=-0.129910831984+-141.7178308
-#C_9=-0.258134624641+-215.767430098
+#C_0=217.876478953+-2522.42314999
+#C_1=0.30164617344+-226.644556603
+#C_2=7.15694574584+-43.4919208221
+#C_3=-0.79690812244+-25.8031634616
+#C_4=0.948333847403+-11.3794625443
+#C_5=-0.751612279429+-8.64083336005
+#C_6=-0.0796650783501+-5.73820483268
+#C_7=-0.518109559543+-3.56660274557
+#C_8=-0.129910831984+-1.99841316045
+#C_9=-0.258134624641+-0.756509301101
 <r>  <Pr>  <dPr>
 0  0  0
-1.6  22.981  20.5137
-3.2  91.986  78.6612
-4.8  207.187  164.741
-6.4  368.819  264.129
-8  577.117  359.427
-9.6  832.224  433.162
-11.2  1134.11  470.89
-12.8  1482.46  464.814
-14.4  1876.63  419.086
-16  2315.57  360.693
-17.6  2797.81  353.245
-19.2  3321.46  448.357
-20.8  3884.25  611.481
-22.4  4483.6  786.464
-24  5116.72  934.559
-25.6  5780.7  1030.78
-27.2  6472.66  1062.21
-28.8  7189.82  1030.28
-30.4  7929.59  955.681
-32  8689.7  883.97
-33.6  9468.17  878.871
-35.2  10263.4  977.553
-36.8  11073.9  1152
-38.4  11898.7  1342.59
-40  12737  1498.3
-41.6  13587.7  1587.08
-43.2  14450.2  1596.95
-44.8  15323.5  1537.99
-46.4  16206.5  1445.58
-48  17098  1379.39
-49.6  17996.3  1401.86
-51.2  18899.8  1532.06
-52.8  19806.5  1728.86
-54.4  20714.3  1927.37
-56  21620.8  2073.39
-57.6  22524  2135.21
-59.2  23421.6  2106.97
-60.8  24311.7  2011.28
-62.4  25192.3  1900.75
-64  26061.8  1848.67
-65.6  26918.7  1912.83
-67.2  27761.7  2089.96
-68.8  28589.6  2319.48
-70.4  29401.3  2527.71
-72  30195.6  2658.76
-73.6  30971.4  2685.17
-75.2  31727.2  2611.68
-76.8  32461.5  2478.05
-78.4  33172.7  2357.19
-80  33858.8  2333.81
-81.6  34517.8  2453.43
-83.2  35147.6  2684.61
-84.8  35745.9  2945.13
-86.4  36310.5  3152.18
-88  36839.5  3249.59
-89.6  37330.9  3217.04
-91.2  37783.1  3074.61
-92.8  38194.6  2886.09
-94.4  38564.3  2752.58
-96  38891.1  2773.03
-97.6  39174.1  2973.02
-99.2  39412.6  3281.62
-100.8  39605.8  3587.59
-102.4  39752.6  3795.95
-104  39851.9  3850.56
-105.6  39902.4  3741.32
-107.2  39902.5  3510.03
-108.8  39850.1  3255.58
-110.4  39743.2  3120.25
-112  39579.6  3218.61
-113.6  39356.7  3539.26
-115.2  39072.3  3954.92
-116.8  38724.1  4317.73
-118.4  38310.3  4517.08
-120  37829  4494.86
-121.6  37278.9  4251.44
-123.2  36659  3856.58
-124.8  35968.4  3465.78
-126.4  35206.6  3306.22
-128  34373  3535.97
-129.6  33466.9  4074.35
-131.2  32487.5  4693.64
-132.8  31433.4  5186.14
-134.4  30303.1  5411.6
-136  29094.1  5295.75
-137.6  27803.8  4830.98
-139.2  26428.9  4092.02
-140.8  24965.9  3287.78
-142.4  23411.1  2858.06
-144  21761.3  3246.41
-145.6  20013.4  4237.67
-147.2  18165.4  5345.59
-148.8  16216.4  6256.27
-150.4  14166.9  6786.26
-152  12019  6829.04
-153.6  9776.61  6338.91
-155.2  7445.34  5326.31
-156.8  5032.57  3854
-158.4  2547.22  2030.4
+1.6  22.981  7.48836
+3.2  91.986  29.737
+4.8  207.187  66.1302
+6.4  368.819  115.749
+8  577.117  177.508
+9.6  832.224  250.315
+11.2  1134.11  333.205
+12.8  1482.46  425.442
+14.4  1876.63  526.549
+16  2315.57  636.272
+17.6  2797.81  754.499
+19.2  3321.46  881.145
+20.8  3884.25  1016.06
+22.4  4483.6  1158.97
+24  5116.72  1309.47
+25.6  5780.7  1467.05
+27.2  6472.66  1631.18
+28.8  7189.82  1801.33
+30.4  7929.59  1977.02
+32  8689.7  2157.88
+33.6  9468.17  2343.55
+35.2  10263.4  2533.72
+36.8  11073.9  2728.06
+38.4  11898.7  2926.2
+40  12737  3127.75
+41.6  13587.7  3332.27
+43.2  14450.2  3539.34
+44.8  15323.5  3748.54
+46.4  16206.5  3959.49
+48  17098  4171.83
+49.6  17996.3  4385.18
+51.2  18899.8  4599.19
+52.8  19806.5  4813.46
+54.4  20714.3  5027.56
+56  21620.8  5241.07
+57.6  22524  5453.53
+59.2  23421.6  5664.53
+60.8  24311.7  5873.68
+62.4  25192.3  6080.65
+64  26061.8  6285.13
+65.6  26918.7  6486.85
+67.2  27761.7  6685.53
+68.8  28589.6  6880.91
+70.4  29401.3  7072.73
+72  30195.6  7260.72
+73.6  30971.4  7444.63
+75.2  31727.2  7624.23
+76.8  32461.5  7799.3
+78.4  33172.7  7969.61
+80  33858.8  8134.91
+81.6  34517.8  8294.9
+83.2  35147.6  8449.24
+84.8  35745.9  8597.54
+86.4  36310.5  8739.34
+88  36839.5  8874.2
+89.6  37330.9  9001.64
+91.2  37783.1  9121.2
+92.8  38194.6  9232.39
+94.4  38564.3  9334.72
+96  38891.1  9427.68
+97.6  39174.1  9510.74
+99.2  39412.6  9583.35
+100.8  39605.8  9644.98
+102.4  39752.6  9695.14
+104  39851.9  9733.38
+105.6  39902.4  9759.36
+107.2  39902.5  9772.74
+108.8  39850.1  9773.24
+110.4  39743.2  9760.52
+112  39579.6  9734.22
+113.6  39356.7  9693.89
+115.2  39072.3  9639.07
+116.8  38724.1  9569.28
+118.4  38310.3  9484.1
+120  37829  9383.18
+121.6  37278.9  9266.24
+123.2  36659  9133
+124.8  35968.4  8983.13
+126.4  35206.6  8816.15
+128  34373  8631.45
+129.6  33466.9  8428.29
+131.2  32487.5  8205.96
+132.8  31433.4  7963.92
+134.4  30303.1  7701.95
+136  29094.1  7420.23
+137.6  27803.8  7119.29
+139.2  26428.9  6799.69
+140.8  24965.9  6461.63
+142.4  23411.1  6104.48
+144  21761.3  5726.35
+145.6  20013.4  5323.88
+147.2  18165.4  4892.46
+148.8  16216.4  4426.74
+150.4  14166.9  3921.46
+152  12019  3372.52
+153.6  9776.61  2777.82
+155.2  7445.34  2138.05
+156.8  5032.57  1456.94
+158.4  2547.22  741.175

pr_inversion/test/utest_invertor.py

@@ -63 +63 @@
             self.x_in[i] = 1.0*(i+1)
 
+    def test_has_bck_flag(self):
+        """
+            Tests the has_bck flag operations
+        """
+        self.assertEqual(self.invertor.has_bck, False)
+        self.invertor.has_bck=True
+        self.assertEqual(self.invertor.has_bck, True)
+        def doit_float():
+            self.invertor.has_bck  = 2.0
+        def doit_str():
+            self.invertor.has_bck  = 'a'
+
+        self.assertRaises(ValueError, doit_float)
+        self.assertRaises(ValueError, doit_str)
+        
 
     def testset_dmax(self):
@@ -177 +192 @@
         # Perform inversion
         out, cov = self.invertor.invert_optimize(10)
+        #out, cov = self.invertor.invert(10)
         # This is a very specific case
         # We should make sure it always passes