nikcleju@67: # -*- coding: utf-8 -*-
nikcleju@67: """
nikcleju@67: Algorithms for exact analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS).
nikcleju@67: Exact reconstruction.
nikcleju@67: 
nikcleju@67: Author: Nicolae Cleju
nikcleju@67: """
nikcleju@67: __author__ = "Nicolae Cleju"
nikcleju@67: __license__ = "GPL"
nikcleju@67: __email__ = "nikcleju@gmail.com"
nikcleju@67: 
nikcleju@67: 
nikcleju@67: import numpy
nikcleju@67: 
nikcleju@67: # Import synthesis solvers from pyCSalgos package
nikcleju@67: import pyCSalgos.BP.l1eq_pd
nikcleju@67: import pyCSalgos.BP.cvxopt_lp
nikcleju@67: import pyCSalgos.OMP.omp_QR
nikcleju@67: import pyCSalgos.SL0.SL0
nikcleju@67: import pyCSalgos.TST.RecommendedTST
nikcleju@67: 
nikcleju@67: 
nikcleju@67: def bp(y,M,Omega,x0, pdtol=1e-3, pdmaxiter=50, cgtol=1e-8, cgmaxiter=200, verbose=False):
nikcleju@67:   """
nikcleju@67:   ABS-exact: Basis Pursuit (based on l1magic toolbox)
nikcleju@67:   """
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67: 
nikcleju@67:   return numpy.dot(D , pyCSalgos.BP.l1eq_pd.l1eq_pd(x0,Atilde,Atilde.T,ytilde, pdtol, pdmaxiter, cgtol, cgmaxiter, verbose))
nikcleju@67: 
nikcleju@67: def bp_cvxopt(y,M,Omega):
nikcleju@67:   """
nikcleju@67:   ABS-exact: Basis Pursuit (based cvxopt)
nikcleju@67:   """
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67: 
nikcleju@67:   return numpy.dot(D , pyCSalgos.BP.cvxopt_lp.cvxopt_lp(ytilde, Atilde))
nikcleju@67: 
nikcleju@67: 
nikcleju@67: def ompeps(y,M,Omega,epsilon):
nikcleju@67:   """
nikcleju@67:   ABS-exact: OMP with stopping criterion residual < epsilon 
nikcleju@67:   """
nikcleju@67:   
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67:   
nikcleju@67:   opts = dict()
nikcleju@67:   opts['stopCrit'] = 'mse'
nikcleju@67:   opts['stopTol'] = epsilon
nikcleju@67:   return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0])
nikcleju@67: 
nikcleju@67: def ompk(y,M,Omega,k):
nikcleju@67:   """
nikcleju@67:   ABS-exact: OMP with stopping criterion fixed number of atoms = k
nikcleju@67:   """
nikcleju@67:   
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67:   
nikcleju@67:   opts = dict()
nikcleju@67:   opts['stopTol'] = k
nikcleju@67:   return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0])
nikcleju@67: 
nikcleju@67: def sl0(y,M,Omega, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, true_s=None):
nikcleju@67:   """
nikcleju@67:   ABS-exact: Smooth L0 (SL0)
nikcleju@67:   """
nikcleju@67:   
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67:   
nikcleju@67:   return numpy.dot(D, pyCSalgos.SL0.SL0.SL0(Atilde,ytilde,sigma_min,sigma_decrease_factor,mu_0,L,true_s))
nikcleju@67: 
nikcleju@67: def tst_recom(y,M,Omega, nsweep=300, tol=0.00001, xinitial=None, ro=None):
nikcleju@67:   """
nikcleju@67:   ABS-exact: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho)
nikcleju@67:   """
nikcleju@67:   
nikcleju@67:   N,n = Omega.shape
nikcleju@67:   D = numpy.linalg.pinv(Omega)
nikcleju@67:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@67:   Aextra = Vt[-(N-n):,:]
nikcleju@67:   
nikcleju@67:   # Create aggregate problem
nikcleju@67:   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
nikcleju@67:   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
nikcleju@67:   
nikcleju@67:   return numpy.dot(D, pyCSalgos.TST.RecommendedTST.RecommendedTST(Atilde, ytilde, nsweep, tol, xinitial, ro))
nikcleju@67: