changeset 51:eb4c66488ddf

Split algos.py and stdparams.py, added nesta to std1, 2, 3, 4

--- a/scripts/ABSapprox.py	Thu Dec 01 22:11:19 2011 +0000
+++ b/scripts/ABSapprox.py	Wed Dec 07 12:44:19 2011 +0000
@@ -10,6 +10,8 @@
 import math
 import os
 
+import stdparams
+
 import pyCSalgos
 import pyCSalgos.GAP.GAP
 import pyCSalgos.BP.l1qc
@@ -19,126 +21,6 @@
 import pyCSalgos.RecomTST.RecommendedTST
 import pyCSalgos.NESTA.NESTA
 
-#==========================
-# Algorithm functions
-#==========================
-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_bp_analysis(y,M,Omega,epsilon):
-  
-  N,n = Omega.shape
-  D = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
-  Aeps = np.dot(M,D)
-  Aexact = Vt[-(N-n):,:]
-  # We don't ned any aggregate matrices anymore
-  
-  x0 = np.zeros(N)
-  return np.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,np.zeros(N-n)))
-
-def run_sl0_analysis(y,M,Omega,epsilon):
-  
-  N,n = Omega.shape
-  D = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
-  Aeps = np.dot(M,D)
-  Aexact = Vt[-(N-n):,:]
-  # We don't ned any aggregate matrices anymore
-  
-  sigmamin = 0.001
-  sigma_decrease_factor = 0.5
-  mu_0 = 2
-  L = 10
-  return np.dot(D , pyCSalgos.SL0.SL0_approx.SL0_approx_analysis(Aeps,Aexact,y,epsilon,sigmamin,sigma_decrease_factor,mu_0,L))
-
-def run_nesta(y,M,Omega,epsilon):
-  
-  U,S,V = np.linalg.svd(M, full_matrices = True)
-  V = V.T         # Make like Matlab
-  m,n = M.shape   # Make like Matlab
-  S = np.hstack((np.diag(S), np.zeros((m,n-m))))  
-
-  opt_muf = 1e-3
-  optsUSV = {'U':U, 'S':S, 'V':V}
-  opts = {'U':Omega, 'Ut':Omega.T.copy(), 'USV':optsUSV, 'TolVar':1e-5, 'Verbose':0}
-  return pyCSalgos.NESTA.NESTA.NESTA(M, None, y, opt_muf, epsilon, opts)[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)
-  
-def run_bp(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)))
-
-  x0 = np.zeros(N)
-  return pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon)
-
-def run_ompeps(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)))
-  
-  opts = dict()
-  opts['stopCrit'] = 'mse'
-  opts['stopTol'] = epsilon**2 / aggy.size
-  return pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]
-  
-def run_tst(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)))
-  
-  nsweep = 300
-  tol = epsilon / np.linalg.norm(aggy)
-  return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep=nsweep, tol=tol)
-
-#==========================
-# Define tuples (algorithm function, name)
-#==========================
-gap = (run_gap, 'GAP')
-sl0 = (run_sl0, 'SL0a')
-sl0analysis = (run_sl0_analysis, 'SL0a2')
-bpanalysis = (run_bp_analysis, 'BPa2')
-nesta = (run_nesta, 'NESTA')
-bp  = (run_bp, 'BP')
-ompeps = (run_ompeps, 'OMPeps')
-tst = (run_tst, 'TST')
 
 #==========================
 # Pool initializer function (multiprocessing)
@@ -151,296 +33,18 @@
     currmodule = sys.modules[__name__]
     currmodule.proccount = share
     currmodule.njobs = njobs
