Changeset 2d06beb in sasview for pr_inversion/invertor.py

Timestamp:

May 5, 2008 2:03:39 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:

b43a009

Parents:

150c04a

Message:

Use simple least-square fit

File:

• pr_inversion/invertor.py (modified) (8 diffs)

pr_inversion/invertor.py

@@ -5 +5 @@
     
     ## Chisqr of the last computation
-    chisqr = 0
+    chi2  = 0
+    ## Time elapsed for last computation
+    elapsed = 0
     
     def __init__(self):
@@ -58 +60 @@
         return None
     
+    def clone(self):
+        """
+            Return a clone of this instance
+        """
+        invertor = Invertor()
+        invertor.chi2    = self.chi2 
+        invertor.elapsed = self.elapsed 
+        invertor.alpha   = self.alpha
+        invertor.d_max   = self.d_max
+        
+        invertor.x = self.x
+        invertor.y = self.y
+        invertor.err = self.err
+        
+        return invertor
+    
     def invert(self, nfunc=5):
         """
@@ -63 +81 @@
         """
         from scipy import optimize
+        import time
         
         # First, check that the current data is valid
@@ -69 +88 @@
         
         p = numpy.ones(nfunc)
+        t_0 = time.time()
         out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True)
         
@@ -78 +98 @@
         
         self.chi2 = chisqr
+
+        # Store computation time
+        self.elapsed = time.time() - t_0
         
         return out, cov_x
@@ -92 +115 @@
         
         p = numpy.ones(nfunc)
+        t_0 = time.time()
         out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True)
         
@@ -102 +126 @@
         self.chisqr = chisqr
         
+        # Store computation time
+        self.elapsed = time.time() - t_0
+
         return out, cov_x
     
@@ -117 +144 @@
 
         return self.get_pr_err(c, c_err, r)
+       
+    def lstsq(self, nfunc=5):
+        import math
+        from scipy.linalg.basic import lstsq
+        
+        # a -- An M x N matrix.
+        # b -- An M x nrhs matrix or M vector.
+        npts = len(self.x)
+        nq   = 20
+        sqrt_alpha = math.sqrt(self.alpha)
+        
+        a = numpy.zeros([npts+nq, nfunc])
+        b = numpy.zeros(npts+nq)
+        err = numpy.zeros(nfunc)
+        
+        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]
+            for i_q in range(nq):
+                r = self.d_max/nq*i_q
+                #a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*math.fabs(math.sin(math.pi*(j+1)*r/self.d_max) + math.pi*(j+1)*r/self.d_max * math.cos(math.pi*(j+1)*r/self.d_max))     
+                a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*(2.0*math.pi*(j+1)/self.d_max*math.cos(math.pi*(j+1)*r/self.d_max) + math.pi**2*(j+1)**2*r/self.d_max**2 * math.sin(math.pi*(j+1)*r/self.d_max))     
+        
+        for i in range(npts):
+            b[i] = self.y[i]/self.err[i]
+            
+        c, chi2, rank, n = lstsq(a, b)
+        self.chi2 = chi2
+                
+        at = numpy.transpose(a)
+        inv_cov = numpy.zeros([nfunc,nfunc])
+        for i in range(nfunc):
+            for j in range(nfunc):
+                inv_cov[i][j] = 0.0
+                for k in range(npts):
+                    inv_cov[i][j] = at[i][k]*a[k][j]
+                    
+        # Compute the reg term size for the output
+        sum_sig = 0.0
+        sum_reg = 0.0
+        for j in range(nfunc):
+            for i in range(npts):
+                sum_sig += (a[i][j])**2
+            for i in range(nq):
+                sum_reg += (a[i_q+npts][j])**2
+                    
+        new_alpha = sum_sig/(sum_reg/self.alpha)
+        print "Suggested alpha =", 0.1*new_alpha
+        
+        try:
+            err = math.fabs(chi2/(npts-nfunc))* inv_cov
+        except:
+            print "Error estimating uncertainties"
+            
+        
+        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
+        
+        
+