Changeset 43c0a8e in sasview for pr_inversion/c_extensions

Timestamp:

May 23, 2008 10:55:39 AM (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:

dfb58f8

Parents:

eb06cbe

Message:

Add a couple of figures of merit

Location:

pr_inversion/c_extensions

Files:

• Cinvertor.c (modified) (3 diffs)
• invertor.c (modified) (3 diffs)
• invertor.h (modified) (1 diff)

pr_inversion/c_extensions/Cinvertor.c

@@ -492 +492 @@
         PyObject *err_obj;
         Py_ssize_t npars2;
-          
+        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);     
@@ -537 +538 @@
 
         return Py_BuildValue("i", count );      
+        
+}
+
+static PyObject * get_positive(Cinvertor *self, PyObject *args) {
+        double *pars;
+        PyObject *data_obj;
+        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 );   
+        
+}
+
+static PyObject * get_pos_err(Cinvertor *self, PyObject *args) {
+        double *pars;
+        double *pars_err;
+        PyObject *data_obj;
+        PyObject *err_obj;
+        Py_ssize_t npars;
+        Py_ssize_t npars2;
+        double fraction;
+        
+        if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL;
+        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 );   
         
 }
@@ -567 +602 @@
                    {"oscillations", (PyCFunction)oscillations, METH_VARARGS, ""},
                    {"get_peaks", (PyCFunction)get_peaks, METH_VARARGS, ""},
+                   {"get_positive", (PyCFunction)get_positive, METH_VARARGS, ""},
+                   {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, ""},
    
    {NULL}

pr_inversion/c_extensions/invertor.c

@@ -100 +100 @@
         func_value = ortho(d_max, i+1, r);
         sum += pars[i] * func_value;
-        sum_err += err[i]*err[i]*func_value*func_value;
+        //sum_err += err[i]*err[i]*func_value*func_value;
+        sum_err += err[i*n_c+i]*func_value*func_value;
     }
     *pr_value = sum;
@@ -115 +116 @@
 double dprdr(double *pars, double d_max, int n_c, double r) {
     double sum = 0.0; 
-        int i;
+        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));
@@ -180 +181 @@
 }
 
+/**
+ * Get the fraction of the integral of P(r) over the whole range
+ * of r that is above zero.
+ * A valid P(r) is define as being positive for all r.
+ */
+double positive_integral(double *pars, double d_max, int n_c, int nslice) {
+    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;
+        value = pr(pars, d_max, n_c, r);
+        if (value>0.0) sum_pos += value;
+        sum += value;
+    }
+    return sum_pos/sum;
+}
+
+/**
+ * Get the fraction of the integral of P(r) over the whole range
+ * of r that is at least one sigma above zero.
+ */
+double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) {
+    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;
+        pr_err(pars, err, d_max, n_c, r, &pr_val, &pr_val_err);
+        if (pr_val>pr_val_err) sum_pos += pr_val;
+        sum += pr_val;
+        
+
+    }
+    return sum_pos/sum;
+}

pr_inversion/c_extensions/invertor.h

@@ -33 +33 @@
                 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);
 
 #endif