%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Compressive Sensing Reconstruction Example (1-d example)    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Coded by Kun Qiu (kqiu@iastate.edu)
%Last updated (Dec. 13, 2009)


clear all
clc
close all

path(path,'./Algorithms');

%Parameter specification
SignalType=0;     %Binary signal=0;Gaussian signal=1;
m=512;                %number of underlying signal length
T=20;                   %number of spikes
N=100;                 %number of measurements
sigma2=1e-5;      %noise variance
TrialNum=1;          %number of trials
visibility=isequal(TrialNum,1);

RunOtherMethod=1;           %Option to run other methods in addition to ExCoV (0=not run; 1=run)

for ll=1:TrialNum
    %Signal and measurement model
    Psi=eye(m);           %the sparsity basis
    Phi=randn(N,m);
    Phi=orth(Phi')';
    H=Phi*Psi;              %effective sensing matrix

    %Sparse signal generation
    s=zeros(m,1);
    q=randperm(m);
    if SignalType==0
        s(q(1:T))=sign(randn(T,1));       %Binary signal
        SignalType_Str='Binary';
    else
        s(q(1:T))=randn(T,1);
        SignalType_Str='Gaussian';
    end
    noi=sqrt(sigma2)*randn(N,1);      %noise

    %Take measurements
    y=Phi*s+noi;
    
    %solve by ExCoV
    tic;
    [x_ExCoV,sigma2_ExCoV,A_index_ExCoV,Count_ExCoV(ll)]=ExCoV(H,y,'Visibility',visibility); 
    t_ExCoV(ll)=toc;
    error_ExCoV(ll)=norm(x_ExCoV-s)^2/m;
    
    
    if RunOtherMethod
        
        %solve by ExCoV
        tic;
        [x_ExCoVapp,sigma2_ExCoVapp,A_index_ExCoVapp,Count_ExCoVapp(ll)]=ExCoVapp(H,y,'Visibility',0);
        t_ExCoVapp(ll)=toc;
        error_ExCoVapp(ll)=norm(x_ExCoVapp-s)^2/m;
        
        %solve by NIHT
        tic;
        [x_NIHT,Count_NIHT(ll)]=hard_l0_Mterm_mod(y,H,m,T,1e-9);
        t_NIHT(ll)=toc;
        error_NIHT(ll)=norm(x_NIHT-s)^2/m;

        %solve by IHT
        tic;
        [x_IHT,Count_IHT(ll)]=hard_l0_Mterm_mod(y,H,m,T,1e-9,'step_size',1);
        t_IHT(ll)=toc;
        error_IHT(ll)=norm(x_IHT-s)^2/m;

        %solve by SBL
        tic;
        [x_SBL,Count_SBL(ll)]=SBL(H,y,1e-9);
        t_SBL(ll)=toc;
        error_SBL(ll)=norm(x_SBL-s)^2/m;
        
        %solve by GPSR
        tic;
        tau=0.01*max(abs(H'*y));
        [x_GPSR,x_GPSRdb]=GPSR_BB(y,H,tau,'Debias',1,'Verbose',0,'ToleranceA',1e-5);
        t_GPSR(ll)=toc;
        if isempty(x_GPSRdb)
            x_GPSRdb=x_GPSR;
        end
        error_GPSR(ll)=norm(x_GPSR-s)^2/m;
        error_GPSRdb(ll)=norm(x_GPSRdb-s)^2/m;

        % solve by SOCP
        x0=Phi'*inv(Phi*Phi')*y;
        epsilon=sqrt(sigma2)*sqrt(N)*sqrt(1+2*sqrt(2)/sqrt(N));
        tic;
        x_SOCP=l1qc_logbarrier(x0,H,[],y,epsilon,1e-3);
        t_SOCP(ll)=toc;
        fprintf(1,'SOCP number of nonzero weights: %d\n',sum(x_SOCP~=0));
        error_SOCP(ll)=norm(x_SOCP-s)^2/m;

        % solve by CoSaMP
        tic;
        [x_CoSaMP,A_index_CoSaMP]=CoSaMP(H,y,T);
        t_CoSaMP(ll)=toc;
        x_CoSaMPdb=zeros(m,1);
        HA=H(:,A_index_CoSaMP);
        x_CoSaMPdb(A_index_CoSaMP)=HA\y;
        error_CoSaMP(ll)=norm(x_CoSaMP-s)^2/m;
        error_CoSaMPdb(ll)=norm(x_CoSaMPdb-s)^2/m;

        % solve by BCS
        initsigma2=std(y)^2/1e2;
        tic;
        [weights,used,sigma2_BCS,errbars]=BCS_fast_rvm(H,y,initsigma2,1e-8);
        t_BCS(ll)=toc;
        x_BCS=zeros(m,1); 
        err = zeros(m,1);
        x_BCS(used)=weights;
        err(used)=errbars;    
        fprintf(1,'BCS number of nonzero weights: %d\n',length(used));
        error_BCS(ll)=norm(x_BCS-s)^2/m;
    end
    
    %solve by clairvoyant estimate
    x_LS=zeros(m,1);
    HA=H(:,q(1:T));
    x_LS(q(1:T))=HA\y;
    error_LS(ll)=norm(x_LS-s)^2/m;
    
    clc
    display(['Trial ',num2str(ll),' completed (total ',num2str(TrialNum),' trials)']);
end

MSE_ExCoV=sum(error_ExCoV)/TrialNum;
TimeAvg_ExCoV=sum(t_ExCoV)/TrialNum;
CountAvg_ExCoV=sum(Count_ExCoV)/TrialNum;

if RunOtherMethod
    MSE_ExCoVapp=sum(error_ExCoVapp)/TrialNum;
    TimeAvg_ExCoVapp=sum(t_ExCoVapp)/TrialNum;
    CountAvg_ExCoVapp=sum(Count_ExCoVapp)/TrialNum;
    
    MSE_NIHT=sum(error_NIHT)/TrialNum;
    TimeAvg_NIHT=sum(t_NIHT)/TrialNum;
    CountAvg_NIHT=sum(Count_NIHT)/TrialNum;

    MSE_IHT=sum(error_IHT)/TrialNum;
    TimeAvg_IHT=sum(t_IHT)/TrialNum;
    CountAvg_IHT=sum(Count_IHT)/TrialNum;

    MSE_SBL=sum(error_SBL)/TrialNum;
    TimeAvg_SBL=sum(t_SBL)/TrialNum;
    CountAvg_SBL=sum(Count_SBL)/TrialNum;

    MSE_GPSR=sum(error_GPSR)/TrialNum;
    MSE_GPSRdb=sum(error_GPSRdb)/TrialNum;
    TimeAvg_GPSR=sum(t_GPSR)/TrialNum;

    MSE_CoSaMP=sum(error_CoSaMP)/TrialNum;
    MSE_CoSaMPdb=sum(error_CoSaMPdb)/TrialNum;
    TimeAvg_CoSaMP=sum(t_CoSaMP)/TrialNum;

    MSE_SOCP=sum(error_SOCP)/TrialNum;
    TimeAvg_SOCP=sum(t_SOCP)/TrialNum;

    MSE_BCS=sum(error_BCS)/TrialNum;
    TimeAvg_BCS=sum(t_BCS)/TrialNum;
end

MSE_LS=sum(error_LS)/TrialNum;

if TrialNum==1
    if RunOtherMethod
        figure(2)
        set(gcf,'color','white');
        subplot(5,2,1);
        stem(s,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('True signal');
        box off
        subplot(5,2,2);
        stem(x_ExCoV,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('ExCoV estimated signal');
        box off
        subplot(5,2,3);
        stem(x_ExCoVapp,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('Approx. ExCoV estimated signal');
        box off    
        subplot(5,2,4);
        stem(x_NIHT,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('Normalized IHT estimated signal');
        box off    
        subplot(5,2,5);
        stem(x_IHT,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('IHT estimated signal');
        box off    
        subplot(5,2,6);
        stem(x_SBL,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('SBL estimated signal');
        box off    
        subplot(5,2,7);
        stem(x_GPSR,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('GPSR estimated signal');
        box off
        subplot(5,2,8);
        stem(x_SOCP,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('SOCP estimated signal');
        box off
        subplot(5,2,9);
        stem(x_CoSaMPdb,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('CoSaMP-DB estimated signal');
        box off
        subplot(5,2,10);
        stem(x_BCS,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('BCS estimated signal');
        box off
    else
        figure(2)
        set(gcf,'color','white');
        subplot(2,1,1);
        stem(s,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('True signal');
        box off
        subplot(2,1,2);
        stem(x_ExCoV,'Marker','none','linewidth',2,'color','k');
        axis([1,m,-1.5,+1.5]);
        title('ExCoV estimated signal');
        box off
    end
else
    clear H HA Phi Psi
    save([SignalType_Str,'_TrialNum',num2str(TrialNum),'_N',num2str(N),'_T',num2str(T),'_sigma2-',num2str(sigma2),'.mat']);
end