-  
-#==========================
-# Standard parameters
-#==========================
-# Standard parameters for quick testing
-# Algorithms: GAP, SL0 and BP
-# d=50, sigma = 2, delta and rho only 3 x 3, lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# Useful for short testing 
-def stdtest():
-  # Define which algorithms to run
-  algosN = nesta,      # tuple of algorithms not depending on lambda
-  #algosL = sl0,bp    # tuple of algorithms depending on lambda (our ABS approach)
-  algosL = ()
-  
-  d = 50.0
-  sigma = 2.0
-  deltas = np.array([0.05, 0.45, 0.95])
-  rhos = np.array([0.05, 0.45, 0.95])
-  #deltas = np.array([0.95])
-  #deltas = np.arange(0.05,1.,0.05)
-  #rhos = np.array([0.05])
-  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.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_stdtest.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_stdtest_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts   
-
-
-# Standard parameters 1
-# All algorithms, 100 vectors
-# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-def std1():
-  # Define which algorithms to run
-  algosN = gap,sl0analysis,bpanalysis               # tuple of algorithms not depending on lambda
-  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 50.0;
-  sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std1.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std1_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-          
-         
-# Standard parameters 2
-# All algorithms, 100 vectors
-# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-def std2():
-  # Define which algorithms to run
-  algosN = gap,sl0analysis,bpanalysis      # tuple of algorithms not depending on lambda
-  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 20.0
-  sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std2.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std2_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-  
-  
-  # Standard parameters 3
-# All algorithms, 100 vectors
-# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# IDENTICAL with 1 but with 10dB SNR noise
-def std3():
-  # Define which algorithms to run
-  algosN = gap,sl0analysis,bpanalysis               # tuple of algorithms not depending on lambda
-  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 50.0;
-  sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # Number of vectors to generate
-  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
-  # Values for lambda
-  #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std3.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std3_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-          
-# Standard parameters 4
-# All algorithms, 100 vectors
-# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# Identical to 2 but with 10dB SNR noise
-def std4():
-  # Define which algorithms to run
-  algosN = gap,sl0analysis,bpanalysis      # tuple of algorithms not depending on lambda
-  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 20.0
-  sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # Number of vectors to generate
-  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
-  # Values for lambda
-  #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std4.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std4_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-          
-# Standard parameters 1nesta
-# Only NESTA, 100 vectors
-# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# Identical to std1 but with only NESTA
-def std1nesta():
-  # Define which algorithms to run
-  algosN = nesta,               # tuple of algorithms not depending on lambda
-  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 50.0;
-  sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std1nesta.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std1nesta_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts    
-          
-# Standard parameters 2nesta
-# Only NESTA, 100 vectors
-# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# Identical with std2, but with only NESTA
-def std2nesta():
-  # Define which algorithms to run
-  algosN = nesta,      # tuple of algorithms not depending on lambda
-  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 20.0
-  sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std2nesta.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std2nesta_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-
-  # Standard parameters 3nesta
-# Only NESTA, 100 vectors
-# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# IDENTICAL with 3 but with only NESTA
-def std3nesta():
-  # Define which algorithms to run
-  algosN = nesta,               # tuple of algorithms not depending on lambda
-  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 50.0;
-  sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # Number of vectors to generate
-  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
-  # Values for lambda
-  #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std3nesta.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std3nesta_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
-          
-# Standard parameters 4nesta
-# Only NESTA, 100 vectors
-# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
-# Do save data, do save plots, don't show plots
-# Identical to 4 but with only NESTA
-def std4nesta():
-  # Define which algorithms to run
-  algosN = nesta,      # tuple of algorithms not depending on lambda
-  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
-  
-  d = 20.0
-  sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
-  numvects = 100; # Number of vectors to generate
-  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
-  # Values for lambda
-  #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_std4nesta.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_std4nesta_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts
           
 #==========================
 # Interface run functions
 #==========================
-def run_mp(std=std2,ncpus=None):
+def run_mp(std=stdparams.std2,ncpus=None):
   
   algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std()
   run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
             doparallel=True, ncpus=ncpus,\
             doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)
 
-def run(std=std2):
+def run(std=stdparams.std2):
   algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std()
   run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
             doparallel=False,\
@@ -726,5 +330,5 @@
 if __name__ == "__main__":
   #import cProfile
   #cProfile.run('mainrun()', 'profile')    
-  run(stdtest)
+  run(stdparams.stdtest)
   #runsingleexampledebug()

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/algos.py	Wed Dec 07 12:44:19 2011 +0000
@@ -0,0 +1,128 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 07 14:06:13 2011
+
+@author: ncleju
+"""
+
+#==========================
+# Algorithm functions
+#==========================
+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_bp_analysis(y,M,Omega,epsilon):
+  
+  N,n = Omega.shape
+  D = np.linalg.pinv(Omega)
+  U,S,Vt = np.linalg.svd(D)
+  Aeps = np.dot(M,D)
+  Aexact = Vt[-(N-n):,:]
+  # We don't ned any aggregate matrices anymore
+  
+  x0 = np.zeros(N)
+  return np.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,np.zeros(N-n)))
+
+def run_sl0_analysis(y,M,Omega,epsilon):
+  
+  N,n = Omega.shape
+  D = np.linalg.pinv(Omega)
+  U,S,Vt = np.linalg.svd(D)
+  Aeps = np.dot(M,D)
+  Aexact = Vt[-(N-n):,:]
+  # We don't ned any aggregate matrices anymore
+  
+  sigmamin = 0.001
+  sigma_decrease_factor = 0.5
+  mu_0 = 2
+  L = 10
+  return np.dot(D , pyCSalgos.SL0.SL0_approx.SL0_approx_analysis(Aeps,Aexact,y,epsilon,sigmamin,sigma_decrease_factor,mu_0,L))
+
+def run_nesta(y,M,Omega,epsilon):
+  
+  U,S,V = np.linalg.svd(M, full_matrices = True)
+  V = V.T         # Make like Matlab
+  m,n = M.shape   # Make like Matlab
+  S = np.hstack((np.diag(S), np.zeros((m,n-m))))  
+
+  opt_muf = 1e-3
+  optsUSV = {'U':U, 'S':S, 'V':V}
+  opts = {'U':Omega, 'Ut':Omega.T.copy(), 'USV':optsUSV, 'TolVar':1e-5, 'Verbose':0}
+  return pyCSalgos.NESTA.NESTA.NESTA(M, None, y, opt_muf, epsilon, opts)[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)
+  
+def run_bp(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)))
+
+  x0 = np.zeros(N)
+  return pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon)
+
+def run_ompeps(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)))
+  
+  opts = dict()
+  opts['stopCrit'] = 'mse'
+  opts['stopTol'] = epsilon**2 / aggy.size
+  return pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]
+  
+def run_tst(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)))
+  
+  nsweep = 300
+  tol = epsilon / np.linalg.norm(aggy)
+  return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep=nsweep, tol=tol)
+
+
+#==========================
+# Define tuples (algorithm function, name)
+#==========================
+gap = (run_gap, 'GAP')
+sl0 = (run_sl0, 'SL0a')
+sl0analysis = (run_sl0_analysis, 'SL0a2')
+bpanalysis = (run_bp_analysis, 'BPa2')
+nesta = (run_nesta, 'NESTA')
+bp  = (run_bp, 'BP')
+ompeps = (run_ompeps, 'OMPeps')
+tst = (run_tst, 'TST')
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/stdparams.py	Wed Dec 07 12:44:19 2011 +0000
@@ -0,0 +1,286 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 07 14:04:40 2011
+
+@author: ncleju
+"""
+
+from algos import *
+
+#==========================
+# Standard parameters
+#==========================
+# Standard parameters for quick testing
+# Algorithms: GAP, SL0 and BP
+# d=50, sigma = 2, delta and rho only 3 x 3, lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Useful for short testing 
+def stdtest():
+  # Define which algorithms to run
+  algosN = nesta,      # tuple of algorithms not depending on lambda
+  #algosL = sl0,bp    # tuple of algorithms depending on lambda (our ABS approach)
+  algosL = ()
+  
+  d = 50.0
+  sigma = 2.0
+  deltas = np.array([0.05, 0.45, 0.95])
+  rhos = np.array([0.05, 0.45, 0.95])
+  #deltas = np.array([0.95])
+  #deltas = np.arange(0.05,1.,0.05)
+  #rhos = np.array([0.05])
+  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.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_stdtest.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_stdtest_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts   
+
+
+# Standard parameters 1
+# All algorithms, 100 vectors
+# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+def std1():
+  # Define which algorithms to run
+  algosN = gap,sl0analysis,bpanalysis,nesta       # tuple of algorithms not depending on lambda
+  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 50.0;
+  sigma = 2.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std1.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std1_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+         
+# Standard parameters 2
+# All algorithms, 100 vectors
+# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+def std2():
+  # Define which algorithms to run
+  algosN = gap,sl0analysis,bpanalysis,nesta      # tuple of algorithms not depending on lambda
+  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 20.0
+  sigma = 10.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std2.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std2_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+  
+  
+  # Standard parameters 3
+# All algorithms, 100 vectors
+# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# IDENTICAL with 1 but with 10dB SNR noise
+def std3():
+  # Define which algorithms to run
+  algosN = gap,sl0analysis,bpanalysis,nesta        # tuple of algorithms not depending on lambda
+  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 50.0;
+  sigma = 2.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std3.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std3_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+# Standard parameters 4
+# All algorithms, 100 vectors
+# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Identical to 2 but with 10dB SNR noise
+def std4():
+  # Define which algorithms to run
+  algosN = gap,sl0analysis,bpanalysis,nesta    # tuple of algorithms not depending on lambda
+  algosL = sl0,bp,ompeps,tst    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 20.0
+  sigma = 10.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std4.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std4_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+# Standard parameters 1nesta
+# Only NESTA, 100 vectors
+# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Identical to std1 but with only NESTA
+def std1nesta():
+  # Define which algorithms to run
+  algosN = nesta,               # tuple of algorithms not depending on lambda
+  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 50.0;
+  sigma = 2.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std1nesta.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std1nesta_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts    
+          
+# Standard parameters 2nesta
+# Only NESTA, 100 vectors
+# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Identical with std2, but with only NESTA
+def std2nesta():
+  # Define which algorithms to run
+  algosN = nesta,      # tuple of algorithms not depending on lambda
+  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 20.0
+  sigma = 10.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # 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.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std2nesta.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std2nesta_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+
+  # Standard parameters 3nesta
+# Only NESTA, 100 vectors
+# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# IDENTICAL with 3 but with only NESTA
+def std3nesta():
+  # Define which algorithms to run
+  algosN = nesta,               # tuple of algorithms not depending on lambda
+  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 50.0;
+  sigma = 2.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std3nesta.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std3nesta_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+# Standard parameters 4nesta
+# Only NESTA, 100 vectors
+# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Identical to 4 but with only NESTA
+def std4nesta():
+  # Define which algorithms to run
+  algosN = nesta,      # tuple of algorithms not depending on lambda
+  algosL = ()    # tuple of algorithms depending on lambda (our ABS approach)
+  
+  d = 20.0
+  sigma = 10.0
+  deltas = np.arange(0.05,1.,0.05)
+  rhos = np.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
+  # Values for lambda
+  #lambdas = [0 10.^linspace(-5, 4, 10)];
+  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  
+  dosavedata = True
+  savedataname = 'approx_pt_std4nesta.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'approx_pt_std4nesta_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
\ No newline at end of file