%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This script is used to obtain the figure 4 in the paper
%%%% "Learning with Dynamic Group Sparsity"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear; close all;
%%% parameters
RunningNum=10;%%% The whole random running number, which is 100 in our paper.
H=512; W=1; n=H*W; P=1; K=64; q=4; 

m_num=round(K*[2.5 2.75 3, 3.25, 3.5, 3.75 4]);

for j=1:length(m_num),
    m=m_num(j)
for iter=1:RunningNum,
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% step 1: randomly create DG(K, q)-sparse 1D-data
    %%% randomly distribute the K nozero items into q groups
    %%% randomly choose location of each group and keep them disjoined
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x=zeros(n,1);[H, W]=size(x);
    r=abs(randn(q,1)); r=r+1; r=round(r*K/sum(r)); r(q)=K-sum(r(1:q-1));
    g=round(r*n/K); g(q)=n-sum(g(1:q-1));
    g_cum=cumsum(g);
    for i=1:q,        
        loc=randperm(g(i)-r(i)-2);        
        x_tmp=zeros(g(i), 1);
        x_tmp(loc(1)+2:loc(1)+1+r(i))=sign(randn(r(i), 1));
        x(g_cum(i)-g(i)+1:g_cum(i), 1)=x_tmp;        
    end
%     figure;plot(x);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% step 2: Creating random projection matrix and measurements
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:P
        A=randn(m,n);item=sqrt(sum(A.*A,2));A=A./repmat(item, [1,n]);
        AA{i}=A; e=randn(m,1); %e=e/norm(e(:));
        y{i}=AA{i}*x(:)+0.01*e;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% step 3: Sparse recovery
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%% Method:GPSR     
    tic;
    for i=1:P,
            tau=0.01*max(abs(AA{i}'*y{i}));tolA = 1.e-2;
            [xest,x_debias,objective,times,debias_start,mses]= ...
            GPSR_BB(y{i},AA{i},tau,'Debias',1,'Monotone',1,'StopCriterion',3,'ToleranceA',tolA,'Continuation',1,'Verbose',0);
            if ~isempty(x_debias) xest=x_debias;  end
            xest1(:,:,i)=reshape(xest, [H, W]);
    end
    time1=toc;    
    
    
    %%%% Method 2: SPGL1-LASSO
    tic;
    opts = spgSetParms('verbosity',0);
%     [xest,r,g,info] = spg_lasso(AA{1},y{1},tau/4,opts );
    [xest,r,g,info] = spgl1(AA{1},y{1},[],[],[],opts);
    xest2   = xest(1:n,1);
    time2=toc;
    
    
    %%%%  Method 3: OMP
    tic;
    [xest, norm_save] = OMP_CS(K,AA{1},y{1});
    xest3=xest(1:n,1);
    time3=toc;

    
    
    %%%%  Method:SP
    tic;
    Rec = DGS_CS(K, AA, y, 1, 0*ones(n,1), H, W, 'None');
    time4=toc; xest=Rec.x_hat; xest4=xest(1:n,1);

    %%%%  Method:DGSR
    tic;
    Rec = DGS_CS(K, AA, y, 1, 0.5*ones(n,1), H, W, 'GridCol');
    time5=toc; xest=Rec.x_hat; xest5=xest(1:n,1);
    
    diff1=norm(xest1(:)-x(:),2)/norm(x(:),2);
    diff2=norm(xest2(:)-x(:),2)/norm(x(:),2);
    diff3=norm(xest3(:)-x(:),2)/norm(x(:),2);
    diff4=norm(xest4(:)-x(:),2)/norm(x(:),2);
    diff5=norm(xest5(:)-x(:),2)/norm(x(:),2);
   

    result(:,iter)=[diff1,diff2,diff3, diff4, diff5, time1 time2,time3, time4, time5];

%     figure; subplot(4,1,1);plot(x); title('(a) Original Signal')
%     subplot(4,1,2); plot(xest1); title('(b) GPSR')
%     subplot(4,1,3); plot(xest3); title('(d) SP')
%     subplot(4,1,4); plot(xest4); title('(e) DGSR')
end
CompareResults(:,j)=mean(result, 2);
CompareResultsStd(:,j)=std(result,0, 2);
end
save CompareResults CompareResults;
save CompareResultsStd CompareResultsStd;
K=64;
figure; hold on;
errorbar(m_num/1, CompareResults(1,:), CompareResultsStd(1,:), 'bs-.', 'linewidth', 2,'MarkerSize',10); 
errorbar(m_num/1, CompareResults(2,:), CompareResultsStd(2,:), 'gd--', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(3,:), CompareResultsStd(3,:), 'cp-.', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(4,:), CompareResultsStd(4,:), 'mv--', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(5,:), CompareResultsStd(5,:), 'rh-', 'linewidth', 2,'MarkerSize',10);
ylabel('Recovery Error'); xlabel('Measurement Number');box on;
legend('GPSR', 'SPGL1-LASSO','OMP','SP', 'DGS')
textobj = findobj('type', 'text');
set(textobj, 'fontsize', 12);
h_xlabel = get(gca,'XLabel');
set(h_xlabel,'FontSize',12); 
h_xlabel = get(gca,'YLabel');
set(h_xlabel,'FontSize',12); 
axis([m_num(1)-K/8 m_num(end)+K/8 -0.2 1.6]);

figure; hold on;
errorbar(m_num/1, CompareResults(6,:), CompareResultsStd(6,:), 'bs-.', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(7,:), CompareResultsStd(7,:), 'gd--', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(8,:), CompareResultsStd(8,:), 'cp-.', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(9,:), CompareResultsStd(9,:), 'mv--', 'linewidth', 2,'MarkerSize',10);
errorbar(m_num/1, CompareResults(10,:), CompareResultsStd(10,:), 'rh-', 'linewidth', 2,'MarkerSize',10);
ylabel('Running Time (Second)'); xlabel('Mesurement Number');box on;
legend('GPSR', 'SPGL1-LASSO','OMP','SP', 'DGS');
textobj = findobj('type', 'text');
set(textobj, 'fontsize', 12);
h_xlabel = get(gca,'XLabel');
set(h_xlabel,'FontSize',12); 
h_xlabel = get(gca,'YLabel');
set(h_xlabel,'FontSize',12); 
axis([m_num(1)-K/8 m_num(end)+K/8 0 2.5]);

