changeset 15:0d66a0aafb39

SL0_approx test working

--- a/matlab/SL0/SL0_approx.m	Sat Nov 05 22:06:11 2011 +0000
+++ b/matlab/SL0/SL0_approx.m	Sat Nov 05 22:10:06 2011 +0000
@@ -113,11 +113,7 @@
 % Initialization
 %s = A\x;
 s = A_pinv*x;
-% Nic:
-s = zeros(size(A,2),1);
-s = randn(size(A,2),1);
-sigma = 10;
-%sigma = 2*max(abs(s));
+sigma = 2*max(abs(s));
 
 % Main Loop
 while sigma>sigma_min

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyCSalgos/SL0/SL0_approx.py	Sat Nov 05 22:10:06 2011 +0000
@@ -0,0 +1,72 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Nov 05 21:29:09 2011
+
+@author: Nic
+"""
+
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Nov 05 18:39:54 2011
+
+@author: Nic
+"""
+
+import numpy as np
+
+#function s=SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s)
+def SL0_approx(A, x, eps, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None):
+  
+  if A_pinv is None:
+    A_pinv = np.linalg.pinv(A)
+  
+  if true_s is not None:
+      ShowProgress = True
+  else:
+      ShowProgress = False
+  
+  # Initialization
+  #s = A\x;
+  s = np.dot(A_pinv,x)
+  sigma = 2.0 * np.abs(s).max()
+  
+  # Main Loop
+  while sigma>sigma_min:
+      for i in np.arange(L):
+          delta = OurDelta(s,sigma)
+          s = s - mu_0*delta
+          # At this point, s no longer exactly satisfies x = A*s
+          # The original SL0 algorithm projects s onto {s | x = As} with
+          # s = s - np.dot(A_pinv,(np.dot(A,s)-x))   # Projection
+          # We want to project s onto {s | |x-As| < eps}
+          # We move onto the direction -A_pinv*(A*s-x), but only with a
+          # smaller step:
+          direction = np.dot(A_pinv,(np.dot(A,s)-x))
+          if (np.linalg.norm(np.dot(A,direction)) >= eps):
+            s = s - (1.0 - eps/np.linalg.norm(np.dot(A,direction))) * direction
+
+          #assert(np.linalg.norm(x - np.dot(A,s)) < eps + 1e-6)          
+      
+      if ShowProgress:
+          #fprintf('     sigma=#f, SNR=#f\n',sigma,estimate_SNR(s,true_s))
+          string = '     sigma=%f, SNR=%f\n' % sigma,estimate_SNR(s,true_s)
+          print string
+      
+      sigma = sigma * sigma_decrease_factor
+  
+  return s
+  
+  
+####################################################################
+#function delta=OurDelta(s,sigma)
+def OurDelta(s,sigma):
+  
+  return s * np.exp( (-np.abs(s)**2) / sigma**2)
+  
+####################################################################
+#function SNR=estimate_SNR(estim_s,true_s)
+def estimate_SNR(estim_s, true_s):
+  
+  err = true_s - estim_s
+  return 10*np.log10((true_s**2).sum()/(err**2).sum())
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/sl0approx_test.py	Sat Nov 05 22:10:06 2011 +0000
@@ -0,0 +1,69 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Nov 05 21:34:49 2011
+
+@author: Nic
+"""
+import numpy as np
+import numpy.linalg
+import scipy.io
+import unittest
+from pyCSalgos.SL0.SL0_approx import SL0_approx
+
+class SL0results(unittest.TestCase):
+  def testResults(self):
+    mdict = scipy.io.loadmat('SL0approxtestdata.mat')
+    
+    # A = system matrix
+    # Y = matrix with measurements (on columns)
+    # sigmamin = vector with sigma_min
+    for A,Y,eps,sigmamin,Xr in zip(mdict['cellA'].squeeze(),mdict['cellY'].squeeze(),mdict['cellEps'].squeeze(),mdict['sigmamin'].squeeze(),mdict['cellXr'].squeeze()):
+      for i in np.arange(Y.shape[1]):
+        
+        # Fix numpy error "LapackError: Parameter a has non-native byte order in lapack_lite.dgesdd"
+        A = A.newbyteorder('=')
+        Y = Y.newbyteorder('=')
+        eps = eps.newbyteorder('=')
+        sigmamin = sigmamin.newbyteorder('=')
+        Xr = Xr.newbyteorder('=')
+        
+        xr = SL0_approx(A, Y[:,i], eps.squeeze()[i], sigmamin)
+        
+        # check if found solution is the same as the correct cslution
+        diff = numpy.linalg.norm(xr - Xr[:,i])
+        self.assertTrue(diff < 1e-12)
+    #        err1 = numpy.linalg.norm(Y[:,i] - np.dot(A,xr))
+    #        err2 = numpy.linalg.norm(Y[:,i] - np.dot(A,Xr[:,i]))
+    #        norm1 = numpy.linalg.norm(xr,1)
+    #        norm2 = numpy.linalg.norm(Xr[:,i],1)
+    #                
+    #        # Make a more robust condition:
+    #        #  OK;    if   solutions are close enough (diff < 1e-6)
+    #        #              or
+    #        #              (
+    #        #               Python solution fulfills the constraint better (or up to 1e-6 worse)
+    #        #                 and
+    #        #               Python solution has l1 norm no more than 1e-6 larger as the reference solution
+    #        #                 (i.e. either norm1 < norm2   or   norm1>norm2 not by more than 1e-6)
+    #        #              )
+    #        #        
+    #        #  ERROR: else        
+    #        differr  = err1 - err2    # intentionately no abs(), since err1` < err2 is good
+    #        diffnorm = norm1 - norm2  # intentionately no abs(), since norm1 < norm2 is good
+    #        if diff < 1e-6 or (differr < 1e-6 and (diffnorm < 1e-6)):
+    #          isok = True
+    #        else:
+    #          isok = False
+    #        self.assertTrue(isok)
+        
+        #diff = numpy.linalg.norm(xr - Xr[:,i])
+        #if diff > 1e-6:
+        #    self.assertTrue(diff < 1e-6)
+
+  
+if __name__ == "__main__":
+    #import cProfile
+    #cProfile.run('unittest.main()', 'profres')
+    unittest.main()    
+    #suite = unittest.TestLoader().loadTestsFromTestCase(CompareResults)
+    #unittest.TextTestRunner(verbosity=2).run(suite)