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@22: import scipy.io
nikcleju@22: import math
nikcleju@29: 
nikcleju@10: import pyCSalgos
nikcleju@19: import pyCSalgos.GAP.GAP
nikcleju@19: import pyCSalgos.SL0.SL0_approx
nikcleju@30: import pyCSalgos.OMP.omp_QR
nikcleju@30: import pyCSalgos.RecomTST.RecommendedTST
nikcleju@10: 
nikcleju@29: #==========================
nikcleju@29: # Algorithm functions
nikcleju@29: #==========================
nikcleju@22: def run_gap(y,M,Omega,epsilon):
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@22:   return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1]))[0]
nikcleju@29: 
nikcleju@22: def run_sl0(y,M,Omega,D,U,S,Vt,epsilon,lbd):
nikcleju@19:   
nikcleju@19:   N,n = Omega.shape
nikcleju@22:   #D = np.linalg.pinv(Omega)
nikcleju@22:   #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@22:   sigmamin = 0.001
nikcleju@22:   sigma_decrease_factor = 0.5
nikcleju@20:   mu_0 = 2
nikcleju@20:   L = 10
nikcleju@22:   return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L)
nikcleju@10: 
nikcleju@27: def run_bp(y,M,Omega,D,U,S,Vt,epsilon,lbd):
nikcleju@27:   
nikcleju@27:   N,n = Omega.shape
nikcleju@27:   #D = np.linalg.pinv(Omega)
nikcleju@27:   #U,S,Vt = np.linalg.svd(D)
nikcleju@27:   aggDupper = np.dot(M,D)
nikcleju@27:   aggDlower = Vt[-(N-n):,:]
nikcleju@27:   aggD = np.concatenate((aggDupper, lbd * aggDlower))
nikcleju@27:   aggy = np.concatenate((y, np.zeros(N-n)))
nikcleju@27:   
nikcleju@27:   sigmamin = 0.001
nikcleju@27:   sigma_decrease_factor = 0.5
nikcleju@27:   mu_0 = 2
nikcleju@27:   L = 10
nikcleju@27:   return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L)
nikcleju@27: 
nikcleju@30: def run_ompeps(y,M,Omega,D,U,S,Vt,epsilon,lbd):
nikcleju@30:   
nikcleju@30:   N,n = Omega.shape
nikcleju@30:   #D = np.linalg.pinv(Omega)
nikcleju@30:   #U,S,Vt = np.linalg.svd(D)
nikcleju@30:   aggDupper = np.dot(M,D)
nikcleju@30:   aggDlower = Vt[-(N-n):,:]
nikcleju@30:   aggD = np.concatenate((aggDupper, lbd * aggDlower))
nikcleju@30:   aggy = np.concatenate((y, np.zeros(N-n)))
nikcleju@30:   
nikcleju@30:   opts = dict()
nikcleju@30:   opts['stopCrit'] = 'mse'
nikcleju@30:   opts['stopTol'] = epsilon**2 / aggy.size
nikcleju@30:   return pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]
nikcleju@30:   
nikcleju@30: def run_tst(y,M,Omega,D,U,S,Vt,epsilon,lbd):
nikcleju@30:   
nikcleju@30:   N,n = Omega.shape
nikcleju@30:   #D = np.linalg.pinv(Omega)
nikcleju@30:   #U,S,Vt = np.linalg.svd(D)
nikcleju@30:   aggDupper = np.dot(M,D)
nikcleju@30:   aggDlower = Vt[-(N-n):,:]
nikcleju@30:   aggD = np.concatenate((aggDupper, lbd * aggDlower))
nikcleju@30:   aggy = np.concatenate((y, np.zeros(N-n)))
nikcleju@30:   
nikcleju@30:   return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep=3000, tol=epsilon / np.linalg.norm(aggy))
nikcleju@30: 
nikcleju@29: #==========================
nikcleju@29: # Define tuples (algorithm function, name)
nikcleju@29: #==========================
nikcleju@22: gap = (run_gap, 'GAP')
nikcleju@30: sl0 = (run_sl0, 'SL0a')
nikcleju@29: bp  = (run_bp, 'BP')
nikcleju@30: ompeps = (run_ompeps, 'OMPeps')
nikcleju@30: tst = (run_tst, 'TST')
nikcleju@10: 
nikcleju@22: # Define which algorithms to run
nikcleju@22: #  1. Algorithms not depending on lambda
nikcleju@22: algosN = gap,   # tuple
nikcleju@22: #  2. Algorithms depending on lambda (our ABS approach)
nikcleju@30: algosL = sl0,bp,ompeps,tst   # tuple
nikcleju@29:   
nikcleju@29: #==========================
nikcleju@29: # Interface functions
nikcleju@29: #==========================
nikcleju@29: def run_multiproc(ncpus=None):
nikcleju@30:   d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = standard_params()
nikcleju@29:   run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
nikcleju@30:             doparallel=True, ncpus=ncpus,\
nikcleju@30:             doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)
nikcleju@22: 
nikcleju@29: def run():
nikcleju@30:   d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = standard_params()
nikcleju@29:   run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
nikcleju@30:             doparallel=False,\
nikcleju@30:             doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)
nikcleju@19:   
nikcleju@29: def standard_params():
nikcleju@29:   #Set up standard experiment parameters
nikcleju@25:   d = 50.0;
nikcleju@22:   sigma = 2.0
nikcleju@27:   #deltas = np.arange(0.05,1.,0.05)
nikcleju@27:   #rhos = np.arange(0.05,1.,0.05)
nikcleju@30:   #deltas = np.array([0.05, 0.45, 0.95])
nikcleju@30:   #rhos = np.array([0.05, 0.45, 0.95])
nikcleju@30:   deltas = np.array([0.05])
nikcleju@30:   rhos = np.array([0.05])
nikcleju@22:   #delta = 0.8;
nikcleju@22:   #rho   = 0.15;
nikcleju@27:   numvects = 100; # Number of vectors to generate
nikcleju@20:   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
nikcleju@22:   # Values for lambda
nikcleju@22:   #lambdas = [0 10.^linspace(-5, 4, 10)];
nikcleju@25:   #lambdas = np.concatenate((np.array([0]), 10**np.linspace(-5, 4, 10)))
nikcleju@25:   lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
nikcleju@29:   
nikcleju@29:   dosavedata = True
nikcleju@30:   savedataname = 'approx_pt_std1.mat'
nikcleju@30:   
nikcleju@30:   doshowplot = False
nikcleju@30:   dosaveplot = True
nikcleju@30:   saveplotbase = 'approx_pt_std1_'
nikcleju@30:   saveplotexts = ('png','pdf','eps')
nikcleju@29:     
nikcleju@29:   
nikcleju@30:   return d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
nikcleju@30:           doshowplot,dosaveplot,saveplotbase,saveplotexts
nikcleju@29:   
nikcleju@29: #==========================
nikcleju@29: # Main functions  
nikcleju@29: #==========================
nikcleju@29: def run_multi(algosN, algosL, d, sigma, deltas, rhos, lambdas, numvects, SNRdb,
nikcleju@29:             doparallel=False, ncpus=None,\
nikcleju@29:             doshowplot=False, dosaveplot=False, saveplotbase=None, saveplotexts=None,\
nikcleju@29:             dosavedata=False, savedataname=None):
nikcleju@30: 
nikcleju@30:   print "This is analysis recovery ABS approximation script by Nic"
nikcleju@30:   print "Running phase transition ( run_multi() )"
nikcleju@29:   
nikcleju@29:   if doparallel:
nikcleju@29:     from multiprocessing import Pool
nikcleju@29:     
nikcleju@30:   if dosaveplot or doshowplot:
nikcleju@30:     try:
nikcleju@30:       import matplotlib
nikcleju@30:       if doshowplot:
nikcleju@30:         print "Importing matplotlib with default (GUI) backend... ",
nikcleju@30:       else:
nikcleju@30:         print "Importing matplotlib with \"Cairo\" backend... ",
nikcleju@30:         matplotlib.use('Cairo')
nikcleju@30:       import matplotlib.pyplot as plt
nikcleju@30:       import matplotlib.cm as cm
nikcleju@30:       print "OK"        
nikcleju@30:     except:
nikcleju@30:       print "FAIL"
nikcleju@30:       print "Importing matplotlib.pyplot failed. No figures at all"
nikcleju@30:       print "Try selecting a different backend"
nikcleju@30:       doshowplot = False
nikcleju@30:       dosaveplot = False
nikcleju@30:   
nikcleju@30:   # Print summary of parameters
nikcleju@30:   print "Parameters:"
nikcleju@30:   if doparallel:
nikcleju@30:     if ncpus is None:
nikcleju@30:       print "  Running in parallel with default threads using \"multiprocessing\" package"
nikcleju@30:     else:
nikcleju@30:       print "  Running in parallel with",ncpus,"threads using \"multiprocessing\" package"
nikcleju@30:   else:
nikcleju@30:     print "Running single thread"
nikcleju@30:   if doshowplot:
nikcleju@30:     print "  Showing figures"
nikcleju@30:   else:
nikcleju@30:     print "  Not showing figures"
nikcleju@30:   if dosaveplot:
nikcleju@30:     print "  Saving figures as "+saveplotbase+"* with extensions ",saveplotexts
nikcleju@30:   else:
nikcleju@30:     print "  Not saving figures"
nikcleju@30:   print "  Running algorithms",[algotuple[1] for algotuple in algosN],[algotuple[1] for algotuple in algosL]
nikcleju@29:   
nikcleju@29:   nalgosN = len(algosN)  
nikcleju@29:   nalgosL = len(algosL)
nikcleju@29:   
nikcleju@22:   meanmatrix = dict()
nikcleju@22:   for i,algo in zip(np.arange(nalgosN),algosN):
nikcleju@22:     meanmatrix[algo[1]]   = np.zeros((rhos.size, deltas.size))
nikcleju@22:   for i,algo in zip(np.arange(nalgosL),algosL):
nikcleju@22:     meanmatrix[algo[1]]   = np.zeros((lambdas.size, rhos.size, deltas.size))
nikcleju@22:   
nikcleju@29:   # Prepare parameters
nikcleju@29:   jobparams = []
nikcleju@30:   print "  (delta, rho) pairs to be run:"
nikcleju@22:   for idelta,delta in zip(np.arange(deltas.size),deltas):
nikcleju@22:     for irho,rho in zip(np.arange(rhos.size),rhos):
nikcleju@22:       
nikcleju@22:       # Generate data and operator
nikcleju@29:       Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb)
nikcleju@22:       
nikcleju@29:       #Save the parameters, and run after
nikcleju@30:       print "    delta = ",delta," rho = ",rho
nikcleju@29:       jobparams.append((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
nikcleju@29: 
nikcleju@30:   print "End of parameters"
nikcleju@30:   
nikcleju@29:   # Run
nikcleju@29:   jobresults = []
nikcleju@29:   if doparallel:
nikcleju@29:     pool = Pool(4)
nikcleju@29:     jobresults = pool.map(run_once_tuple,jobparams)
nikcleju@29:   else:
nikcleju@29:     for jobparam in jobparams:
nikcleju@29:       jobresults.append(run_once(algosN,algosL,Omega,y,lambdas,realnoise,M,x0))
nikcleju@29: 
nikcleju@29:   # Read results
nikcleju@29:   idx = 0
nikcleju@29:   for idelta,delta in zip(np.arange(deltas.size),deltas):
nikcleju@29:     for irho,rho in zip(np.arange(rhos.size),rhos):
nikcleju@29:       mrelerrN,mrelerrL = jobresults[idx]
nikcleju@29:       idx = idx+1
nikcleju@22:       for algotuple in algosN: 
nikcleju@22:         meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[algotuple[1]]
nikcleju@22:         if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]):
nikcleju@22:           meanmatrix[algotuple[1]][irho,idelta] = 0
nikcleju@22:       for algotuple in algosL:
nikcleju@22:         for ilbd in np.arange(lambdas.size):
nikcleju@22:           meanmatrix[algotuple[1]][ilbd,irho,idelta] = 1 - mrelerrL[algotuple[1]][ilbd]
nikcleju@22:           if meanmatrix[algotuple[1]][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][ilbd,irho,idelta]):
nikcleju@22:             meanmatrix[algotuple[1]][ilbd,irho,idelta] = 0
nikcleju@22:    
nikcleju@22:   #  # Prepare matrices to show
nikcleju@22:   #  showmats = dict()
nikcleju@22:   #  for i,algo in zip(np.arange(nalgosN),algosN):
nikcleju@22:   #    showmats[algo[1]]   = np.zeros(rhos.size, deltas.size)
nikcleju@22:   #  for i,algo in zip(np.arange(nalgosL),algosL):
nikcleju@22:   #    showmats[algo[1]]   = np.zeros(lambdas.size, rhos.size, deltas.size)
nikcleju@22: 
nikcleju@22:     # Save
nikcleju@29:     if dosavedata:
nikcleju@29:       tosave = dict()
nikcleju@29:       tosave['meanmatrix'] = meanmatrix
nikcleju@29:       tosave['d'] = d
nikcleju@29:       tosave['sigma'] = sigma
nikcleju@29:       tosave['deltas'] = deltas
nikcleju@29:       tosave['rhos'] = rhos
nikcleju@29:       tosave['numvects'] = numvects
nikcleju@29:       tosave['SNRdb'] = SNRdb
nikcleju@29:       tosave['lambdas'] = lambdas
nikcleju@29:       try:
nikcleju@29:         scipy.io.savemat(savedataname, tosave)
nikcleju@29:       except:
nikcleju@29:         print "Save error"
nikcleju@22:   # Show
nikcleju@29:   if doshowplot or dosaveplot:
nikcleju@27:     for algotuple in algosN:
nikcleju@29:       algoname = algotuple[1]
nikcleju@27:       plt.figure()
nikcleju@29:       plt.imshow(meanmatrix[algoname], cmap=cm.gray, interpolation='nearest',origin='lower')
nikcleju@29:       if dosaveplot:
nikcleju@29:         for ext in saveplotexts:
nikcleju@29:           plt.savefig(saveplotbase + algoname + '.' + ext)
nikcleju@27:     for algotuple in algosL:
nikcleju@29:       algoname = algotuple[1]
nikcleju@27:       for ilbd in np.arange(lambdas.size):
nikcleju@27:         plt.figure()
nikcleju@29:         plt.imshow(meanmatrix[algoname][ilbd], cmap=cm.gray, interpolation='nearest',origin='lower')
nikcleju@29:         if dosaveplot:
nikcleju@29:           for ext in saveplotexts:
nikcleju@30:             plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext)
nikcleju@29:     if doshowplot:
nikcleju@29:       plt.show()
nikcleju@29:     
nikcleju@22:   print "Finished."
nikcleju@22:   
nikcleju@29: def run_once_tuple(t):
nikcleju@29:   return run_once(*t)
nikcleju@10: 
nikcleju@29: def run_once(algosN,algosL,Omega,y,lambdas,realnoise,M,x0):
nikcleju@22:   
nikcleju@22:   d = Omega.shape[1]  
nikcleju@22:   
nikcleju@22:   nalgosN = len(algosN)  
nikcleju@22:   nalgosL = len(algosL)
nikcleju@10:   
nikcleju@19:   xrec = dict()
nikcleju@19:   err = dict()
nikcleju@19:   relerr = dict()
nikcleju@22: 
nikcleju@22:   # Prepare storage variables for algorithms non-Lambda
nikcleju@22:   for i,algo in zip(np.arange(nalgosN),algosN):
nikcleju@22:     xrec[algo[1]]   = np.zeros((d, y.shape[1]))
nikcleju@22:     err[algo[1]]    = np.zeros(y.shape[1])
nikcleju@22:     relerr[algo[1]] = np.zeros(y.shape[1])
nikcleju@22:   # Prepare storage variables for algorithms with Lambda    
nikcleju@22:   for i,algo in zip(np.arange(nalgosL),algosL):
nikcleju@22:     xrec[algo[1]]   = np.zeros((lambdas.size, d, y.shape[1]))
nikcleju@22:     err[algo[1]]    = np.zeros((lambdas.size, y.shape[1]))
nikcleju@22:     relerr[algo[1]] = np.zeros((lambdas.size, y.shape[1]))
nikcleju@19:   
nikcleju@22:   # Run algorithms non-Lambda
nikcleju@22:   for iy in np.arange(y.shape[1]):
nikcleju@22:     for algofunc,strname in algosN:
nikcleju@22:       epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
nikcleju@22:       xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega,epsilon)
nikcleju@22:       err[strname][iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
nikcleju@22:       relerr[strname][iy] = err[strname][iy] / np.linalg.norm(x0[:,iy])
nikcleju@22:   for algotuple in algosN:
nikcleju@22:     print algotuple[1],' : avg relative error = ',np.mean(relerr[strname])  
nikcleju@22: 
nikcleju@22:   # Run algorithms with Lambda
nikcleju@19:   for ilbd,lbd in zip(np.arange(lambdas.size),lambdas):
nikcleju@19:     for iy in np.arange(y.shape[1]):
nikcleju@22:       D = np.linalg.pinv(Omega)
nikcleju@22:       U,S,Vt = np.linalg.svd(D)
nikcleju@22:       for algofunc,strname in algosL:
nikcleju@19:         epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
nikcleju@22:         gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
nikcleju@22:         xrec[strname][ilbd,:,iy] = np.dot(D,gamma)
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:     print 'Lambda = ',lbd,' :'
nikcleju@22:     for algotuple in algosL:
nikcleju@22:       print '   ',algotuple[1],' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
nikcleju@10: 
nikcleju@22:   # Prepare results
nikcleju@22:   mrelerrN = dict()
nikcleju@22:   for algotuple in algosN:
nikcleju@22:     mrelerrN[algotuple[1]] = np.mean(relerr[algotuple[1]])
nikcleju@22:   mrelerrL = dict()
nikcleju@22:   for algotuple in algosL:
nikcleju@22:     mrelerrL[algotuple[1]] = np.mean(relerr[algotuple[1]],1)
nikcleju@22:   
nikcleju@22:   return mrelerrN,mrelerrL
nikcleju@29: 
nikcleju@29: def generateData(d,sigma,delta,rho,numvects,SNRdb):
nikcleju@29: 
nikcleju@29:   # Process parameters
nikcleju@29:   noiselevel = 1.0 / (10.0**(SNRdb/10.0));
nikcleju@29:   p = round(sigma*d);
nikcleju@29:   m = round(delta*d);
nikcleju@29:   l = round(d - rho*m);
nikcleju@29:   
nikcleju@29:   # Generate Omega and data based on parameters
nikcleju@29:   Omega = pyCSalgos.GAP.GAP.Generate_Analysis_Operator(d, p);
nikcleju@29:   # Optionally make Omega more coherent
nikcleju@29:   U,S,Vt = np.linalg.svd(Omega);
nikcleju@29:   Sdnew = S * (1+np.arange(S.size)) # Make D coherent, not Omega!
nikcleju@29:   Snew = np.vstack((np.diag(Sdnew), np.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
nikcleju@29:   Omega = np.dot(U , np.dot(Snew,Vt))
nikcleju@29: 
nikcleju@29:   # Generate data  
nikcleju@29:   x0,y,M,Lambda,realnoise = pyCSalgos.GAP.GAP.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
nikcleju@29:   
nikcleju@29:   return Omega,x0,y,M,realnoise
nikcleju@22:   
nikcleju@19: # Script main
nikcleju@19: if __name__ == "__main__":
nikcleju@27:   #import cProfile
nikcleju@27:   #cProfile.run('mainrun()', 'profile')    
nikcleju@29:   run()