Changeset ffca8f2 in sasview for pr_inversion

Timestamp:

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

a8b6364

Parents:

896abb3

Message:

More docs and a little cleaning up

Location:

pr_inversion

Files:

• __init__.py (modified) (2 diffs)
• invertor.py (modified) (10 diffs)
• test/utest_invertor.py (modified) (2 diffs)

pr_inversion/__init__.py

@@ -24 +24 @@
 # [5] S. Hansen and J. Skov Pedersen, J.Appl. Cryst (1991) 24, 541-548.
 #
+## \subsection class Class Diagram:
+# The following shows a partial class diagram with the main attributes and methods of the invertor.
+#
+# \image html architecture.png
 #
 # \section install_sec Installation
@@ -77 +81 @@
 # matrix for those coefficients, respectively.
 #
+# To get P(r):
+# \verbatim
+#    r = 10.0
+#    pr = invertor.pr(c_out, r)
+# \endverbatim
+# Alternatively, one can get P(r) with the error on P(r):
+# \verbatim
+#    r = 10.0
+#    pr, dpr = invertor.pr_err(c_out, c_cov, r)
+# \endverbatim
+#
+# To get the output I(q) from the set of coefficients found:
+# \verbatim
+#    q = 0.001
+#    iq = invertor.iq(c_out, q)
+# \endverbatim
+#
 # Examples are available as unit tests under sans.pr_inversion.test.
 #

pr_inversion/invertor.py

@@ -35 +35 @@
         Invertor class to perform P(r) inversion
         
-        TODO: explain the maths
-        
+        The problem is solved by posing the problem as  Ax = b,
+        where x is the set of coefficients we are looking for.
+        
+        Npts is the number of points.
+        
+        In the following i refers to the ith base function coefficient.
+        The matrix has its entries j in its first Npts rows set to
+            A[i][j] = (Fourier transformed base function for point j) 
+            
+        We them choose a number of r-points, n_r, to evaluate the second
+        derivative of P(r) at. This is used as our regularization term.
+        For a vector r of length n_r, the following n_r rows are set to
+            A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])
+            
+        The vector b has its first Npts entries set to
+            b[j] = (I(q) observed for point j)
+            
+        The following n_r entries are set to zero.
+        
+        The result is found by using scipy.linalg.basic.lstsq to invert
+        the matrix and find the coefficients x.
         
         Methods inherited from Cinvertor:
@@ -144 +163 @@
         return invertor
     
-    def invert(self, nfunc=5):
+    def invert(self, nfunc=10, nr=20):
         """
             Perform inversion to P(r)
-        """
+            
+            The problem is solved by posing the problem as  Ax = b,
+            where x is the set of coefficients we are looking for.
+            
+            Npts is the number of points.
+            
+            In the following i refers to the ith base function coefficient.
+            The matrix has its entries j in its first Npts rows set to
+                A[i][j] = (Fourier transformed base function for point j) 
+                
+            We them choose a number of r-points, n_r, to evaluate the second
+            derivative of P(r) at. This is used as our regularization term.
+            For a vector r of length n_r, the following n_r rows are set to
+                A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])
+                
+            The vector b has its first Npts entries set to
+                b[j] = (I(q) observed for point j)
+                
+            The following n_r entries are set to zero.
+            
+            The result is found by using scipy.linalg.basic.lstsq to invert
+            the matrix and find the coefficients x.
+            
+            @param nfunc: number of base functions to use.
+            @param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
+            @return: c_out, c_cov - the coefficients with covariance matrix 
+        """
+        #TODO: call the pyhton implementation for now. In the future, translate this to C.
+        return self.lstsq(nfunc, nr=nr)
+    
+    def invert_optimize(self, nfunc=10, nr=20):
+        """
+            Slower version of the P(r) inversion that uses scipy.optimize.leastsq.
+            
+            This probably produce more reliable results, but is much slower.
+            The minimization function is set to sum_i[ (I_obs(q_i) - I_theo(q_i))/err**2 ] + alpha * reg_term,
+            where the reg_term is given by Svergun: it is the integral of the square of the first derivative
+            of P(r), d(P(r))/dr, integrated over the full range of r.
+            
+            @param nfunc: number of base functions to use.
+            @param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
+            @return: c_out, c_cov - the coefficients with covariance matrix 
+        """
+        
         from scipy import optimize
         import time
@@ -175 +237 @@
     def pr_fit(self, nfunc=5):
         """
