Changeset eca05c8 in sasview for pr_inversion/c_extensions

Timestamp:

Apr 30, 2008 4:29:36 PM (17 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:

34ae302

Parents:

47f695c9

Message:

Added P(r) fit test

Location:

pr_inversion/c_extensions

Files:

• Cinvertor.c (modified) (8 diffs)
• invertor.c (modified) (3 diffs)
• invertor.h (modified) (2 diffs)

pr_inversion/c_extensions/Cinvertor.c

@@ -222 +222 @@
 }
 
+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);  
+}
+
+/**
+ * Gets the maximum distance
+ */
+static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("d", self->params.alpha);  
+}
+
 /**
  * Gets the number of x points
@@ -246 +261 @@
 /**
  * Function to call to evaluate the residuals
- * @param args: input q or [q,phi]
- * @return: function value
+ * @param args: input parameters
+ * @return: list of residuals
  */
 static PyObject * residuals(Cinvertor *self, PyObject *args) {
@@ -257 +272 @@
         double residual, diff;
         // Regularization factor
-        double lmult = 0.0;
+        double regterm = 0.0;
+        double tmp = 0.0;
         
         PyObject *data_obj;
@@ -266 +282 @@
         OUTVECTOR(data_obj,pars,npars);
                 
-        //printf("npars: %i", npars);
     // PyList of residuals
         // Should create this list only once and refill it
     residuals = PyList_New(self->params.npoints);
+
+    regterm = reg_term(pars, self->params.d_max, npars);
     
     for(i=0; i<self->params.npoints; i++) {
-        //printf("\n%i: ", i);
         diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
         residual = diff*diff / (self->params.err[i]*self->params.err[i]);
-        //printf("%i %g\n", i, residual);
+        tmp = residual;
+        
+        // regularization term
+        residual += self->params.alpha * regterm;
+        
         if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
             PyErr_SetString(CinvertorError, 
@@ -286 +306 @@
         return residuals;
 }
+/**
+ * Function to call to evaluate the residuals
+ * for P(r) minimization (for testing purposes)
+ * @param args: input parameters
+ * @return: list of residuals
+ */
+static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
+        double *pars;
+        PyObject* residuals;
+        PyObject* temp;
+        double *res;
+        int i;
+        double residual, diff;
+        // Regularization factor
+        double regterm = 0.0;
+        double tmp = 0.0;
+        
+        PyObject *data_obj;
+        Py_ssize_t npars;
+          
+        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);
+
+    regterm = reg_term(pars, self->params.d_max, npars);
+
+    
+    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, 
+                "Cinvertor.residuals: error setting residual.");
+                return NULL;
+        };
+                
+    }
+    
+        return residuals;
+}
 
 /**
@@ -323 +391 @@
 }
 
+/**
+ * Function to call to evaluate P(r) with errors
+ * @param args: c-parameters and r
+ * @return: P(r)
+ */
+static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
+        double *pars;
+        double *pars_err;
+        double pr_err_value;
+        double r, pr_value;
+        PyObject *data_obj;
+        Py_ssize_t npars;
+        PyObject *err_obj;
+        Py_ssize_t npars2;
+          
+        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
+        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);     
+}
+
 
 static PyMethodDef Cinvertor_methods[] = {
                    {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
+                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"},
                    {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
                    {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
@@ -334 +426 @@
                    {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
                    {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
+                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
+                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
                    {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
                    {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
@@ -339 +433 @@
                    {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
                    {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
+                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
                    {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
    

pr_inversion/c_extensions/invertor.c

@@ -40 +40 @@
 double ortho(double d_max, int n, double r) {
         return 2.0*r*sin(pi*n*r/d_max);
-        
 }
 
@@ -76 +75 @@
  */
 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++) {
@@ -84 +83 @@
 }
 
+/**
+ * P(r) calculated from the expansion, with errors
+ */
+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_err = 0.0;
+    double func_value;
+        int i;
+    for (i=0; i<n_c; i++) {
+        func_value = ortho(d_max, i+1, r);
+        sum += pars[i] * func_value;
+        sum_err += err[i]*err[i]*func_value*func_value;
+    }
+    *pr_value = sum;
+    if (sum_err>0) {
+        *pr_value_err = sqrt(sum_err);
+    } else {
+        *pr_value_err = sum;
+    }
+} 
 
+/**
+ * dP(r)/dr calculated from the expansion.
+ */
+double dprdr(double *pars, double d_max, int n_c, double r) {
+    double sum = 0.0; 
+        int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max));
+    }
+    return sum;
+}
+
+/**
+ * regularization term calculated from the expansion.
+ */
+double reg_term(double *pars, double d_max, int n_c) {
+    double sum = 0.0; 
+    double r;
+    double deriv;
+        int i;
+    for (i=0; i<25; i++) {
+        r = d_max/25.0*i;
+        deriv = dprdr(pars, d_max, n_c, r);
+        sum += deriv*deriv;
+    }
+    return sum/25.0*d_max;
+}
+
+

pr_inversion/c_extensions/invertor.h

@@ -10 +10 @@
     int npoints;    
     int ny;    
-    int nerr;    
-} Invertor_params;
+    int nerr;  
+    double alpha;
+} Invertor_params; 
 
 void invertor_dealloc(Invertor_params *pars);
@@ -24 +25 @@
 double iq(double *c, double d_max, int n_c, double q);
 double pr(double *c, double d_max, int n_c, double r);
+double dprdr(double *pars, double d_max, int n_c, double r);
+double reg_term(double *pars, 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);
 
 #endif