% COMPARE_FXN 	Compares real-valued FBMP with other popular sparse reconstruction
%               algorithms, collecting NMSE, runtime, and sparsity data
%
% SYNTAX:   [NMSEs Times Sparse] = compare_fxn(y, A, sig2s, p1, mus, ...
%                                              sig2w, D, x)
%
% INPUTS:
%           y:  An M-by-1 column vector of observations
%           A: An M-by-N measurement matrix with unit-norm columns, (N>>M)
%           sig2s:  2-by-1 vector of sparse coefficient variances,  ([off; on]).  
%                   Note that the variances of all Q-1 active Gaussian mixture
%                   densities will be equal
%           p1:	Prior probability of active taps (i.e. Pr{s_n = q} = p1/(Q-1),
%               q > 0)
%           mus:    Q-by-1 vector of means for the Q-ary Gaussian mixture densities,
%                   ([off; on_1; on_2; ...; on_Q-1])
%           sig2w:  Noise variance scalar
%           D:  Maximum number of RGS iterations
%           x:  An N-by-1 vector of sparse coefficients for MSE comparison
%
% OUTPUTS:
%           NMSEs: A structure containing the following elements:
%                   .NMSE_fbmp: NMSE btwn x and xmmse
%                   .NMSE_bayes: NMSE btwn x and X_bayes
%                   .NMSE_omp: NMSE btwn x and X_omp
%                   .NMSE_stomp: NMSE btwn x and X_stomp
%                   .NMSE_gpsr: NMSE btwn x and X_gpsr (best of 2 options)
%                   .NMSE_bcs: NMSE btwn x and X_bcs
%                   .NMSE_vb_bcs: NMSE btwn x and X_vb_bcs
%           Times: A structure containing the following elements:
%                   .fbmp_time: Runtime for FBMP algorithm
%                   .bayes_time: Runtime for SparseBayes algorithm
%                   .omp_time: Runtime for OMP algorithm
%                   .stomp_time: Runtime for StOMP algorithm
%                   .gpsr_time: Runtime for GPSR algorithm
%                   .bcs_time: Runtime for BCS algorithm
%                   .vb_bcs_time: Runtime for VB-BCS algorithm
%           Sparse: A structure containing the following elements:
%                   .fbmp_count: Active-tap count for FBMP algorithm
%                   .bayes_count: Active-tap count for SparseBayes algorithm
%                   .omp_count: Active-tap count for OMP algorithm
%                   .stomp_count: Active-tap count for StOMP algorithm
%                   .gpsr_count: Active-tap count for GPSR algorithm
%                   .bcs_count: Active-tap count for BCS algorithm
%                   .vb_bcs_count: Active-tap count for VB-BCS algorithm
%

%
% Coded by: Justin Ziniel, Ohio State Univ.
% Email: ziniel.1@osu.edu
% Last modified: April 2, 2009
% Change summary: Added NMSE of FBMP's MAP estimator, & sparsity (4/2/09)
%                 Added VB-BCS to the comparison (4/2/09)
%                 Adjusted parameters for most algorithms (4/2/09)
% FBMP version 1.2 
% Copyright (c) Philip Schniter, Lee C. Potter, and Justin Ziniel, 2008
%

function [NMSEs Times Sparse] = compare_fxn(y, A, sig2s, p1, mus, sig2w, ...
                                            D, x)
                                     

%% Variables needed by some algorithms

[M, N] = size(A);


%% FBMP-RGS Recovery
% Algorithm: "Fast Bayesian Matching Pursuit", 
% Philip Schniter, Lee C. Potter, and Justin Ziniel

stop = 0;   % Stop at zero std. devs. below E[nu|s,y]

try
    t = cputime;

    [xmmse, xmmse_star] = fbmpr_fxn_reduced(y, A, p1, sig2w, sig2s, mus, ...
                                           D, stop);

    Times.fbmp_time = cputime - t;
    NMSEs.NMSE_fbmp(1) = norm(x - xmmse)^2/(norm(x)^2); % NMSE-MMSE
    NMSEs.NMSE_fbmp(2) = norm(x - xmmse_star{1})^2/(norm(x)^2); % NMSE-MAP
    Sparse.fbmp_count(1) = sum(xmmse ~= 0);      % MMSE sparsity
    Sparse.fbmp_count(2) = sum(xmmse_star{1} ~= 0);    % MAP sparsity
catch
    Times.fbmp_refine_time = NaN;
    NMSEs.NMSE_fbmp_refine = [NaN, NaN];
    Sparse.fbmp_count = [NaN, NaN];
end


%%  SparseBayes Recovery
% Algorithm: "SparseBayes", 
% Michael E. Tipping

setEnvironment('InfoLevel', 0);     % Least verbose output

alpha = 1;      % Initial hyperparameter (suggested default)
beta = -1/sig2w;  % Noise variance (fixed)
maxiter = N;    % Maximum number of iterations

% Run the Sparse Bayesian Learning Algorithm
try
    t = cputime;
    [weights, index] = sbl_estimate(A, y, alpha, beta, maxiter);
    Times.bayes_time = cputime - t;

    % Reconstruct the output signal
    X_bayes = zeros(N,1);
    X_bayes(index) = weights;

    NMSEs.NMSE_bayes = norm(x - X_bayes)^2/(norm(x)^2);
    Sparse.bayes_count = sum(X_bayes ~= 0);
catch
    Times.bayes_time = NaN;
    NMSEs.NMSE_bayes = NaN;
    Sparse.bayes_count = NaN;
