Mercurial > hg > absrec
changeset 17:7fdf964f4edd
Added docstrings to files and functions
| author | Nic Cleju <nikcleju@gmail.com> |
|---|---|
| date | Tue, 03 Apr 2012 16:27:18 +0300 |
| parents | 23e9b536ba71 |
| children | 4a967f4f18a0 |
| files | ABSexact.py ABSlambda.py ABSmixed.py AnalysisGenerate.py SL0_pt_noise.py algos.py runbatch.py stdparams_approx.py stdparams_exact.py test_approx.py test_exact.py test_exact_old.py utils.py |
| diffstat | 13 files changed, 302 insertions(+), 492 deletions(-) [+] |
line wrap: on
line diff
--- a/ABSexact.py Tue Apr 03 10:18:11 2012 +0300 +++ b/ABSexact.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,19 +1,29 @@ # -*- coding: utf-8 -*- """ -Created on Fri Mar 09 14:06:13 2012 +Algorithms for exact analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Exact reconstruction. -@author: ncleju +Author: Nicolae Cleju """ +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + import numpy + +# Import synthesis solvers from pyCSalgos package import pyCSalgos.BP.l1eq_pd import pyCSalgos.BP.cvxopt_lp import pyCSalgos.OMP.omp_QR import pyCSalgos.SL0.SL0 import pyCSalgos.TST.RecommendedTST + def bp(y,M,Omega,x0, pdtol=1e-3, pdmaxiter=50, cgtol=1e-8, cgmaxiter=200, verbose=False): - + """ + ABS-exact: Basis Pursuit (based on l1magic toolbox) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -26,7 +36,9 @@ return numpy.dot(D , pyCSalgos.BP.l1eq_pd.l1eq_pd(x0,Atilde,Atilde.T,ytilde, pdtol, pdmaxiter, cgtol, cgmaxiter, verbose)) def bp_cvxopt(y,M,Omega): - + """ + ABS-exact: Basis Pursuit (based cvxopt) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -40,6 +52,9 @@ def ompeps(y,M,Omega,epsilon): + """ + ABS-exact: OMP with stopping criterion residual < epsilon + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) @@ -56,6 +71,9 @@ return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0]) def ompk(y,M,Omega,k): + """ + ABS-exact: OMP with stopping criterion fixed number of atoms = k + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) @@ -71,6 +89,9 @@ return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0]) def sl0(y,M,Omega, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, true_s=None): + """ + ABS-exact: Smooth L0 (SL0) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) @@ -84,6 +105,9 @@ return numpy.dot(D, pyCSalgos.SL0.SL0.SL0(Atilde,ytilde,sigma_min,sigma_decrease_factor,mu_0,L,true_s)) def tst_recom(y,M,Omega, nsweep=300, tol=0.00001, xinitial=None, ro=None): + """ + ABS-exact: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega)
--- a/ABSlambda.py Tue Apr 03 10:18:11 2012 +0300 +++ b/ABSlambda.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,18 +1,27 @@ # -*- coding: utf-8 -*- """ -Created on Fri Mar 09 14:06:13 2012 +Algorithms for approximate analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Approximate reconstruction, ABS-lambda. -@author: ncleju +Author: Nicolae Cleju """ +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + import numpy + +# Import synthesis solvers from pyCSalgos package import pyCSalgos.BP.l1qc import pyCSalgos.SL0.SL0_approx import pyCSalgos.OMP.omp_QR import pyCSalgos.TST.RecommendedTST def sl0(y,M,Omega,epsilon,lbd,sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None): - + """ + ABS-lambda: Smooth L0 (SL0) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -24,7 +33,9 @@ return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s)) def bp(y,M,Omega,epsilon,lbd, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False): - + """ + ABS-lambda: Basis Pursuit (based on l1magic toolbox) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -36,7 +47,9 @@ return numpy.dot(D, pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon, lbtol, mu, cgtol, cgmaxiter, verbose)) def ompeps(y,M,Omega,epsilon,lbd): - + """ + ABS-lambda: OMP with stopping criterion residual < epsilon + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -51,7 +64,9 @@ return numpy.dot(D, pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]) def tst_recom(y,M,Omega,epsilon,lbd, nsweep=300, xinitial=None, ro=None): - + """ + ABS-lambda: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D)
--- a/ABSmixed.py Tue Apr 03 10:18:11 2012 +0300 +++ b/ABSmixed.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,16 +1,25 @@ # -*- coding: utf-8 -*- """ -Created on Fri Mar 09 14:06:13 2012 +Algorithms for approximate analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS). +Approximate reconstruction, ABS-mixed. -@author: ncleju +Author: Nicolae Cleju """ +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + import numpy + +# Import synthesis solvers from pyCSalgos package import pyCSalgos.BP.l1qec import pyCSalgos.SL0.SL0_approx def bp(y,M,Omega,epsilon, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False): - + """ + ABS-mixed: Basis Pursuit (based on l1magic toolbox) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D) @@ -20,7 +29,9 @@ return numpy.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,numpy.zeros(N-n), lbtol, mu, cgtol, cgmaxiter, verbose)) def sl0(y,M,Omega,epsilon, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, Aeps_pinv=None, Aexact_pinv=None, true_s=None): - + """ + ABS-mixed: Smooth L0 (SL0) + """ N,n = Omega.shape D = numpy.linalg.pinv(Omega) U,S,Vt = numpy.linalg.svd(D)
--- a/AnalysisGenerate.py Tue Apr 03 10:18:11 2012 +0300 +++ b/AnalysisGenerate.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,9 +1,14 @@ # -*- coding: utf-8 -*- """ -Created on Thu Dec 08 10:56:56 2011 +Generate analysis operator and cosparse data. +Based on Matlab functions writteh by Sangnam Nam, INRIA Rennes -@author: ncleju +Author: Nicolae Cleju """ +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" + import numpy import numpy.linalg @@ -12,7 +17,9 @@ rng = RandomState() def Generate_Analysis_Operator(d, p): - # generate random tight frame with equal column norms + """ + Generate random tight frame with equal column norms + """ if p == d: T = rng.randn(d,d); [Omega, discard] = numpy.qr(T); @@ -37,23 +44,23 @@ return Omega +#function [x0,y,M,LambdaMat] = Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr) def Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr): - #function [x0,y,M,LambdaMat] = Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr) - - # Building an analysis problem, which includes the ingredients: - # - Omega - the analysis operator of size p*d - # - M is anunderdetermined measurement matrix of size m*d (m<d) - # - x0 is a vector of length d that satisfies ||Omega*x0||=p-k - # - Lambda is the true location of these k zeros in Omega*x0 - # - a measurement vector y0=Mx0 is computed - # - noise contaminated measurement vector y is obtained by - # y = y0 + n where n is an additive gaussian noise with norm(n,2)/norm(y0,2) = noiselevel - # Added by Nic: - # - Omega = analysis operator - # - normstr: if 'l0', generate l0 sparse vector (unchanged). If 'l1', - # generate a vector of Laplacian random variables (gamma) and - # pseudoinvert to find x - + """ + Building an analysis problem, which includes the ingredients: + - Omega - the analysis operator of size p*d + - M is anunderdetermined measurement matrix of size m*d (m<d) + - x0 is a vector of length d that satisfies ||Omega*x0||=p-k + - Lambda is the true location of these k zeros in Omega*x0 + - a measurement vector y0=Mx0 is computed + - noise contaminated measurement vector y is obtained by + y = y0 + n where n is an additive gaussian noise with norm(n,2)/norm(y0,2) = noiselevel + Added by Nicolae Cleju: + - Omega = analysis operator + - normstr: if 'l0', generate l0 sparse vector (unchanged). If 'l1', + generate a vector of Laplacian random variables (gamma) and + pseudoinvert to find x + """ # Omega is known as input parameter #Omega=Generate_Analysis_Operator(d, p); # Omega = randn(p,d);
--- a/SL0_pt_noise.py Tue Apr 03 10:18:11 2012 +0300 +++ b/SL0_pt_noise.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,8 +1,8 @@ # -*- coding: utf-8 -*- """ -Created on Sat Nov 05 18:08:40 2011 +Test SL0 behavior in the presence of strong noise -@author: Nic +Author: Nicolae Cleju """ from algos import *
--- a/algos.py Tue Apr 03 10:18:11 2012 +0300 +++ b/algos.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,24 +1,28 @@ # -*- coding: utf-8 -*- """ -Define simple wrappers for algorithms +Define simple wrappers for algorithms, with similar header. +Specific algorithm parameters are defined inside here. -Created on Wed Dec 07 14:06:13 2011 +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" -@author: ncleju -""" import numpy + +# Module with algorithms implemented in Python import pyCSalgos import pyCSalgos.GAP.GAP + +# Analysis by Synthesis - exact algorithms import ABSexact +# Analysis by Synthesis - mixed algorithms import ABSmixed +# Analysis by Synthesis - lambda algorithms import ABSlambda -#import pyCSalgos.BP.l1qc -#import pyCSalgos.BP.l1qec -#import pyCSalgos.SL0.SL0_approx -#import pyCSalgos.OMP.omp_QR -#import pyCSalgos.RecomTST.RecommendedTST -#import pyCSalgos.NESTA.NESTA + ###--------------------------------- ### Exact reconstruction algorithms @@ -161,7 +165,8 @@ return ABSlambda.tst_recom(y,M,Omega,epsilon,lbd, nsweep) -### Algorithm tuples +### Define algorithm tuples: (function, name) +### Will be used in stdparams and in test scripts ## Exact recovery exact_gap = (run_exact_gap, 'GAP') exact_bp = (run_exact_bp, 'ABSexact_BP') @@ -181,124 +186,3 @@ lambda_bp = (run_lambda_bp, 'ABSlambda_BP') lambda_ompeps = (run_lambda_ompeps, 'ABSlambda_OMPeps') lambda_tst = (run_lambda_tst, 'ABSlambda_TST') - -##========================== -## 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,numpy.zeros(Omega.shape[1]))[0] -# -#def run_bp_analysis(y,M,Omega,epsilon): -# -# N,n = Omega.shape -# D = numpy.linalg.pinv(Omega) -# U,S,Vt = numpy.linalg.svd(D) -# Aeps = numpy.dot(M,D) -# Aexact = Vt[-(N-n):,:] -# # We don't ned any aggregate matrices anymore -# -# x0 = numpy.zeros(N) -# return numpy.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,numpy.zeros(N-n))) -# -#def run_sl0_analysis(y,M,Omega,epsilon): -# -# N,n = Omega.shape -# D = numpy.linalg.pinv(Omega) -# U,S,Vt = numpy.linalg.svd(D) -# Aeps = numpy.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 numpy.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 = numpy.linalg.svd(M, full_matrices = True) -# V = V.T # Make like Matlab -# m,n = M.shape # Make like Matlab -# S = numpy.hstack((numpy.diag(S), numpy.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 = numpy.linalg.pinv(Omega) -# #U,S,Vt = numpy.linalg.svd(D) -# aggDupper = numpy.dot(M,D) -# aggDlower = Vt[-(N-n):,:] -# aggD = numpy.concatenate((aggDupper, lbd * aggDlower)) -# aggy = numpy.concatenate((y, numpy.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 = numpy.linalg.pinv(Omega) -# #U,S,Vt = numpy.linalg.svd(D) -# aggDupper = numpy.dot(M,D) -# aggDlower = Vt[-(N-n):,:] -# aggD = numpy.concatenate((aggDupper, lbd * aggDlower)) -# aggy = numpy.concatenate((y, numpy.zeros(N-n))) -# -# x0 = numpy.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 = numpy.linalg.pinv(Omega) -# #U,S,Vt = numpy.linalg.svd(D) -# aggDupper = numpy.dot(M,D) -# aggDlower = Vt[-(N-n):,:] -# aggD = numpy.concatenate((aggDupper, lbd * aggDlower)) -# aggy = numpy.concatenate((y, numpy.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 = numpy.linalg.pinv(Omega) -# #U,S,Vt = numpy.linalg.svd(D) -# aggDupper = numpy.dot(M,D) -# aggDlower = Vt[-(N-n):,:] -# aggD = numpy.concatenate((aggDupper, lbd * aggDlower)) -# aggy = numpy.concatenate((y, numpy.zeros(N-n))) -# -# nsweep = 300 -# tol = epsilon / numpy.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
--- a/runbatch.py Tue Apr 03 10:18:11 2012 +0300 +++ b/runbatch.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- """ -Created on Sat Nov 05 18:08:40 2011 +Trying to write a common framework for simulations +Not finished yet, better be ignored @author: Nic """
--- a/stdparams_approx.py Tue Apr 03 10:18:11 2012 +0300 +++ b/stdparams_approx.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,25 +1,22 @@ # -*- coding: utf-8 -*- """ -Created on Wed Dec 07 14:04:40 2011 +Defines standard parameters for approximate reconstruction simulation -@author: ncleju +Author: Nicolae Cleju """ +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" import numpy + +# Solver algorithms to run 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 +# Test parameters paramstest = dict() paramstest['algosN'] = nesta, # tuple of algorithms not depending on lambda - #algosL = sl0,bp # tuple of algorithms depending on lambda (our ABS approach) -paramstest['algosL'] = lambda_sl0, +paramstest['algosL'] = lambda_sl0, # tuple of algorithms depending on lambda (ABS-lambda) paramstest['d'] = 50.0 paramstest['sigma'] = 2.0 paramstest['deltas'] = numpy.array([0.05, 0.45, 0.95]) @@ -27,7 +24,7 @@ #deltas = numpy.array([0.95]) #deltas = numpy.arange(0.05,1.,0.05) #rhos = numpy.array([0.05]) -paramstest['numvects'] = 10; # Number of vectors to generate +paramstest['numvects'] = 10; # Number of vectors to generate paramstest['SNRdb'] = 20.; # This is norm(signal)/norm(noise), so power, not energy # Values for lambda #lambdas = [0 10.^linspace(-5, 4, 10)]; @@ -39,240 +36,73 @@ # 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,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda - algosL = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (our ABS approach) - - d = 50.0; - sigma = 2.0 - deltas = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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') +# d = 50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000 +# Noise 40db +params1 = dict() +params1['algosN'] = gap,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda +params1['algosL'] = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (ABS-lambda) +params1['d'] = 50.0 +params1['sigma'] = 2.0 +params1['deltas'] = numpy.arange(0.05,1.,0.05) +params1['rhos'] = numpy.arange(0.05,1.,0.05) +params1['numvects'] = 10; # Number of vectors to generate +params1['SNRdb'] = 40.; # This is norm(signal)/norm(noise), so power, not energy +params1['lambdas'] = numpy.array([0., 0.0001, 0.01, 1, 100, 10000]) +params1['savedataname'] = 'approx_pt_stdtest.mat' +params1['saveplotbase'] = 'approx_pt_stdtest_' +params1['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,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda - algosL = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (our ABS approach) +# d = 20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000 +# Noise 40db +params2 = dict() +params2['algosN'] = gap,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda +params2['algosL'] = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (ABS-lambda) +params2['d'] = 20.0 +params2['sigma'] = 10.0 +params2['deltas'] = numpy.arange(0.05,1.,0.05) +params2['rhos'] = numpy.arange(0.05,1.,0.05) +params2['numvects'] = 10; # Number of vectors to generate +params2['SNRdb'] = 40.; # This is norm(signal)/norm(noise), so power, not energy +params2['lambdas'] = numpy.array([0., 0.0001, 0.01, 1, 100, 10000]) +params2['savedataname'] = 'approx_pt_stdtest.mat' +params2['saveplotbase'] = 'approx_pt_stdtest_' +params2['saveplotexts'] = ('png','pdf','eps') - d = 20.0 - sigma = 10.0 - deltas = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 +# 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,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda - algosL = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (our ABS approach) - - d = 50.0; - sigma = 2.0 - deltas = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 +# d = 50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000 +# Identical with 1 but with 20dB SNR noise +params3 = dict() +params3['algosN'] = gap,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda +params3['algosL'] = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (ABS-lambda) +params3['d'] = 50.0 +params3['sigma'] = 2.0 +params3['deltas'] = numpy.arange(0.05,1.,0.05) +params3['rhos'] = numpy.arange(0.05,1.,0.05) +params3['numvects'] = 10; # Number of vectors to generate +params3['SNRdb'] = 20.; # This is norm(signal)/norm(noise), so power, not energy +params3['lambdas'] = numpy.array([0., 0.0001, 0.01, 1, 100, 10000]) +params3['savedataname'] = 'approx_pt_stdtest.mat' +params3['saveplotbase'] = 'approx_pt_stdtest_' +params3['saveplotexts'] = ('png','pdf','eps') # 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,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda - algosL = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (our ABS approach) - - d = 20.0 - sigma = 10.0 - deltas = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 = numpy.arange(0.05,1.,0.05) - rhos = numpy.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 = numpy.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 +# d = 20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000 +# Identical to 2 but with 20dB SNR noise +params4 = dict() +params4['algosN'] = gap,mixed_sl0,mixed_bp,nesta # tuple of algorithms not depending on lambda +params4['algosL'] = lambda_sl0,lambda_bp,lambda_ompeps,lambda_tst # tuple of algorithms depending on lambda (ABS-lambda) +params4['d'] = 20.0 +params4['sigma'] = 10.0 +params4['deltas'] = numpy.arange(0.05,1.,0.05) +params4['rhos'] = numpy.arange(0.05,1.,0.05) +params4['numvects'] = 10; # Number of vectors to generate +params4['SNRdb'] = 20.; # This is norm(signal)/norm(noise), so power, not energy +params4['lambdas'] = numpy.array([0., 0.0001, 0.01, 1, 100, 10000]) +params4['savedataname'] = 'approx_pt_stdtest.mat' +params4['saveplotbase'] = 'approx_pt_stdtest_' +params4['saveplotexts'] = ('png','pdf','eps')
--- a/stdparams_exact.py Tue Apr 03 10:18:11 2012 +0300 +++ b/stdparams_exact.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,13 +1,20 @@ # -*- coding: utf-8 -*- """ -Created on Wed Dec 07 14:04:40 2011 +Defines standard parameters for exact reconstruction simulation +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" -@author: ncleju -""" import numpy + +# Solver algorithms to run from algos import * + +# Test parameters paramstest = dict() paramstest['algos'] = exact_gap,exact_sl0,exact_bp,exact_ompeps,exact_tst # tuple of algorithms #paramstest['algos'] = exact_bp_cvxopt, # tuple of algorithms @@ -26,11 +33,11 @@ # 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 +# d = 200, sigma = 1.2, delta and rho full resolution (0.05 step) +# Virtually no noise (100db) params1 = dict() params1['algos'] = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst # tuple of algorithms -params1['d'] = 50.0; +params1['d'] = 200.0; params1['sigma'] = 1.2 params1['deltas'] = numpy.arange(0.05,1.,0.05) params1['rhos'] = numpy.arange(0.05,1.,0.05) @@ -43,8 +50,8 @@ # 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 +# d = 20, sigma = 10, delta and rho full resolution (0.05 step) +# Virtually no noise (100db) params2 = dict() params2['algos'] = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst # tuple of algorithms params2['d'] = 20.0
--- a/test_approx.py Tue Apr 03 10:18:11 2012 +0300 +++ b/test_approx.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,16 +1,22 @@ # -*- coding: utf-8 -*- """ -Created on Sat Nov 05 18:08:40 2011 +Main script for approximate reconstruction tests. +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" -@author: Nic -""" import numpy import scipy.io import math import os import time +import multiprocessing +import sys +# Try to do smart importing of matplotlib try: import matplotlib if os.name == 'nt': @@ -26,14 +32,15 @@ print "Importing matplotlib.pyplot failed. No figures at all" print "Try selecting a different backend" -import multiprocessing -import sys currmodule = sys.modules[__name__] -# Lock for printing in a thread-safe way -printLock = None +printLock = None # Lock for printing in a thread-safe way + # Thread-safe variable to store number of finished tasks proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array) +# Contains pre-defined simulation parameters import stdparams_approx + +# Analysis operator and data generation functions import AnalysisGenerate # For exceptions @@ -41,6 +48,8 @@ import pyCSalgos.BP.l1qc import pyCSalgos.NESTA.NESTA + + def initProcess(share, ntasks, printLock): """ Pool initializer function (multiprocessing) @@ -55,7 +64,7 @@ def generateTaskParams(globalparams): """ - Generate a list of task parameters + Generate a list of task parameters (for parallel running) """ taskparams = [] SNRdb = globalparams['SNRdb'] @@ -208,7 +217,8 @@ def plot(savedataname): """ - Plot results + Plot results from a mat file. + The files are saved in the current folder. """ params, procresults = loadSim(savedataname) meanmatrix = procresults['meanmatrix'] @@ -234,9 +244,12 @@ # Main function #========================== def run(params): - + """ + Run simulation with given parameters + """ + print "This is analysis recovery ABS approximation script by Nic" - print "Running phase transition ( run_multi() )" + print "Running simulation" algosN = params['algosN'] algosL = params['algosL'] @@ -295,16 +308,30 @@ print "Finished." def run_once_tuple(t): + """ + Wrapper for run_once() that explodes the tuple argument t and shows + the number of finished / remaining tasks + """ + + # Call run_once() here results = run_once(*t) - import sys - currmodule = sys.modules[__name__] - currmodule.proccount.value = currmodule.proccount.value + 1 - print "================================" - print "Finished task",currmodule.proccount.value,"of",currmodule.ntasks - print "================================" + + if currmodule.printLock: + currmodule.printLock.acquire() + + currmodule.proccount.value = currmodule.proccount.value + 1 + print "================================" + print "Finished task",currmodule.proccount.value,"of",currmodule.ntasks + print "================================" + + currmodule.printLock.release() + return results def run_once(algosN,algosL,Omega,y,lambdas,realnoise,M,x0): + """ + Run single task (i.e. task function) + """ d = Omega.shape[1] @@ -374,39 +401,20 @@ return mrelerrN,mrelerrL,elapsed - -def runsingleexampledebug(): - d = 50.0; - sigma = 2.0 - delta = 0.9 - rho = 0.05 - numvects = 20; # Number of vectors to generate - SNRdb = 7.; # This is norm(signal)/norm(noise), so power, not energy - lbd = 10000 - - Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb) - D = numpy.linalg.pinv(Omega) - U,S,Vt = numpy.linalg.svd(D) +def generateFig(): + """ + Generates figures. + Figures are saved in the current folder. + """ + run(stdparams_exact.std1) + plot(stdparams_exact.std1['savedataname']) - xrec = numpy.zeros((d, y.shape[1])) - err = numpy.zeros((y.shape[1])) - relerr = numpy.zeros((y.shape[1])) - - for iy in numpy.arange(y.shape[1]): - epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy]) - gamma = run_sl0(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd) - xrec[:,iy] = numpy.dot(D,gamma) - err[iy] = numpy.linalg.norm(x0[:,iy] - xrec[:,iy]) - relerr[iy] = err[iy] / numpy.linalg.norm(x0[:,iy]) - - print "Finished runsingleexampledebug()" - # Script main if __name__ == "__main__": - #import cProfile - #cProfile.run('mainrun()', 'profile') - #runsingleexampledebug() + + generateFig() - stdparams_approx.paramstest['ncpus'] = 1 - run(stdparams_approx.paramstest) - plot(stdparams_approx.paramstest['savedataname']) + # Run test parameters + #stdparams_approx.paramstest['ncpus'] = 1 + #run(stdparams_approx.paramstest) + #plot(stdparams_approx.paramstest['savedataname'])
--- a/test_exact.py Tue Apr 03 10:18:11 2012 +0300 +++ b/test_exact.py Tue Apr 03 16:27:18 2012 +0300 @@ -1,16 +1,22 @@ # -*- coding: utf-8 -*- """ -Created on Sat Nov 05 18:08:40 2011 +Main script for exact reconstruction tests. +Author: Nicolae Cleju +""" +__author__ = "Nicolae Cleju" +__license__ = "GPL" +__email__ = "nikcleju@gmail.com" -@author: Nic -""" import numpy import scipy.io import math import os import time +import multiprocessing +import sys +# Try to do smart importing of matplotlib try: import matplotlib if os.name == 'nt': @@ -26,15 +32,15 @@ print "Importing matplotlib.pyplot failed. No figures at all" print "Try selecting a different backend" -import multiprocessing -import sys currmodule = sys.modules[__name__] -# Lock for printing in a thread-safe way -printLock = None -# Thread-safe variable to store number of finished tasks +printLock = None # Lock for printing in a thread-safe way + # Thread-safe variable to store number of finished tasks currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array) +# Contains pre-defined simulation parameters import stdparams_exact + +# Analysis operator and data generation functions import AnalysisGenerate # For exceptions @@ -42,6 +48,7 @@ import pyCSalgos.NESTA.NESTA + def initProcess(share, ntasks, printLock): """ Pool initializer function (multiprocessing) @@ -57,7 +64,7 @@ def generateTaskParams(globalparams): """ - Generate a list of task parameters + Generate a list of task parameters (for parallel running) """ taskparams = [] SNRdb = globalparams['SNRdb'] @@ -184,7 +191,8 @@ def plot(savedataname): """ - Plot results + Plot results from a mat file. + The files are saved in the current folder. """ params, procresults = loadSim(savedataname) meanmatrix = procresults['meanmatrix'] @@ -204,11 +212,11 @@ #========================== def run(params): """ - Run with given parameters + Run simulation with given parameters """ - print "This is analysis recovery ABS approximation script by Nic" - print "Running phase transition ( run_multi() )" + print "This is analysis recovery ABS exact script by Nic" + print "Running simulation algos = params['algos'] d = params['d'] @@ -265,6 +273,12 @@ print "Finished." def run_once_tuple(t): + """ + Wrapper for run_once() that explodes the tuple argument t and shows + the number of finished / remaining tasks + """ + + # Call run_once() here results = run_once(*t) if currmodule.printLock: @@ -281,7 +295,7 @@ def run_once(algos,Omega,y,M,x0): """ - Run single task (task function) + Run single task (i.e. task function) """ d = Omega.shape[1] @@ -326,7 +340,6 @@ # mrelerr[algotuple[1]] = numpy.mean(relerr[algotuple[1]]) #return mrelerr,elapsed - # Should return number of reconstructions with error < threshold, not average error exactthr = 1e-6 mrelerr = dict() for algotuple in algos: @@ -335,20 +348,30 @@ def testMatlab(): + """ + For debugging only. + Load parameters from a mat file saved by Matlab. + """ mdict = scipy.io.loadmat("E:\\CS\\Ale mele\\Analysis_ExactRec\\temp.mat") algos = stdparams_exact.std1()[0] res = run_once(algos, mdict['Omega'].byteswap().newbyteorder(),mdict['y'],mdict['M'],mdict['x0']) def generateFig(): + """ + Generates figures from paper "Analysis-based sparse reconstruction with synthesis-based solvers". + The figures are saved in the current folder. + """ run(stdparams_exact.std1) + plot(stdparams_exact.std1['savedataname']) # Script main if __name__ == "__main__": - #import cProfile - #cProfile.run('mainrun()', 'profile') - #run_mp(stdparams_exact.stdtest) - #runsingleexampledebug() - stdparams_exact.paramstest['ncpus'] = 1 - run(stdparams_exact.paramstest) - plot(stdparams_exact.paramstest['savedataname']) + # Set the number of cpus for paraller running (or comment to leave default = max) + #stdparams_exact.paramstest['ncpus'] = 1 + + generateFig() + + # Run test parameters + #run(stdparams_exact.paramstest) + #plot(stdparams_exact.paramstest['savedataname'])