nikcleju@10: # -*- coding: utf-8 -*-
nikcleju@10: """
nikcleju@17: Algorithms for approximate analysis recovery based on synthesis solvers (a.k.a. Analysis by Synthesis, ABS).
nikcleju@17: Approximate reconstruction, ABS-lambda.
nikcleju@10: 
nikcleju@17: Author: Nicolae Cleju
nikcleju@10: """
nikcleju@17: __author__ = "Nicolae Cleju"
nikcleju@17: __license__ = "GPL"
nikcleju@17: __email__ = "nikcleju@gmail.com"
nikcleju@17: 
nikcleju@10: 
nikcleju@10: import numpy
nikcleju@17: 
nikcleju@17: # Import synthesis solvers from pyCSalgos package
nikcleju@10: import pyCSalgos.BP.l1qc
nikcleju@13: import pyCSalgos.SL0.SL0_approx
nikcleju@10: import pyCSalgos.OMP.omp_QR
nikcleju@10: import pyCSalgos.TST.RecommendedTST
nikcleju@10: 
nikcleju@10: def sl0(y,M,Omega,epsilon,lbd,sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, A_pinv=None, true_s=None):
nikcleju@17:   """
nikcleju@17:   ABS-lambda: Smooth L0 (SL0)
nikcleju@17:   """    
nikcleju@10:   N,n = Omega.shape
nikcleju@10:   D = numpy.linalg.pinv(Omega)
nikcleju@10:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@10:   aggDupper = numpy.dot(M,D)
nikcleju@10:   aggDlower = Vt[-(N-n):,:]
nikcleju@13:   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
nikcleju@10:   aggy = numpy.concatenate((y, numpy.zeros(N-n)))
nikcleju@10:   
nikcleju@21:   #return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s))
nikcleju@21:   return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx_dai(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s))
nikcleju@10:   
nikcleju@13: def bp(y,M,Omega,epsilon,lbd, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False):
nikcleju@17:   """
nikcleju@17:   ABS-lambda: Basis Pursuit (based on l1magic toolbox)
nikcleju@17:   """    
nikcleju@10:   N,n = Omega.shape
nikcleju@10:   D = numpy.linalg.pinv(Omega)
nikcleju@10:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@10:   aggDupper = numpy.dot(M,D)
nikcleju@10:   aggDlower = Vt[-(N-n):,:]
nikcleju@13:   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
nikcleju@10:   aggy = numpy.concatenate((y, numpy.zeros(N-n)))
nikcleju@10: 
nikcleju@16:   return numpy.dot(D, pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon, lbtol, mu, cgtol, cgmaxiter, verbose))
nikcleju@10: 
nikcleju@10: def ompeps(y,M,Omega,epsilon,lbd):
nikcleju@17:   """
nikcleju@17:   ABS-lambda: OMP with stopping criterion residual < epsilon 
nikcleju@17:   """
nikcleju@10:   N,n = Omega.shape
nikcleju@10:   D = numpy.linalg.pinv(Omega)
nikcleju@10:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@10:   aggDupper = numpy.dot(M,D)
nikcleju@10:   aggDlower = Vt[-(N-n):,:]
nikcleju@21:   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
nikcleju@10:   aggy = numpy.concatenate((y, numpy.zeros(N-n)))
nikcleju@10:   
nikcleju@10:   opts = dict()
nikcleju@10:   opts['stopCrit'] = 'mse'
nikcleju@10:   opts['stopTol'] = epsilon**2 / aggy.size
nikcleju@16:   return numpy.dot(D, pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0])
nikcleju@10:   
nikcleju@10: def tst_recom(y,M,Omega,epsilon,lbd, nsweep=300, xinitial=None, ro=None):
nikcleju@17:   """
nikcleju@17:   ABS-lambda: Two Stage Thresholding (TST) with optimized parameters (see Maleki & Donoho)
nikcleju@17:   """  
nikcleju@10:   N,n = Omega.shape
nikcleju@10:   D = numpy.linalg.pinv(Omega)
nikcleju@10:   U,S,Vt = numpy.linalg.svd(D)
nikcleju@10:   aggDupper = numpy.dot(M,D)
nikcleju@10:   aggDlower = Vt[-(N-n):,:]
nikcleju@13:   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
nikcleju@10:   aggy = numpy.concatenate((y, numpy.zeros(N-n)))
nikcleju@10:   
nikcleju@10:   tol = epsilon / numpy.linalg.norm(aggy)
nikcleju@21:   return numpy.dot(D, pyCSalgos.TST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep, tol, xinitial, ro))
nikcleju@10: