
clear all; 

% load the dictionary matrix
load Dic;
Dtn = Dic;
% problem dimension
[N,M] = size(Dtn);
D = 3;              % diversity (number of active sources)
SNR = 25;           % SNR
iterNum = 1000;      % number of experiments

 

for it = 1 : iterNum
  
    % dictionary matrix with columns draw uniformly from the surface of a unit hypersphere
    Phi = Dtn./(ones(N,1)*sqrt(sum(Dtn.^2)));

    % nonzero Rows
    nonzeroW = sign(randn(D,1)).* ( rand(D,1)*0.5 + 0.5 );
     

    % select active sources at random locations
    ind = randperm(M);
    indice = ind(1:D);
    Wgen = zeros(M,1);
    Wgen(indice,:) = nonzeroW;

    % noiseless signal
    signal = Phi * Wgen;

    % observation noise
    stdnoise = std(signal)*10^(-SNR/20);
    noise = randn(N,1).*(ones(N,1)*stdnoise);                   
    

    % noisy signal
    T = signal + noise;
    
    fprintf('\nTrial: %d\n',it);

    %============================ TMSBL ==========================
    X_tsbl = TMSBL(Phi, T, 'noise','mild');

    % Failure Rate
    [F1,P1,R1,ind1] = perfSupp(X_tsbl,indice,'firstlargest', D);
    fail_TMSBL(it) = (F1~=1);

    % MSE
    mse_TMSBL(it) = (norm(Wgen - X_tsbl,'fro')/norm(Wgen,'fro'))^2;

    fprintf(' TMSBL : Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_TMSBL(it),mse_TMSBL(it),mean(fail_TMSBL),mean(mse_TMSBL));

    %============================== ExCoV ========================
    X_excov = ExCoV(Phi,T,'Visibility',0);

    % Failure Rate
    [F2] = perfSupp(X_excov,indice,'firstlargest', D);
    fail_EXCOV(it) = (F2~=1);

    % MSE
    mse_EXCOV(it) = (norm(Wgen - X_excov,'fro')/norm(Wgen,'fro'))^2;

    fprintf(' ExCov : Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_EXCOV(it),mse_EXCOV(it),mean(fail_EXCOV),mean(mse_EXCOV));
    
    
    %============================== CoSaMP ========================
    X_cosamp = cosamp(Phi,T,D,1e-5);
    
    % Failure Rate
    F3 = perfSupp(X_cosamp,indice,'firstlargest', D);
    fail_cosamp(it) = (F3~=1);

    % MSE
    mse_cosamp(it) = (norm(Wgen - X_cosamp,'fro')/norm(Wgen,'fro'))^2;

    fprintf(' CoSaMP : Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_cosamp(it),mse_cosamp(it),mean(fail_cosamp),mean(mse_cosamp));
    
    %============================== Subspace Pursuit ========================
    Xsub = CSRec_SP(D,Phi,T);
    
    % Failure Rate
    F4 = perfSupp(Xsub.x_hat,indice,'firstlargest', D);
    fail_subspace(it) = (F4~=1);

    % MSE
    mse_subspace(it) = (norm(Wgen - Xsub.x_hat,'fro')/norm(Wgen,'fro'))^2;

    fprintf('Subspace Pursuit : Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_subspace(it),mse_subspace(it),mean(fail_subspace),mean(mse_subspace));
    
    
    %============================== Approximate Message Passing ========================
    Xamp = reconstructAmp(Phi, T, 500, 1e-3);
    
    % Failure Rate
    F5 = perfSupp(Xamp,indice,'firstlargest', D);
    fail_amp(it) = (F5~=1);

    % MSE
    mse_amp(it) = (norm(Wgen - Xamp,'fro')/norm(Wgen,'fro'))^2;

    fprintf('Appro.Msg.Pass.: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_amp(it),mse_amp(it),mean(fail_amp),mean(mse_amp));
    
    %============================== Bayesian Compressive Sensing ========================
    initsigma2 = std(T)^2/1e2;
    [X_coef,used6] = BCS_fast_rvm(Phi,T,initsigma2,1e-8);
    
    X_bcs = zeros(M,1);
    X_bcs(used6) = X_coef;

    % Failure Rate
    F6 = perfSupp(X_bcs,indice,'firstlargest', D);
    fail_bcs(it) = (F6~=1);

    % MSE
    mse_bcs(it) = (norm(Wgen - X_bcs,'fro')/norm(Wgen,'fro'))^2;

    fprintf('Bayesian CS: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_bcs(it),mse_bcs(it),mean(fail_bcs),mean(mse_bcs));
    
    %============================== EM-SBL ========================
    ini_lambda = std(T)^2/1e2;
    X_emsbl = EMSBL(Phi, T, ini_lambda, 1);
     
    % Failure Rate
    F7 = perfSupp(X_emsbl,indice,'firstlargest', D);
    fail_emsbl(it) = (F7~=1);

    % MSE
    mse_emsbl(it) = (norm(Wgen - X_emsbl,'fro')/norm(Wgen,'fro'))^2;

    fprintf('  EM-SBL: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_emsbl(it),mse_emsbl(it),mean(fail_emsbl),mean(mse_emsbl));
    
    %============================== Magic L1 ========================  
    x0 = Phi'*T;
    epsilon =  stdnoise*sqrt(D)*sqrt(1 + 2*sqrt(2)/sqrt(D));                 
    X_l1 = l1qc_logbarrier(x0, Phi, [], T, epsilon, 1e-3);
 
    % Failure Rate
    F8 = perfSupp(X_l1,indice,'firstlargest', D);
    fail_l1(it) = (F8~=1);

    % MSE
    mse_l1(it) = (norm(Wgen - X_l1,'fro')/norm(Wgen,'fro'))^2;

    fprintf('Magic L1: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_l1(it),mse_l1(it),mean(fail_l1),mean(mse_l1));

    %============================== Hard Thresholding Pursuit algorithm ========================  
    % implemented as a noiseless case (have removed the noise)
    [X_htp] = HTP(signal,Phi,D); 
 
    % Failure Rate
    F10 = perfSupp(X_htp,indice,'firstlargest', D);
    fail_htp(it) = (F10~=1);

    % MSE
    mse_htp(it) = (norm(Wgen - X_htp,'fro')/norm(Wgen,'fro'))^2;

    fprintf('HTP: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_htp(it),mse_htp(it),mean(fail_htp),mean(mse_htp));

 
    %========================= Fast Bayesian Matching Pursuit (know SNR) ========================  
    
    p1 = D/M;
    sig2w = norm(Phi*Wgen)^2/N*10^(-SNR/10);
    sig2s = [0; 1];
    mus = [0; 0];
    LP = 30;
    stop = 0;
    [X_fbmp, xmmse_star, psy_star, nu_star, T_star, d_tot, d_max] = ...
    fbmpr_fxn_reduced(T, Phi, p1, sig2w, sig2s, mus, LP, stop); 
 
    % Failure Rate
    F11 = perfSupp(X_fbmp,indice,'firstlargest', D);
    fail_fbmp(it) = (F11~=1);

    % MSE
    mse_fbmp(it) = (norm(Wgen - X_fbmp,'fro')/norm(Wgen,'fro'))^2;

    fprintf('FBMP: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_fbmp(it),mse_fbmp(it),mean(fail_fbmp),mean(mse_fbmp));

    
    %========================= FOCUSS (know SNR) ========================  
    
    sig2 = norm(Phi*Wgen)^2/N*10^(-SNR/10);
    [X_focuss] = MFOCUSS(Phi, T, sig2, 'p',0.8);
 
    % Failure Rate
    F12 = perfSupp(X_focuss,indice,'firstlargest', D);
    fail_focuss(it) = (F12~=1);

    % MSE
    mse_focuss(it) = (norm(Wgen - X_focuss,'fro')/norm(Wgen,'fro'))^2;

    fprintf('FOCUSS: Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_focuss(it),mse_focuss(it),mean(fail_focuss),mean(mse_focuss));
    
    
    % ========================= SL0 (removed noise) ==============
    sigma_min = 0.001;
    sigma_decrease_factor = 0.5;
    X_sl0 = SL0(Phi, signal, sigma_min, sigma_decrease_factor);
    % Failure Rate
    F13 = perfSupp(X_sl0,indice,'firstlargest', D);
    fail_sl0(it) = (F13~=1);

    % MSE
    mse_sl0(it) = (norm(Wgen - X_sl0,'fro')/norm(Wgen,'fro'))^2;

    fprintf(' SL0 : Failure Rate = %4.3f; MSE = %5.4f; Ave-Failure Rate = %5.4f; Ave-MSE = %5.4f\n',...
        fail_sl0(it),mse_sl0(it),mean(fail_sl0),mean(mse_sl0));
    
    
    
    save result;
end