end
    

%% OMP Recovery
% Algorithm: "Orthogonal Matching Pursuit", 
% Stanford's SparseLab software package

stop_tol = sqrt(2)*sqrt(sig2w)*gamma((M+1)/2)/gamma(M/2);
stop_tol = stop_tol/norm(y);    % Modification to impose the tolerance
                                % to the value computed on the previous
                                % line (OMP multiplies threshold by
                                % norm(y))

try
    t = cputime;
    [X_omp] = SolveOMP(A, y, N, M, 0, 0, 0, stop_tol);
    Times.omp_time = cputime - t;

    NMSEs.NMSE_omp = norm(x - X_omp)^2/(norm(x)^2);
    Sparse.omp_count = sum(X_omp ~= 0);
catch
    Times.omp_time = NaN;
    NMSEs.NMSE_omp = NaN;
    Sparse.omp_count = NaN;
end


%% StOMP Recovery
% Algorithm: "Stagewise Orthogonal Matching Pursuit", 
% Stanford's SparseLab software package

iters = 10;     % # of iteration stages (default)

% % Threshold parameter (same formula as the one in examples used to test StOMP)
thresh = (M/N)*(1 - (N - max(find(cumsum(sort(x.^2, 'ascend')) <= norm(y - A*x)^2)))/M)/iters;

stop_tol = sqrt(2)*sqrt(sig2w)*gamma((M+1)/2)/gamma(M/2);
stop_tol = stop_tol/norm(y);    % Modification to impose the tolerance
                                % to the value computed on the previous
                                % line (StOMP multiplies threshold by
                                % norm(y))
% q = 0.35;       % False discovery rate (empirically seems to work well)

try
    t = cputime;
    [X_stomp] = SolveStOMP(A, y, N, 'FAR', 0.5*thresh, iters, 0, stop_tol);
%     [X_stomp] = SolveStOMP(A, y, N, 'FDR', q, iters, 0, stop_tol);
    Times.stomp_time = cputime - t;

    NMSEs.NMSE_stomp = norm(x - X_stomp)^2/(norm(x)^2);
    Sparse.stomp_count = sum(X_stomp ~= 0);
catch
    Times.stomp_time = NaN;
    NMSEs.NMSE_stomp = NaN;
    Sparse.stomp_count = NaN;
end


%% GPSR Recovery
% Algorithm: "Gradient Projection for Sparse Reconstruction", 
% Mario Figueiredo, Robert D. Nowak, and Stephen J. Wright

tau = 0.1*max(abs(A'*y));       % See ref [35] in the GPSR paper for justification

try
    t = cputime;
    [X_gpsr, X_gpsr_db] = GPSR_Basic(y, A, tau, 'Debias', 1, 'Verbose', 0);
    Times.gpsr_time = cputime - t;

    % Calculate the NMSE using the better of two recoveries (debiased and
    % un-debiased)
    [NMSEs.NMSE_gpsr] = min([norm(x - X_gpsr)^2/(norm(x)^2), ...
                                norm(x - X_gpsr_db)^2/(norm(x)^2)]);
    Sparse.gpsr_count = sum(X_gpsr ~= 0);
catch
    Times.gpsr_time = NaN;
    [NMSEs.NMSE_gpsr] = NaN;
    Sparse.gpsr_count = NaN;
end


%% Bayesian Compressive Sensing (BCS) Recovery
% Algorithm: "Bayesian Compressive Sensing", 
% Shihao Ji, Ya Xue, and Lawrence Carin

eta = 1e-4;     % Stopping threshold

try
    t = cputime;
    [weights, index] = BCS_fast_rvm(A, y, sig2w, eta);
    Times.bcs_time = cputime - t;

    % Reconstruct the output signal
    X_bcs = zeros(N,1);
    X_bcs(index) = weights;
    
    NMSEs.NMSE_bcs = norm(x - X_bcs)^2/(norm(x)^2);
    Sparse.bcs_count = sum(X_bcs ~= 0);
catch
    Times.bcs_time = NaN;
    NMSEs.NMSE_bcs = NaN;
    Sparse.bcs_count = NaN;
end


%% Variational Bayesian - Bayesian Compressive Sensing (VB-BCS) Recovery
% Algorithm: "Variational Bayesian - Bayesian Compressive Sensing", 
% Shihao Ji, Ya Xue, and Lawrence Carin

try
    % Use default suggested values for most optional parameters
    hyperpara.a0 = 1e-6;    % Noise precision 1
    hyperpara.b0 = 1e-6;    % Noise precision 2
    hyperpara.c0 = 1e-6;    % Non-zero coefficients precision 1
    hyperpara.d0 = 1e-6;    % Non-zero coefficients precision 2
	niter = 40;             % Maximum number of VB iterations
    tol = 1e-3;             % Tolerance of relative change of VB lower 
                            % bound for two consecutive steps
    plotflag = 0;           % Non-verbose
    
    t = cputime;
    [X_vb_bcs, posterior, bound] = BCSvb(A, y, hyperpara, niter, tol, ...
                                         plotflag);
    Times.vb_bcs_time = cputime - t;

    NMSEs.NMSE_vb_bcs = norm(x - X_vb_bcs)^2/(norm(x)^2);
    Sparse.vb_bcs_count = sum(X_vb_bcs ~= 0);
catch
    Times.vb_bcs_time = NaN;
    NMSEs.NMSE_vb_bcs = NaN;
    Sparse.vb_bcs_count = NaN;
end
