Changeset 43c0a8e in sasview for pr_inversion

Timestamp:

May 23, 2008 8: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

Files:

• c_extensions/Cinvertor.c (modified) (3 diffs)
• c_extensions/invertor.c (modified) (3 diffs)
• c_extensions/invertor.h (modified) (1 diff)
• invertor.py (modified) (3 diffs)
• test/utest_invertor.py (modified) (2 diffs)

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

pr_inversion/invertor.py

@@ -6 +6 @@
     """
         Provide general online help text
+        Future work: extend this function to allow topic selection
     """
     info_txt  = "The inversion approach is based on Moore, J. Appl. Cryst. (1980) 13, 168-175.\n\n"
@@ -31 +32 @@
         
         TODO: explain the maths
+        
+        
+        Methods inherited from Cinvertor:
+        - get_peaks: returns the number of P(r) peaks
     """
     ## Chisqr of the last computation
@@ -189 +194 @@
     
     def pr_err(self, c, c_cov, r):
-        import math
-        c_err = numpy.zeros(len(c))
-        for i in range(len(c)):
-            try:
-                c_err[i] = math.sqrt(math.fabs(c_cov[i][i]))
-            except:
-                import sys
-                print sys.exc_value
-                print "oups", c_cov[i][i]
-                c_err[i] = c[i]
+        #import math
+        #c_err = numpy.zeros(len(c))
+        #for i in range(len(c)):
+        #    try:
+        #        c_err[i] = math.sqrt(math.fabs(c_cov[i][i]))
+        #    except:
+        #        import sys
+        #        print sys.exc_value
+        #        print "oups", c_cov[i][i]
+        #        c_err[i] = c[i]
 
-        return self.get_pr_err(c, c_err, r)
+        return self.get_pr_err(c, c_cov, r)
        
     def _accept_q(self, q):

pr_inversion/test/utest_invertor.py

@@ -11 +11 @@
 from sans.pr.invertor import Invertor
         
+class TestFiguresOfMerit(unittest.TestCase):
+            
+    def setUp(self):
+        self.invertor = Invertor()
+        self.invertor.d_max = 100.0
+        
+        # Test array
+        self.ntest = 5
+        self.x_in = numpy.ones(self.ntest)
+        for i in range(self.ntest):
+            self.x_in[i] = 1.0*(i+1)
+       
+        x, y, err = load("sphere_80.txt")
+
+        # Choose the right d_max...
+        self.invertor.d_max = 160.0
+        # Set a small alpha
+        self.invertor.alpha = .0007
+        # Set data
+        self.invertor.x   = x
+        self.invertor.y   = y
+        self.invertor.err = err
+        # Perform inversion
+        #out, cov = self.invertor.invert(10)
+        
+        self.out, self.cov = self.invertor.lstsq(10)
+
+    def test_positive(self):
+        """
+            Test whether P(r) is positive
+        """
+        self.assertEqual(self.invertor.get_positive(self.out), 1)
+        
+    def test_positive_err(self):
+        """
+            Test whether P(r) is at least 1 sigma greater than zero
+            for all r-values
+        """
+        self.assertTrue(self.invertor.get_pos_err(self.out, self.cov)>0.9)
+       
 class TestBasicComponent(unittest.TestCase):
     
@@ -236 +276 @@
             print "chi2(P(r)) =", chi2/51.0
             raise
+        
+        # Test the number of peaks
+        self.assertEqual(self.invertor.get_peaks(out), 1)
             
     def test_q_zero(self):