view algos.py @ 13:a2d881253324

In working, not debugged yet

# -*- coding: utf-8 -*-
"""
Define simple wrappers for algorithms

Created on Wed Dec 07 14:06:13 2011

@author: ncleju
"""

import numpy
import pyCSalgos
import pyCSalgos.GAP.GAP
import ABSexact
import ABSmixed
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
###---------------------------------
def run_exact_gap(y,M,Omega):
  """
  Wrapper for GAP algorithm for exact analysis recovery
  """
  gapparams = {"num_iteration" : 1000,\
                   "greedy_level" : 0.9,\
                   "stopping_coefficient_size" : 1e-4,\
                   "l2solver" : 'pseudoinverse',\
                   "noise_level": 1e-10}
  return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,numpy.zeros(Omega.shape[1]))[0]
    
def run_exact_bp(y,M,Omega):
  """
  Wrapper for BP algorithm for exact analysis recovery
  Algorithm implementation is l1eq_pd() from l1-magic toolbox
  """
  return ABSexact.bp(y,M,Omega,numpy.zeros(Omega.shape[0]))

def run_exact_ompeps(y,M,Omega):
  """
  Wrapper for OMP algorithm for exact analysis recovery, with stopping criterion = epsilon
  """
  return ABSexact.ompeps(y,M,Omega,1e-9)
  
#def run_exact_ompk(y,M,Omega)
#  """
#  Wrapper for OMP algorithm for exact analysis recovery, with stopping criterion = fixed no. of atoms
#  """

def run_exact_sl0(y,M,Omega):
  """
  Wrapper for SL0 algorithm for exact analysis recovery
  """
  sigma_min = 1e-12
  sigma_decrease_factor = 0.5
  mu_0 = 2
  L = 20
  return ABSexact.sl0(y,M,Omega, sigma_min, sigma_decrease_factor, mu_0, L)

def run_exact_tst(y,M,Omega):
  """
  Wrapper for TST algorithm (with default optimized params) for exact analysis recovery
  """
  nsweep = 300
  tol = 1e-5
  return ABSexact.tst_recom(y,M,Omega, nsweep, tol)


###---------------------------------------
### Approximate reconstruction algorithms
###---------------------------------------
#  1. Native

def run_gap(y,M,Omega,epsilon):
  """
  Wrapper for GAP algorithm for approximate analysis recovery
  """
  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_nesta(y,M,Omega,epsilon):
  """
  Wrapper for NESTA algorithm for approximate analysis recovery
  """  
  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]

#  2. ABS-mixed

def run_mixed_sl0(y,M,Omega,epsilon):
  """
  Wrapper for SL0-mixed algorithm for approximate analysis recovery
  """ 
  sigma_min = 0.001
  sigma_decrease_factor = 0.5
  mu_0 = 2
  L = 10
  return ABSmixed.sl0(y,M,Omega,epsilon,sigma_min, sigma_decrease_factor, mu_0, L)

def run_mixed_bp(y,M,Omega,epsilon):
  """
  Wrapper for BP-mixed algorithm for approximate analysis recovery
  """   
  return ABSmixed.bp(y,M,Omega,epsilon, numpy.zeros(Omega.shape[0]))

#  3. ABS-lambda

def run_lambda_sl0(y,M,Omega,epsilon,lbd):
  """
  Wrapper for SL0 algorithm within ABS-lambda approach for approximate analysis recovery
  """    
  sigma_min = 0.001
  sigma_decrease_factor = 0.5
  mu_0 = 2
  L = 10
  return ABSlambda.sl0(y,M,Omega,epsilon, lbd, sigma_min, sigma_decrease_factor, mu_0, L)

def run_lambda_bp(y,M,Omega,epsilon,lbd):
  """
  Wrapper for BP algorithm within ABS-lambda approach for approximate analysis recovery
  """     
  return ABSlambda.bp(y,M,Omega,epsilon,lbd,numpy.zeros(Omega.shape[0]))

def run_lambda_ompeps(y,M,Omega,epsilon,lbd):
  """
  Wrapper for OMP algorithm, with stopping criterion = epsilon,
  for approximate analysis recovery within ABS-lambda approach
  """     
  return ABSlambda.ompeps(y,M,Omega,epsilon,lbd)

def run_lambda_tst(y,M,Omega,epsilon,lbd):
  """
  Wrapper for TST algorithm (with default optimized params)
  for approximate analysis recovery within ABS-lambda approach
  """
  nsweep = 300
  return ABSlambda.tst_recom(y,M,Omega,epsilon,lbd, nsweep)
  
  
### Algorithm tuples
## Exact recovery
exact_gap = (run_exact_gap, 'GAP')
exact_bp = (run_exact_bp, 'ABS-exact: BP')
exact_ompeps = (run_exact_ompeps, 'ABS-exact: OMP-eps')
exact_sl0 = (run_exact_sl0, 'ABS-exact: SL0')
exact_tst = (run_exact_tst, 'ABS-exact: TST')
## Approximate recovery
# Native
gap = (run_gap, 'GAP')
nesta = (run_nesta, 'NESTA')
# ABS-mixed
mixed_sl0 = (run_mixed_sl0, 'ABS-mixed: SL0')
mixed_bp = (run_mixed_bp, 'ABS-mixed: BP')
# ABS-lambda
lambda_sl0 = (run_lambda_sl0, 'ABS-lambda: SL0')
lambda_bp = (run_lambda_bp, 'ABS-lambda: BP')
lambda_ompeps = (run_lambda_ompeps, 'ABS-lambda: OMP-eps')
lambda_tst = (run_lambda_tst, 'ABS-lambda: 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')