Changeset 9a23253e in sasview for pr_inversion/c_extensions

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/c_extensions

Files:

• Cinvertor.c (modified) (43 diffs)
• invertor.c (modified) (15 diffs)
• invertor.h (modified) (3 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