-            Perform inversion to P(r)
+            This is a direct fit to a given P(r). It assumes that the y data
+            is set to some P(r) distribution that we are trying to reproduce
+            with a set of base functions.
+            
+            This method is provided as a test. 
         """
         from scipy import optimize
@@ -209 +275 @@
             @param r: r-value to evaluate P(r) at
             @return: P(r)
-            
         """
         return self.get_pr_err(c, c_cov, r)
@@ -223 +288 @@
         return True
        
-    def lstsq(self, nfunc=5):
+    def lstsq(self, nfunc=5, nr=20):
         #TODO: do this on the C side
         #
@@ -230 +295 @@
         # ... before passing it to C
         """
-            TODO: Document this
+            The problem is solved by posing the problem as  Ax = b,
+            where x is the set of coefficients we are looking for.
+            
+            Npts is the number of points.
+            
+            In the following i refers to the ith base function coefficient.
+            The matrix has its entries j in its first Npts rows set to
+                A[i][j] = (Fourier transformed base function for point j) 
+                
+            We them choose a number of r-points, n_r, to evaluate the second
+            derivative of P(r) at. This is used as our regularization term.
+            For a vector r of length n_r, the following n_r rows are set to
+                A[i+Npts][j] = (2nd derivative of P(r), d**2(P(r))/d(r)**2, evaluated at r[j])
+                
+            The vector b has its first Npts entries set to
+                b[j] = (I(q) observed for point j)
+                
+            The following n_r entries are set to zero.
+            
+            The result is found by using scipy.linalg.basic.lstsq to invert
+            the matrix and find the coefficients x.
+            
+            @param nfunc: number of base functions to use.
+            @param nr: number of r points to evaluate the 2nd derivative at for the reg. term.
         """
         import math
@@ -239 +327 @@
         # b -- An M x nrhs matrix or M vector.
         npts = len(self.x)
-        nq   = 20
+        nq   = nr
         sqrt_alpha = math.sqrt(self.alpha)
         
@@ -250 +338 @@
                 if self._accept_q(self.x[i]):
                     a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i]
+                    
+            #TODO: refactor this: i_q should really be i_r
             for i_q in range(nq):
                 r = self.d_max/nq*i_q
@@ -298 +388 @@
         return c, err
         
-    def svd(self, nfunc=5):
-        import math, time
-        # Ac - b = 0
-        
-        A = numpy.zeros([nfunc, nfunc])
-        y = numpy.zeros(nfunc)
-        
-        t_0 = time.time()
-        for i in range(nfunc):
-            # A
-            for j in range(nfunc):
-                A[i][j] = 0.0
-                for k in range(len(self.x)):
-                    err = self.err[k]
-                    A[i][j] += 1.0/err/err*self.basefunc_ft(self.d_max, j+1, self.x[k]) \
-                            *self.basefunc_ft(self.d_max, i+1, self.x[k]);
-                #print A[i][j]
-                #A[i][j] -= self.alpha*(math.cos(math.pi*(i+j)) - math.cos(math.pi*(i-j)));
-                if i==j:
-                    A[i][j] += -1.0*self.alpha
-                elif i-j==1 or i-j==-1:
-                    A[i][j] += 1.0*self.alpha
-                #print "   ",A[i][j]
-            # y
-            y[i] = 0.0
-            for k in range(len(self.x)):
-                y[i] = self.y[k]/self.err[k]/self.err[k]*self.basefunc_ft(self.d_max, i+1, self.x[k])
-            
-        print time.time()-t_0, 'secs'
-        
-        # use numpy.pinv(A)
-        #inv_A = numpy.linalg.inv(A)
-        #c = y*inv_A
-        print y
-        c = numpy.linalg.solve(A, y)
-        
-        
-        print c
-        
-        err = numpy.zeros(len(c))
-        return c, err
-        
     def estimate_alpha(self, nfunc):
         """
@@ -392 +440 @@
                     
                     peaks = pr.get_peaks(out)
-                    print pr.alpha, peaks
                     if peaks>1:
                         found = True

pr_inversion/test/utest_invertor.py

@@ -176 +176 @@
         self.invertor.err = err
         # Perform inversion
-        out, cov = self.invertor.invert(10)
+        out, cov = self.invertor.invert_optimize(10)
         # This is a very specific case
         # We should make sure it always passes
@@ -404 +404 @@
     def test_qmin(self):
         self.invertor.q_min = 1.0
-        print self.invertor.q_min
         self.assertEqual(self.invertor.q_min, 1.0)