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