changeset 22:2dd78e37b23a

ABS approx script is working
Started working on parallel

--- a/scripts/ABSapprox.py	Tue Nov 08 14:45:35 2011 +0000
+++ b/scripts/ABSapprox.py	Wed Nov 09 00:11:14 2011 +0000
@@ -6,61 +6,133 @@
 """
 
 import numpy as np
+import scipy.io
+import math
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
 import pyCSalgos
 import pyCSalgos.GAP.GAP
 import pyCSalgos.SL0.SL0_approx
 
 # Define functions that prepare arguments for each algorithm call
-def gap_paramsetup(y,M,Omega,epsilon,lbd):
+def run_gap(y,M,Omega,epsilon):
   gapparams = {"num_iteration" : 1000,\
                    "greedy_level" : 0.9,\
                    "stopping_coefficient_size" : 1e-4,\
                    "l2solver" : 'pseudoinverse',\
                    "noise_level": epsilon}
-  return y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1])
-def sl0_paramsetup(y,M,Omega,epsilon,lbd):
+  return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1]))[0]
+ 
+def run_sl0(y,M,Omega,D,U,S,Vt,epsilon,lbd):
   
   N,n = Omega.shape
-  D = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
+  #D = np.linalg.pinv(Omega)
+  #U,S,Vt = np.linalg.svd(D)
   aggDupper = np.dot(M,D)
   aggDlower = Vt[-(N-n):,:]
   aggD = np.concatenate((aggDupper, lbd * aggDlower))
   aggy = np.concatenate((y, np.zeros(N-n)))
   
-  sigmamin = 0.01
-  sigma_decrease_factor = 0.8
+  sigmamin = 0.001
+  sigma_decrease_factor = 0.5
   mu_0 = 2
   L = 10
-  return aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L
-
-def post_multiply_with_D(D,gamma):
-    return np.dot(D,gamma)
-def post_do_nothing(D,gamma):
-    return gamma
+  return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L)
 
 # Define tuples (algorithm setup function, algorithm function, name)
-gap = (gap_paramsetup, pyCSalgos.GAP.GAP.GAP, post_do_nothing, 'GAP')
-sl0 = (sl0_paramsetup, pyCSalgos.SL0.SL0_approx.SL0_approx, post_multiply_with_D, 'SL0_approx')
-#sl0 = (sl0_paramsetup, lambda x: np.dot(x[0],x[1]()), 'SL0_approx')
+gap = (run_gap, 'GAP')
+sl0 = (run_sl0, 'SL0_approx')
 
-# Main function
+# Define which algorithms to run
+#  1. Algorithms not depending on lambda
+algosN = gap,   # tuple
+#  2. Algorithms depending on lambda (our ABS approach)
+algosL = sl0,   # tuple
+
 def mainrun():
-
-  # Define which algorithms to run
-  algos = (gap, sl0)
-  numalgos = len(algos)
   
-  # Set up experiment parameters
-  sigma = 2.0;
-  delta = 0.8;
-  rho   = 0.15;
+  nalgosN = len(algosN)  
+  nalgosL = len(algosL)
+  
+  #Set up experiment parameters
+  d = 50;
+  sigma = 2.0
+  #deltas = np.arange(0.05,0.95,0.05)
+  #rhos = np.arange(0.05,0.95,0.05)
+  deltas = np.array([0.05,0.95])
+  rhos = np.array([0.05,0.95])
+  #deltas = np.array([0.05])
+  #rhos = np.array([0.05])
+  #delta = 0.8;
+  #rho   = 0.15;
   numvects = 10; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.concatenate((np.array([0]), 10**np.linspace(-5, 4, 10)))
+
+  meanmatrix = dict()
+  for i,algo in zip(np.arange(nalgosN),algosN):
+    meanmatrix[algo[1]]   = np.zeros((rhos.size, deltas.size))
+  for i,algo in zip(np.arange(nalgosL),algosL):
+    meanmatrix[algo[1]]   = np.zeros((lambdas.size, rhos.size, deltas.size))
+  
+  for idelta,delta in zip(np.arange(deltas.size),deltas):
+    for irho,rho in zip(np.arange(rhos.size),rhos):
+      
+      # Generate data and operator
+      Omega,x0,y,M,realnoise = genData(d,sigma,delta,rho,numvects,SNRdb)
+      
+      # Run algorithms
+      mrelerrN,mrelerrL = runonce(algosN,algosL,Omega,y,lambdas,realnoise,M,x0)
+      
+      for algotuple in algosN: 
+        meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[algotuple[1]]
+        if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]):
+          meanmatrix[algotuple[1]][irho,idelta] = 0
+      for algotuple in algosL:
+        for ilbd in np.arange(lambdas.size):
+          meanmatrix[algotuple[1]][ilbd,irho,idelta] = 1 - mrelerrL[algotuple[1]][ilbd]
+          if meanmatrix[algotuple[1]][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][ilbd,irho,idelta]):
+            meanmatrix[algotuple[1]][ilbd,irho,idelta] = 0
+   
+  #  # Prepare matrices to show
+  #  showmats = dict()
+  #  for i,algo in zip(np.arange(nalgosN),algosN):
+  #    showmats[algo[1]]   = np.zeros(rhos.size, deltas.size)
+  #  for i,algo in zip(np.arange(nalgosL),algosL):
+  #    showmats[algo[1]]   = np.zeros(lambdas.size, rhos.size, deltas.size)
+
+    # Save
+    tosave = dict()
+    tosave['meanmatrix'] = meanmatrix
+    tosave['d'] = d
+    tosave['sigma'] = sigma
+    tosave['deltas'] = deltas
+    tosave['rhos'] = rhos
+    tosave['numvects'] = numvects
+    tosave['SNRdb'] = SNRdb
+    tosave['lambdas'] = lambdas
+    try:
+      scipy.io.savemat('ABSapprox.mat',tosave)
+    except TypeError:
+      print "Oops, Type Error"
+      raise    
+  # Show
+  for algotuple in algosN:
+    plt.figure()
+    plt.imshow(meanmatrix[algotuple[1]], cmap=cm.gray, interpolation='nearest')
+  for algotuple in algosL:
+    for ilbd in np.arange(lambdas.size):
+      plt.figure()
+      plt.imshow(meanmatrix[algotuple[1]][ilbd], cmap=cm.gray, interpolation='nearest')
+  plt.show()
+  print "Finished."
+  
+def genData(d,sigma,delta,rho,numvects,SNRdb):
 
   # Process parameters
   noiselevel = 1.0 / (10.0**(SNRdb/10.0));
-  d = 50;
   p = round(sigma*d);
   m = round(delta*d);
   l = round(d - rho*m);
@@ -68,43 +140,73 @@
   # Generate Omega and data based on parameters
   Omega = pyCSalgos.GAP.GAP.Generate_Analysis_Operator(d, p);
   # Optionally make Omega more coherent
-  #[U, S, Vt] = np.linalg.svd(Omega);
-  #Sdnew = np.diag(S) * (1+np.arange(np.diag(S).size)); % Make D coherent, not Omega!
-  #Snew = [diag(Sdnew); zeros(size(S,1) - size(S,2), size(S,2))];
-  #Omega = U * Snew * V';
+  U,S,Vt = np.linalg.svd(Omega);
+  Sdnew = S * (1+np.arange(S.size)) # Make D coherent, not Omega!
+  Snew = np.vstack((np.diag(Sdnew), np.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
+  Omega = np.dot(U , np.dot(Snew,Vt))
 
   # Generate data  
   x0,y,M,Lambda,realnoise = pyCSalgos.GAP.GAP.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
+  
+  return Omega,x0,y,M,realnoise
 
-  # Values for lambda
-  #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.concatenate((np.array([0]), 10**np.linspace(-5, 4, 10)))
+def runonce(algosN,algosL,Omega,y,lambdas,realnoise,M,x0):
+  
+  d = Omega.shape[1]  
+  
+  nalgosN = len(algosN)  
+  nalgosL = len(algosL)
   
   xrec = dict()
   err = dict()
   relerr = dict()
-  for i,algo in zip(np.arange(numalgos),algos):
-    xrec[algo[3]]   = np.zeros((lambdas.size, d, y.shape[1]))
-    err[algo[3]]    = np.zeros((lambdas.size, y.shape[1]))
-    relerr[algo[3]] = np.zeros((lambdas.size, y.shape[1]))
+
+  # Prepare storage variables for algorithms non-Lambda
+  for i,algo in zip(np.arange(nalgosN),algosN):
+    xrec[algo[1]]   = np.zeros((d, y.shape[1]))
+    err[algo[1]]    = np.zeros(y.shape[1])
+    relerr[algo[1]] = np.zeros(y.shape[1])
+  # Prepare storage variables for algorithms with Lambda    
+  for i,algo in zip(np.arange(nalgosL),algosL):
+    xrec[algo[1]]   = np.zeros((lambdas.size, d, y.shape[1]))
+    err[algo[1]]    = np.zeros((lambdas.size, y.shape[1]))
+    relerr[algo[1]] = np.zeros((lambdas.size, y.shape[1]))
   
+  # Run algorithms non-Lambda
+  for iy in np.arange(y.shape[1]):
+    for algofunc,strname in algosN:
+      epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+      xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega,epsilon)
+      err[strname][iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
+      relerr[strname][iy] = err[strname][iy] / np.linalg.norm(x0[:,iy])
+  for algotuple in algosN:
+    print algotuple[1],' : avg relative error = ',np.mean(relerr[strname])  
+
+  # Run algorithms with Lambda
   for ilbd,lbd in zip(np.arange(lambdas.size),lambdas):
     for iy in np.arange(y.shape[1]):
-      for algosetupfunc,algofunc,algopostfunc,strname in algos:
+      D = np.linalg.pinv(Omega)
+      U,S,Vt = np.linalg.svd(D)
+      for algofunc,strname in algosL:
         epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
-        
-        inparams = algosetupfunc(y[:,iy],M,Omega,epsilon,lbd)
-        xrec[strname][ilbd,:,iy] = algopostfunc(algofunc(*inparams)[0])
-        
+        gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
+        xrec[strname][ilbd,:,iy] = np.dot(D,gamma)
         err[strname][ilbd,iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy])
         relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / np.linalg.norm(x0[:,iy])
-        
     print 'Lambda = ',lbd,' :'
-    for strname in relerr:
-      print '   ',strname,' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
+    for algotuple in algosL:
+      print '   ',algotuple[1],' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
 
-
-
+  # Prepare results
+  mrelerrN = dict()
+  for algotuple in algosN:
+    mrelerrN[algotuple[1]] = np.mean(relerr[algotuple[1]])
+  mrelerrL = dict()
+  for algotuple in algosL:
+    mrelerrL[algotuple[1]] = np.mean(relerr[algotuple[1]],1)
+  
+  return mrelerrN,mrelerrL
+  
 # Script main
 if __name__ == "__main__":
   mainrun()
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/ABSapproxPP.py	Wed Nov 09 00:11:14 2011 +0000
@@ -0,0 +1,229 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Nov 05 18:08:40 2011
+
+@author: Nic
+"""
+
+import numpy as np
+import scipy.io
+import math
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import pp
+import pyCSalgos
+import pyCSalgos.GAP.GAP
+import pyCSalgos.SL0.SL0_approx
+
+# Define functions that prepare arguments for each algorithm call
+def run_gap(y,M,Omega,epsilon):
+  gapparams = {"num_iteration" : 1000,\
+                   "greedy_level" : 0.9,\
+                   "stopping_coefficient_size" : 1e-4,\
+                   "l2solver" : 'pseudoinverse',\
+                   "noise_level": epsilon}
+  return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1]))[0]
+ 
+def run_sl0(y,M,Omega,D,U,S,Vt,epsilon,lbd):
+  
+  N,n = Omega.shape
+  #D = np.linalg.pinv(Omega)
+  #U,S,Vt = np.linalg.svd(D)
+  aggDupper = np.dot(M,D)
+  aggDlower = Vt[-(N-n):,:]
+  aggD = np.concatenate((aggDupper, lbd * aggDlower))
+  aggy = np.concatenate((y, np.zeros(N-n)))
+  
+  sigmamin = 0.001
+  sigma_decrease_factor = 0.5
+  mu_0 = 2
+  L = 10
+  return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L)
+
+# Define tuples (algorithm setup function, algorithm function, name)
+gap = (run_gap, 'GAP')
+sl0 = (run_sl0, 'SL0_approx')
+
+# Define which algorithms to run
+#  1. Algorithms not depending on lambda
+algosN = gap,   # tuple
+#  2. Algorithms depending on lambda (our ABS approach)
+algosL = sl0,   # tuple
+
+def mainrun():
+  
+  nalgosN = len(algosN)  
+  nalgosL = len(algosL)
+  
+  #Set up experiment parameters
+  d = 50;
+  sigma = 2.0
+  #deltas = np.arange(0.05,0.95,0.05)
+  #rhos = np.arange(0.05,0.95,0.05)
+  deltas = np.array([0.05,0.95])
+  rhos = np.array([0.05,0.95])
+  #deltas = np.array([0.05])
+  #rhos = np.array([0.05])
+  #delta = 0.8;
+  #rho   = 0.15;
+  numvects = 10; # Number of vectors to generate
+  SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.concatenate((np.array([0]), 10**np.linspace(-5, 4, 10)))
+
+  meanmatrix = dict()
+  for i,algo in zip(np.arange(nalgosN),algosN):
+    meanmatrix[algo[1]]   = np.zeros((rhos.size, deltas.size))
+  for i,algo in zip(np.arange(nalgosL),algosL):
+    meanmatrix[algo[1]]   = np.zeros((lambdas.size, rhos.size, deltas.size))
+
+  # PP: start job server  
+  job_server = pp.Server(ncpus = 1) 
+  idx = 0
+  jobparams = []
+  for idelta,delta in zip(np.arange(deltas.size),deltas):
+    for irho,rho in zip(np.arange(rhos.size),rhos):
+      
+      # Generate data and operator
+      Omega,x0,y,M,realnoise = genData(d,sigma,delta,rho,numvects,SNRdb)
+      
+      jobparams.append((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
+      
+      idx = idx + 1
+      
+  # Run algorithms
+  jobs = [job_server.submit(runonce, jobparam, (run_gap,run_sl0), ('numpy',)) for jobparam in jobparams]
+  #funcarray[idelta,irho] = job_server.submit(runonce,(algosN,algosL, Omega,y,lambdas,realnoise,M,x0), (run_gap,run_sl0))
+      #mrelerrN,mrelerrL = runonce(algosN,algosL,Omega,y,lambdas,realnoise,M,x0)
+
+  # Get data from jobs
+  idx = 0
+  for idelta,delta in zip(np.arange(deltas.size),deltas):
+    for irho,rho in zip(np.arange(rhos.size),rhos):
+      mrelerrN,mrelerrL = jobs[idx]()
+      for algotuple in algosN: 
+        meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[algotuple[1]]
+        if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]):
+          meanmatrix[algotuple[1]][irho,idelta] = 0
+      for algotuple in algosL:
+        for ilbd in np.arange(lambdas.size):
+          meanmatrix[algotuple[1]][ilbd,irho,idelta] = 1 - mrelerrL[algotuple[1]][ilbd]
+          if meanmatrix[algotuple[1]][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][ilbd,irho,idelta]):
+            meanmatrix[algotuple[1]][ilbd,irho,idelta] = 0
+      idx = idx + 1
+   
+  #  # Prepare matrices to show
+  #  showmats = dict()
+  #  for i,algo in zip(np.arange(nalgosN),algosN):
+  #    showmats[algo[1]]   = np.zeros(rhos.size, deltas.size)
+  #  for i,algo in zip(np.arange(nalgosL),algosL):
+  #    showmats[algo[1]]   = np.zeros(lambdas.size, rhos.size, deltas.size)
+
+    # Save
+    tosave = dict()
+    tosave['meanmatrix'] = meanmatrix
+    tosave['d'] = d
+    tosave['sigma'] = sigma
+    tosave['deltas'] = deltas
+    tosave['rhos'] = rhos
+    tosave['numvects'] = numvects
+    tosave['SNRdb'] = SNRdb
+    tosave['lambdas'] = lambdas
+    try:
+      scipy.io.savemat('ABSapprox.mat',tosave)
+    except TypeError:
+      print "Oops, Type Error"
+      raise    
+  # Show
+  for algotuple in algosN:
+    plt.figure()
+    plt.imshow(meanmatrix[algotuple[1]], cmap=cm.gray, interpolation='nearest')
+  for algotuple in algosL:
+    for ilbd in np.arange(lambdas.size):
+      plt.figure()
+      plt.imshow(meanmatrix[algotuple[1]][ilbd], cmap=cm.gray, interpolation='nearest')
+  plt.show()
+  print "Finished."
+  
+def genData(d,sigma,delta,rho,numvects,SNRdb):
+
+  # Process parameters
+  noiselevel = 1.0 / (10.0**(SNRdb/10.0));
+  p = round(sigma*d);
+  m = round(delta*d);
+  l = round(d - rho*m);
+  
+  # Generate Omega and data based on parameters
+  Omega = pyCSalgos.GAP.GAP.Generate_Analysis_Operator(d, p);
+  # Optionally make Omega more coherent
+  U,S,Vt = np.linalg.svd(Omega);
+  Sdnew = S * (1+np.arange(S.size)) # Make D coherent, not Omega!
+  Snew = np.vstack((np.diag(Sdnew), np.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
+  Omega = np.dot(U , np.dot(Snew,Vt))
+
+  # Generate data  
+  x0,y,M,Lambda,realnoise = pyCSalgos.GAP.GAP.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
+  
+  return Omega,x0,y,M,realnoise
+
+def runonce(algosN,algosL,Omega,y,lambdas,realnoise,M,x0):
+  
+  d = Omega.shape[1]  
+  
+  nalgosN = len(algosN)  
+  nalgosL = len(algosL)
+  
+  xrec = dict()
+  err = dict()
+  relerr = dict()
+
+  # Prepare storage variables for algorithms non-Lambda
+  for i,algo in zip(np.arange(nalgosN),algosN):
+    xrec[algo[1]]   = np.zeros((d, y.shape[1]))
+    err[algo[1]]    = np.zeros(y.shape[1])
+    relerr[algo[1]] = np.zeros(y.shape[1])
+  # Prepare storage variables for algorithms with Lambda    
+  for i,algo in zip(np.arange(nalgosL),algosL):
+    xrec[algo[1]]   = np.zeros((lambdas.size, d, y.shape[1]))
+    err[algo[1]]    = np.zeros((lambdas.size, y.shape[1]))
+    relerr[algo[1]] = np.zeros((lambdas.size, y.shape[1]))
+  
+  # Run algorithms non-Lambda
+  for iy in np.arange(y.shape[1]):
+    for algofunc,strname in algosN:
+      epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+      xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega,epsilon)
+      err[strname][iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
+      relerr[strname][iy] = err[strname][iy] / np.linalg.norm(x0[:,iy])
+  for algotuple in algosN:
+    print algotuple[1],' : avg relative error = ',np.mean(relerr[strname])  
+
+  # Run algorithms with Lambda
+  for ilbd,lbd in zip(np.arange(lambdas.size),lambdas):
+    for iy in np.arange(y.shape[1]):
+      D = np.linalg.pinv(Omega)
+      U,S,Vt = np.linalg.svd(D)
+      for algofunc,strname in algosL:
+        epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+        gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
+        xrec[strname][ilbd,:,iy] = np.dot(D,gamma)
+        err[strname][ilbd,iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy])
+        relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / np.linalg.norm(x0[:,iy])
+    print 'Lambda = ',lbd,' :'
+    for algotuple in algosL:
+      print '   ',algotuple[1],' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
+
+  # Prepare results
+  mrelerrN = dict()
+  for algotuple in algosN:
+    mrelerrN[algotuple[1]] = np.mean(relerr[algotuple[1]])
+  mrelerrL = dict()
+  for algotuple in algosL:
+    mrelerrL[algotuple[1]] = np.mean(relerr[algotuple[1]],1)
+  
+  return mrelerrN,mrelerrL
+  
+# Script main
+if __name__ == "__main__":
+  mainrun()
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/study_analysis_rec_algos_noisy.m	Wed Nov 09 00:11:14 2011 +0000
@@ -0,0 +1,417 @@
+% File: study_analysis_rec_algos
+% Run global experiment to compare algorithms used for analysis-based reconstruction
+% and plot phast transition graphs
+
+clear all
+close all
+
+% =================================
+% Set up experiment parameters
+%==================================
+% Which form factor, delta and rho we want
+sigmas = 1.2;
+deltas = 0.05:0.05:0.95;
+rhos   = 0.05:0.05:0.95;
+% deltas = [0.95];
+% rhos   = [0.1];
+%deltas = 0.3:0.3:0.9;
+%rhos   = 0.3:0.3:0.9;
+
+% Number of vectors to generate each time
+numvects = 100;
+
+% Add noise 
+% This is norm(signal)/norm(noise), so power, not energy
+SNRdb = 20; % virtually no noise
+
+% Show progressbar ? (not recommended when running on parallel threads)
+do_progressbar = 0;
+
+% Value of lambda
+lambda = 1e-2;
+
+% What algos to run
+do_abs_ompk = 1;
+do_abs_ompeps = 1;
+do_abs_tst = 1;
+do_abs_bp = 1;
+do_gap = 1;
+do_nesta = 0;
+
+% =================================
+% Processing the parameters
+%==================================
+% Set up parameter structure
+count = 0;
+for isigma = 1:sigmas
+    for idelta = 1:numel(deltas)
+        for irho = 1:numel(rhos)
+            sigma = sigmas(isigma);
+            delta = deltas(idelta);
+            rho = rhos(irho);
+        
+            d = 200;
+            p = round(sigma*d);
+            m = round(delta*d);
+            l = round(d - rho*m);
+            
+            params(count+1).d = d;
+            params(count+1).p = p;
+            params(count+1).m = m;
+            params(count+1).l = l;
+            
+            count = count + 1;
+        end
+    end
+end
+
+% Compute noiselevel from db
+noiselevel = 1 / (10^(SNRdb/10));
+
+%load study_analysis_init
+
+% Generate an analysis operator Omega
+Omega = Generate_Analysis_Operator(d, p);
+
+% Progressbar
+if do_progressbar
+    progressbar('Total', 'Current parameters');
+end
+
+% Init times
+elapsed_abs_ompk = 0;
+elapsed_abs_ompeps = 0;
+elapsed_abs_tst = 0;
+elapsed_abs_bp = 0;
+elapsed_gap = 0;
+elapsed_nesta = 0;
+
+% Init results structure
+results = [];
+
+% Prepare progressbar reduction variables
+% if do_progressbar
+%     incr2 = 1/numvects;
+%     incr1 = 1/numvects/count;
+%     frac2 = 0;
+%     frac1 = 0;
+% end 
+
+% ========
+% Run
+% ========
+parfor iparam = 1:numel(params)
+    
+    % Read current parameters
+    d = params(iparam).d;
+    p = params(iparam).p;
+    m = params(iparam).m;
+    l = params(iparam).l;
+    
+    % Init stuff
+    xrec_abs_ompk   = zeros(d, numvects);
+    xrec_abs_ompeps = zeros(d, numvects);
+    xrec_abs_tst    = zeros(d, numvects);
+    xrec_abs_bp     = zeros(d, numvects);
+    xrec_gap        = zeros(d, numvects);
+    xrec_nesta      = zeros(d, numvects);
+    %
+       err_abs_ompk   = zeros(numvects,1);
+    relerr_abs_ompk   = zeros(numvects,1);
+       err_abs_ompeps = zeros(numvects,1);
+    relerr_abs_ompeps = zeros(numvects,1);
+       err_abs_tst    = zeros(numvects,1);
+    relerr_abs_tst    = zeros(numvects,1);
+       err_abs_bp     = zeros(numvects,1);
+    relerr_abs_bp     = zeros(numvects,1);
+       err_gap        = zeros(numvects,1);
+    relerr_gap        = zeros(numvects,1);
+       err_nesta      = zeros(numvects,1);
+    relerr_nesta      = zeros(numvects,1);
+
+    % Generate data based on parameters
+    [x0,y,M,Lambda] = Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
+    
+    % For every generated signal do
+    for iy = 1:size(y,2)
+        
+        % Compute epsilon 
+        epsilon = noiselevel * norm(y(:,iy));
+        
+        %--------------------------------
+        % Reconstruct (and measure delay)
+        % Compute reconstruction error
+        %--------------------------------
+        % ABS-OMPk
+        if (do_abs_ompk)
+            timer_abs_ompk = tic;
+            xrec_abs_ompk(:,iy) = ABS_OMPk_approx(y(:,iy), Omega, M, p-l, lambda);
+            elapsed_abs_ompk = elapsed_abs_ompk + toc(timer_abs_ompk);
+            %
+            err_abs_ompk(iy)    = norm(x0(:,iy) - xrec_abs_ompk(:,iy));
+            relerr_abs_ompk(iy) = err_abs_ompk(iy) / norm(x0(:,iy));   
+        end
+        % ABS-OMPeps
+        if (do_abs_ompeps)
+            timer_abs_ompeps = tic;
+            xrec_abs_ompeps(:,iy) = ABS_OMPeps_approx(y(:,iy), Omega, M, epsilon, lambda);
+            elapsed_abs_ompeps = elapsed_abs_ompeps + toc(timer_abs_ompeps);
+            %
+            err_abs_ompeps(iy)    = norm(x0(:,iy) - xrec_abs_ompeps(:,iy));
+            relerr_abs_ompeps(iy) = err_abs_ompeps(iy) / norm(x0(:,iy));
+        end
+        % ABS-TST
+        if (do_abs_tst)
+            timer_abs_tst = tic;
+            xrec_abs_tst(:,iy) = ABS_TST_approx(y(:,iy), Omega, M, epsilon, lambda);
+            elapsed_abs_tst = elapsed_abs_tst + toc(timer_abs_tst);
+            %
+            err_abs_tst(iy)     = norm(x0(:,iy) - xrec_abs_tst(:,iy));
+            relerr_abs_tst(iy)  = err_abs_tst(iy) / norm(x0(:,iy));
+        end
+        % ABS-BP
+        if (do_abs_bp)
+            timer_abs_bp = tic;
+            xrec_abs_bp(:,iy)  = ABS_BP_approx(y(:,iy), Omega, M, epsilon, lambda);
+            elapsed_abs_bp = elapsed_abs_bp + toc(timer_abs_bp);
+            %
+            err_abs_bp(iy)     = norm(x0(:,iy) - xrec_abs_bp(:,iy));
+            relerr_abs_bp(iy)  = err_abs_bp(iy) / norm(x0(:,iy));
+        end
+        % GAP
+        if (do_gap)
+            gapparams = [];
+            gapparams.num_iteration = 40;
+            gapparams.greedy_level = 0.9;
+            gapparams.stopping_coefficient_size = 1e-4;
+            gapparams.l2solver = 'pseudoinverse';
+            gapparams.noise_level = noiselevel;
+            timer_gap = tic;
+            xrec_gap(:,iy) = GAP(y(:,iy), M, M', Omega, Omega', gapparams, zeros(d,1));
+            elapsed_gap = elapsed_gap + toc(timer_gap);
+            %
+            err_gap(iy)     = norm(x0(:,iy) - xrec_gap(:,iy));
+            relerr_gap(iy)  = err_gap(iy) / norm(x0(:,iy));
+        end
+        % NESTA
+        if (do_nesta)
+            try
+                timer_nesta = tic;
+                xrec_nesta(:,iy) = do_nesta_DemoNonProjector(x0(:,iy), M, Omega', 0);
+                elapsed_nesta = elapsed_nesta + toc(timer_nesta);
+            catch err
+                disp('*****ERROR: NESTA throwed error *****');
+                xrec_nesta(:,iy) = zeros(size(x0(:,iy)));
+            end
+            %
+            err_nesta(iy)       = norm(x0(:,iy) - xrec_nesta(:,iy));
+            relerr_nesta(iy)    = err_nesta(iy) / norm(x0(:,iy)); 
+        end
+
+        % Update progressbar
+%         if do_progressbar
+%             %frac2 = iy/numvects;
+%             %frac1 = ((iparam-1) + frac2) / count;
+%             if norm(frac2 - 1) < 1e-6
+%                 frac2 = 0;
+%             end
+%             frac2 = frac2 + incr2;
+%             frac1 = frac1 + incr1;
+%             progressbar(frac1, frac2);
+%         end
+    end
+    
+    %--------------------------------
+    % Save results in big stucture & display
+    %--------------------------------
+    % Save reconstructed signals
+    % Save rel & abs errors
+    % Display error
+    disp(['Simulation no. ' num2str(iparam)]);
+    if (do_abs_ompk)
+        results(iparam).xrec_abs_ompk   = xrec_abs_ompk;
+        results(iparam).err_abs_ompk    = err_abs_ompk;
+        results(iparam).relerr_abs_ompk = relerr_abs_ompk;
+        disp(['  ABS_OMPk:   avg relative error = ' num2str(mean(relerr_abs_ompk))]);
+    end
+    if (do_abs_ompeps)
+        results(iparam).xrec_abs_ompeps   = xrec_abs_ompeps;
+        results(iparam).err_abs_ompeps    = err_abs_ompeps;
+        results(iparam).relerr_abs_ompeps = relerr_abs_ompeps;   
+        disp(['  ABS_OMPeps: avg relative error = ' num2str(mean(relerr_abs_ompeps))]);
+    end
+    if (do_abs_tst)
+        results(iparam).xrec_abs_tst   = xrec_abs_tst;
+        results(iparam).err_abs_tst    = err_abs_tst;
+        results(iparam).relerr_abs_tst = relerr_abs_tst;
+        disp(['  ABS_TST:    avg relative error = ' num2str(mean(relerr_abs_tst))]);
+    end
+    if (do_abs_bp)
+        results(iparam).xrec_abs_bp   = xrec_abs_bp;
+        results(iparam).err_abs_bp    = err_abs_bp;
+        results(iparam).relerr_abs_bp = relerr_abs_bp;
+        disp(['  ABS_BP:     avg relative error = ' num2str(mean(relerr_abs_bp))]);
+    end
+    if (do_gap)
+        results(iparam).xrec_gap   = xrec_gap;
+        results(iparam).err_gap    = err_gap;
+        results(iparam).relerr_gap = relerr_gap;
+        disp(['  GAP:        avg relative error = ' num2str(mean(relerr_gap))]);
+    end
+    if (do_nesta)
+        results(iparam).xrec_nesta   = xrec_nesta;
+        results(iparam).err_nesta    = err_nesta;
+        results(iparam).relerr_nesta = relerr_nesta;
+        disp(['  NESTA:      avg relative error = ' num2str(mean(relerr_nesta))]);
+    end
+end
+
+% =================================
+% Save
+% =================================
+save mat/avgerr_SNR20_lbd1e-2
+
+% =================================
+% Plot phase transition
+% =================================
+%--------------------------------
+% Prepare
+%--------------------------------
+%d = 200;
+%m = 190;
+%exactthr = d/m * noiselevel;
+%sigma = 1.2;
+iparam = 1;
+for idelta = 1:numel(deltas)
+    for irho = 1:numel(rhos)
+        % Create exact recovery count matrix 
+%         nexact_abs_bp  (irho, idelta)    = sum(results(iparam).relerr_abs_bp < exactthr);
+%         nexact_abs_ompk (irho, idelta)   = sum(results(iparam).relerr_abs_ompk < exactthr);
+%         nexact_abs_ompeps (irho, idelta) = sum(results(iparam).relerr_abs_ompeps < exactthr);
+%         nexact_gap (irho, idelta)        = sum(results(iparam).relerr_gap < exactthr);
+%         nexact_abs_tst (irho, idelta)    = sum(results(iparam).relerr_abs_tst < exactthr);
+% %         nexact_nesta(irho, idelta)       = sum(results(iparam).relerr_nesta < exactthr);
+
+        % Get histogram (for a single param set only!)
+%         hist_abs_ompk   = hist(results(iparam).relerr_abs_ompk, 0.001:0.001:0.1);
+%         hist_abs_ompeps = hist(results(iparam).relerr_abs_ompeps, 0.001:0.001:0.1);
+%         hist_abs_tst    = hist(results(iparam).relerr_abs_tst, 0.001:0.001:0.1);
+%         hist_abs_bp     = hist(results(iparam).relerr_abs_bp, 0.001:0.001:0.1);
+%         hist_gap        = hist(results(iparam).relerr_gap, 0.001:0.001:0.1);
+        
+        % Compute average error value
+        if (do_abs_ompk)
+            avgerr_abs_ompk(irho, idelta)    = 1 - mean(results(iparam).relerr_abs_ompk);
+            avgerr_abs_ompk(avgerr_abs_ompk < 0) = 0;
+        end
+        if (do_abs_ompeps)
+            avgerr_abs_ompeps(irho, idelta)  = 1 - mean(results(iparam).relerr_abs_ompeps);
+            avgerr_abs_ompeps(avgerr_abs_ompeps < 0) = 0;
+        end
+        if (do_abs_tst)
+            avgerr_abs_tst(irho, idelta)     = 1 - mean(results(iparam).relerr_abs_tst);
+            avgerr_abs_tst(avgerr_abs_tst < 0) = 0;
+        end
+        if (do_abs_bp)
+            avgerr_abs_bp(irho, idelta)      = 1 - mean(results(iparam).relerr_abs_bp);
+            avgerr_abs_bp(avgerr_abs_bp < 0) = 0;
+        end
+        if (do_gap)
+            avgerr_gap(irho, idelta)         = 1 - mean(results(iparam).relerr_gap);
+            avgerr_gap(avgerr_gap < 0) = 0;
+        end
+        if (do_nesta)
+            avgerr_nesta(irho, idelta)       = 1 - mean(results(iparam).relerr_nesta);
+            avgerr_nesta(avgerr_nesta < 0) = 0;
+        end
+        
+        iparam = iparam + 1;
+    end
+end
+
+%--------------------------------
+% Plot
+%--------------------------------
+show_phasetrans = @show_phasetrans_win;
+iptsetpref('ImshowAxesVisible', 'on');
+close all
+figbase = 'figs/avgerr_SNR20_lbd1e-2_';
+do_save = 1;
+%
+if (do_abs_ompk)
+    figure;
+    %h = show_phasetrans(nexact_abs_ompk, numvects);
+    %bar(0.001:0.001:0.1, hist_abs_ompk);
+    h = show_phasetrans(avgerr_abs_ompk, 1);
+    if do_save
+        figname = [figbase 'ABS_OMPk'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
+%
+if (do_abs_ompeps)
+    figure;
+    %h = show_phasetrans(nexact_abs_ompeps, numvects);
+    %bar(0.001:0.001:0.1, hist_abs_ompeps);
+    h = show_phasetrans(avgerr_abs_ompeps, 1);
+    if do_save
+        figname = [figbase 'ABS_OMPeps'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
+%
+if (do_abs_tst)
+    figure;
+    %h = show_phasetrans(nexact_abs_tst, numvects);
+    %bar(0.001:0.001:0.1, hist_abs_tst);
+    h = show_phasetrans(avgerr_abs_tst, 1);
+    if do_save
+        figname = [figbase 'ABS_TST'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
+%
+if (do_abs_bp)
+    figure;
+    %h = show_phasetrans(nexact_abs_bp, numvects);
+    %bar(0.001:0.001:0.1, hist_abs_bp);
+    h = show_phasetrans(avgerr_abs_bp, 1);
+    if do_save
+        figname = [figbase 'ABS_BP'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
+%
+if (do_gap)
+    figure;
+    %h = show_phasetrans(nexact_gap, numvects);
+    %bar(0.001:0.001:0.1, hist_gap);
+    h = show_phasetrans(avgerr_gap, 1);
+    if do_save
+        figname = [figbase 'GAP'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
+%
+if (do_nesta)
+    figure;
+    %h = show_phasetrans(nexact_nesta, numvects);
+    %bar(0.001:0.001:0.1, hist_nesta);
+    h = show_phasetrans(avgerr_nesta, 1);
+    if do_save
+        figname = [figbase 'NESTA'];
+        saveas(h, [figname '.fig']);
+        saveas(h, [figname '.png']);
+        saveTightFigure(h, [figname '.pdf']);
+    end
+end
\ No newline at end of file