nikcleju@10: # -*- coding: utf-8 -*-
nikcleju@10: """
nikcleju@10: Created on Sat Nov 05 18:08:40 2011
nikcleju@10: 
nikcleju@10: @author: Nic
nikcleju@10: """
nikcleju@10: 
nikcleju@19: import numpy as np
nikcleju@10: import pyCSalgos
nikcleju@19: import pyCSalgos.GAP.GAP
nikcleju@19: import pyCSalgos.SL0.SL0_approx
nikcleju@10: 
nikcleju@19: # Define functions that prepare arguments for each algorithm call
nikcleju@10: def gap_paramsetup(y,M,Omega,epsilon,lbd):
nikcleju@19:   gapparams = {"num_iteration" : 1000,\
nikcleju@19:                    "greedy_level" : 0.9,\
nikcleju@19:                    "stopping_coefficient_size" : 1e-4,\
nikcleju@19:                    "l2solver" : 'pseudoinverse',\
nikcleju@19:                    "noise_level": epsilon}
nikcleju@19:   return y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1])
nikcleju@19: def sl0_paramsetup(y,M,Omega,epsilon,lbd):
nikcleju@19:   
nikcleju@19:   N,n = Omega.shape
nikcleju@19:   D = np.linalg.pinv(Omega)
nikcleju@19:   U,S,Vt = np.linalg.svd(D)
nikcleju@19:   aggDupper = np.dot(M,D)
nikcleju@19:   aggDlower = Vt[-(N-n):,:]
nikcleju@19:   aggD = np.concatenate((aggDupper, lbd * aggDlower))
nikcleju@19:   aggy = np.concatenate((y, np.zeros(N-n)))
nikcleju@19:   
nikcleju@19:   sigmamin = 0.1
nikcleju@19:   return aggD,aggy,epsilon,sigmamin
nikcleju@10: 
nikcleju@19: # Define tuples (algorithm setup function, algorithm function, name)
nikcleju@19: gap = (gap_paramsetup, pyCSalgos.GAP.GAP.GAP, 'GAP')
nikcleju@19: sl0 = (sl0_paramsetup, pyCSalgos.SL0.SL0_approx.SL0_approx, 'SL0_approx')
nikcleju@10: 
nikcleju@19: # Main function
nikcleju@19: def mainrun():
nikcleju@10: 
nikcleju@19:   # Define which algorithms to run
nikcleju@19:   algos = (gap, sl0)
nikcleju@19:   numalgos = len(algos)
nikcleju@19:   
nikcleju@19:   # Set up experiment parameters
nikcleju@19:   sigma = 2.0;
nikcleju@19:   delta = 0.8;
nikcleju@19:   rho   = 0.15;
nikcleju@19:   numvects = 2; # Number of vectors to generate
nikcleju@19:   SNRdb = 20;    # This is norm(signal)/norm(noise), so power, not energy
nikcleju@10: 
nikcleju@19:   # Process parameters
nikcleju@19:   noiselevel = 1 / (10^(SNRdb/10));
nikcleju@19:   d = 50;
nikcleju@19:   p = round(sigma*d);
nikcleju@19:   m = round(delta*d);
nikcleju@19:   l = round(d - rho*m);
nikcleju@19:   
nikcleju@19:   # Generate Omega and data based on parameters
nikcleju@19:   Omega = pyCSalgos.GAP.GAP.Generate_Analysis_Operator(d, p);
nikcleju@19:   # Optionally make Omega more coherent
nikcleju@19:   #[U, S, Vt] = np.linalg.svd(Omega);
nikcleju@19:   #Sdnew = np.diag(S) * (1+np.arange(np.diag(S).size)); % Make D coherent, not Omega!
nikcleju@19:   #Snew = [diag(Sdnew); zeros(size(S,1) - size(S,2), size(S,2))];
nikcleju@19:   #Omega = U * Snew * V';
nikcleju@10: 
nikcleju@19:   # Generate data  
nikcleju@19:   x0,y,M,Lambda,realnoise = pyCSalgos.GAP.GAP.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
nikcleju@19: 
nikcleju@19:   # Values for lambda
nikcleju@19:   #lambdas = [0 10.^linspace(-5, 4, 10)];
nikcleju@19:   lambdas = np.concatenate((np.array([0]), 10**np.linspace(-5, 4, 10)))
nikcleju@10:   
nikcleju@19:   xrec = dict()
nikcleju@19:   err = dict()
nikcleju@19:   relerr = dict()
nikcleju@19:   for i,algo in zip(np.arange(numalgos),algos):
nikcleju@19:     xrec[algo[2]]   = np.zeros((lambdas.size, d, y.shape[1]))
nikcleju@19:     err[algo[2]]    = np.zeros((lambdas.size, y.shape[1]))
nikcleju@19:     relerr[algo[2]] = np.zeros((lambdas.size, y.shape[1]))
nikcleju@19:   
nikcleju@19:   for ilbd,lbd in zip(np.arange(lambdas.size),lambdas):
nikcleju@19:     for iy in np.arange(y.shape[1]):
nikcleju@19:       for algosetupfunc,algofunc,strname in algos:
nikcleju@19:         epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
nikcleju@19:         
nikcleju@19:         inparams = algosetupfunc(y[:,iy],M,Omega,epsilon,lbd)
nikcleju@19:         xrec[strname][ilbd,:,iy] = algofunc(*inparams)[0]
nikcleju@19:         
nikcleju@19:         err[strname][ilbd,iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy])
nikcleju@19:         relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / np.linalg.norm(x0[:,iy])
nikcleju@19:         
nikcleju@19:     print 'Lambda = ',lbd,' :'
nikcleju@19:     for strname in relerr:
nikcleju@19:       print '   ',strname,' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
nikcleju@10: 
nikcleju@10: 
nikcleju@10: 
nikcleju@19: # Script main
nikcleju@19: if __name__ == "__main__":
nikcleju@19:   mainrun()