Changeset f71287f4 in sasview

Timestamp:

May 9, 2008 6:09:10 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:

634f1cf

Parents:

32dffae4

Message:

Allow user to set q_min/q_max

Location:

pr_inversion

Files:

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

pr_inversion/c_extensions/Cinvertor.c

@@ -266 +266 @@
         return Py_BuildValue("d", self->params.d_max);  
 }
+
+/**
+ * Sets the minimum q
+ */
+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);  
+}
+
+/**
+ * Gets the minimum q
+ */
+static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("d", self->params.q_min);  
+}
+
+
+/**
+ * Sets the maximum q
+ */
+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);  
+}
+
+/**
+ * Gets the maximum q
+ */
+static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
+        return Py_BuildValue("d", self->params.q_max);  
+}
+
 
 static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
@@ -513 +551 @@
                    {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
                    {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
+                   {"set_qmin", (PyCFunction)set_qmin, METH_VARARGS, ""},
+                   {"get_qmin", (PyCFunction)get_qmin, METH_VARARGS, ""},
+                   {"set_qmax", (PyCFunction)set_qmax, METH_VARARGS, ""},
+                   {"get_qmax", (PyCFunction)get_qmax, METH_VARARGS, ""},
                    {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
                    {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},

pr_inversion/c_extensions/invertor.c

@@ -18 +18 @@
 void invertor_init(Invertor_params *pars) {
         pars->d_max = 180;
+        pars->q_min = -1.0;
+        pars->q_max = -1.0;
 }
 

pr_inversion/c_extensions/invertor.h

@@ -12 +12 @@
     int nerr;  
     double alpha;
+    double q_min;
+    double q_max;
 } Invertor_params; 
 

pr_inversion/invertor.py

@@ -1 +1 @@
 from sans.pr.core.pr_inversion import Cinvertor
 import numpy
+import sys
 
 class Invertor(Cinvertor):
@@ -10 +11 @@
     ## Alpha to get the reg term the same size as the signal
     suggested_alpha = 0
+    ## Last number of base functions used
+    nfunc = 0
+    ## Last output values
+    out = None
+    ## Last errors on output values
+    cov = None
     
     def __init__(self):
@@ -30 +37 @@
         elif name=='d_max':
             return self.set_dmax(value)
+        elif name=='q_min':
+            if value==None:
+                return self.set_qmin(-1.0)
+            return self.set_qmin(value)
+        elif name=='q_max':
+            if value==None:
+                return self.set_qmax(-1.0)
+            return self.set_qmax(value)
         elif name=='alpha':
             return self.set_alpha(value)
@@ -56 +71 @@
         elif name=='d_max':
             return self.get_dmax()
+        elif name=='q_min':
+            qmin = self.get_qmin()
+            if qmin<0:
+                return None
+            return qmin
+        elif name=='q_max':
+            qmax = self.get_qmax()
+            if qmax<0:
+                return None
+            return qmax
         elif name=='alpha':
             return self.get_alpha()
@@ -71 +96 @@
         invertor.alpha   = self.alpha
         invertor.d_max   = self.d_max
+        invertor.q_min   = self.q_min
+        invertor.q_max   = self.q_max
         
         invertor.x = self.x
@@ -85 +112 @@
         import time
         
+        self.nfunc = nfunc
         # First, check that the current data is valid
         if self.is_valid()<=0:
@@ -147 +175 @@
         return self.get_pr_err(c, c_err, r)
        
+    def _accept_q(self, q):
+        """
+            Check q-value against user-defined range
+        """
+        if not self.q_min==None and q<self.q_min:
+            return False
+        if not self.q_max==None and q>self.q_max:
+            return False
+        return True
+       
     def lstsq(self, nfunc=5):
+        #TODO: do this on the C side
+        #
+        # To make sure an array is contiguous:
+        # blah = numpy.ascontiguousarray(blah_original)
+        # ... before passing it to C
+        
         import math
         from scipy.linalg.basic import lstsq
         
+        self.nfunc = nfunc
         # a -- An M x N matrix.
         # b -- An M x nrhs matrix or M vector.
@@ -163 +208 @@
         for j in range(nfunc):
             for i in range(npts):
-                a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
+                if self._accept_q(self.x[i]):
+                    a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
             for i_q in range(nq):
                 r = self.d_max/nq*i_q
@@ -170 +216 @@
         
         for i in range(npts):
-            b[i] = self.y[i]/self.err[i]
+            if self._accept_q(self.x[i]):
+                b[i] = self.y[i]/self.err[i]
             
         c, chi2, rank, n = lstsq(a, b)
@@ -181 +228 @@
                 inv_cov[i][j] = 0.0
                 for k in range(npts):
-                    inv_cov[i][j] = at[i][k]*a[k][j]
+                    if self._accept_q(self.x[i]):
+                        inv_cov[i][j] = at[i][k]*a[k][j]
                     
         # Compute the reg term size for the output
@@ -188 +236 @@
         for j in range(nfunc):
             for i in range(npts):
-                sum_sig += (a[i][j])**2
+                if self._accept_q(self.x[i]):
+                    sum_sig += (a[i][j])**2
             for i in range(nq):
-                sum_reg += (a[i_q+npts][j])**2
+                sum_reg += (a[i+npts][j])**2
                     
         if math.fabs(self.alpha)>0:
@@ -203 +252 @@
             print "Error estimating uncertainties"
             
+        # Keep a copy of the last output
+        self.out = c
+        self.cov = err
         
         return c, err
@@ -248 +300 @@
         return c, err
         
-        
-        
-        
+    def estimate_alpha(self, nfunc):
+        """
+            Returns a reasonable guess for the
+            regularization constant alpha
+            
+            @return: alpha, message, elapsed
+            
+            where alpha is the estimate for alpha,
+            message is a message for the user,
+            elapsed is the computation time
+        """
+        import time
+        try:            
+            pr = self.clone()
+            
+            # T_0 for computation time
+            starttime = time.time()
+            
+            # If the current alpha is zero, try
+            # another value
+            if pr.alpha<=0:
+                pr.alpha = 0.0001
+                 
+            # Perform inversion to find the largest alpha
+            out, cov = pr.lstsq(nfunc)
+            elapsed = time.time()-starttime
+            initial_alpha = pr.alpha
+            initial_peaks = pr.get_peaks(out)
+    
+            # Try the inversion with the estimated alpha
+            pr.alpha = pr.suggested_alpha
+            out, cov = pr.lstsq(nfunc)
+    
+            npeaks = pr.get_peaks(out)
+            # if more than one peak to start with
+            # just return the estimate
+            if npeaks>1:
+                message = "Your P(r) is not smooth, please check your inversion parameters"
+                return pr.suggested_alpha, message, elapsed
+            else:
+                
+                # Look at smaller values
+                # We assume that for the suggested alpha, we have 1 peak
+                # if not, send a message to change parameters
+                alpha = pr.suggested_alpha
+                best_alpha = pr.suggested_alpha
+                found = False
+                for i in range(10):
+                    pr.alpha = (0.33)**(i+1)*alpha
+                    out, cov = pr.lstsq(nfunc)
+                    
+                    peaks = pr.get_peaks(out)
+                    print pr.alpha, peaks
+                    if peaks>1:
+                        found = True
+                        break
+                    best_alpha = pr.alpha
+                    
+                # If we didn't find a turning point for alpha and
+                # the initial alpha already had only one peak,
+                # just return that
+                if not found and initial_peaks==1 and initial_alpha<best_alpha:
+                    best_alpha = initial_alpha
+                    
+                # Check whether the size makes sense
+                message=''
+                
+                if not found:
+                    message = "None"
+                elif best_alpha>=0.5*pr.suggested_alpha:
+                    # best alpha is too big, return a 
+                    # reasonable value
+                    message  = "The estimated alpha for your system is too large. "
+                    message += "Try increasing your maximum distance."
+                
+                return best_alpha, message, elapsed
+    
+        except:
+            message = "Invertor.estimate_alpha: %s" % sys.exc_value
+            return 0, message, elapsed
+    
+        
+    def to_file(self, path, npts=100):
+        """
+            Save the state to a file that will be readable
+            by SliceView.
+            @param path: path of the file to write
+            @param npts: number of P(r) points to be written
+        """
+        import pylab
+        
+        file = open(path, 'w')
+        file.write("#d_max=%g\n" % self.d_max)
+        file.write("#nfunc=%g\n" % self.nfunc)
+        file.write("#alpha=%g\n" % self.alpha)
+        file.write("#chi2=%g\n" % self.chi2)
+        file.write("#elapsed=%g\n" % self.elapsed)
+        file.write("#alpha_estimate=%g\n" % self.suggested_alpha)
+        if not self.out==None:
+            if len(self.out)==len(self.cov):
+                for i in range(len(self.out)):
+                    file.write("#C_%i=%s+-%s\n" % (i, str(self.out[i]), str(self.cov[i][i])))
+        file.write("<r>  <Pr>  <dPr>\n")
+        r = pylab.arange(0.0, self.d_max, self.d_max/npts)
+        
+        for r_i in r:
+            (value, err) = self.pr_err(self.out, self.cov, r_i)
+            file.write("%g  %g  %g\n" % (r_i, value, err))
+    
+        file.close()
+        
+    def from_file(self, path):
+        """
+            Load the state of the Invertor from a file,
+            to be able to generate P(r) from a set of
+            parameters.
+            @param path: path of the file to load
+        """
+        import os
+        import re
+        if os.path.isfile(path):
+            try:
+                fd = open(path, 'r')
+                
+                buff    = fd.read()
+                lines   = buff.split('\n')
+                for line in lines:
+                    if line.startswith('#d_max='):
+                        toks = line.split('=')
+                        self.d_max = float(toks[1])
+                    elif line.startswith('#nfunc='):
+                        toks = line.split('=')
+                        self.nfunc = int(toks[1])
+                        self.out = numpy.zeros(self.nfunc)
+                        self.cov = numpy.zeros([self.nfunc, self.nfunc])
+                    elif line.startswith('#alpha='):
+                        toks = line.split('=')
+                        self.alpha = float(toks[1])
+                    elif line.startswith('#chi2='):
+                        toks = line.split('=')
+                        self.chi2 = float(toks[1])
+                    elif line.startswith('#elapsed='):
+                        toks = line.split('=')
+                        self.elapsed = float(toks[1])
+                    elif line.startswith('#alpha_estimate='):
+                        toks = line.split('=')
+                        self.suggested_alpha = float(toks[1])
+            
+                    # Now read in the parameters
+                    elif line.startswith('#C_'):
+                        toks = line.split('=')
+                        p = re.compile('#C_([0-9]+)')
+                        m = p.search(toks[0])
+                        toks2 = toks[1].split('+-')
+                        i = int(m.group(1))
+                        self.out[i] = float(toks2[0])
+                        
+                        self.cov[i][i] = float(toks2[1])                        
+            
+            except:
+                raise RuntimeError, "Invertor.from_file: corrupted file\n%s" % sys.exc_value
+        else:
+            raise RuntimeError, "Invertor.from_file: '%s' is not a file" % str(path) 
+        
+        
+    
     
 if __name__ == "__main__":

pr_inversion/test/utest_invertor.py

@@ -333 +333 @@
         for i in range(len(self.x_in)):
             self.assertEqual(self.x_in[i], clone.x[i])
+        
+    def test_save(self):
+        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.lstsq(10)
+        
+        # Save
+        self.invertor.to_file("test_output.txt")
+    
+    def test_load(self):
+        self.invertor.from_file("test_output.txt")
+        self.assertEqual(self.invertor.d_max, 160.0)
+        self.assertEqual(self.invertor.alpha, 0.0007)
+        self.assertEqual(self.invertor.chi2, 16654.1)
+        self.assertAlmostEqual(self.invertor.pr(self.invertor.out, 10.0), 8948.22689927, 4)
+        
+    def test_qmin(self):
+        self.invertor.q_min = 1.0
+        print self.invertor.q_min
+        self.assertEqual(self.invertor.q_min, 1.0)
+        
+        self.invertor.q_min = None
+        self.assertEqual(self.invertor.q_min, None)
+        
+                         
+    def test_qmax(self):
+        self.invertor.q_max = 1.0
+        self.assertEqual(self.invertor.q_max, 1.0)
+       
+        self.invertor.q_max = None
+        self.assertEqual(self.invertor.q_max, None)
+
         
 def pr_theory(r, R):