%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This script is used to obtain the figure 5 in the paper
%%%% "Learning with Dynamic Group Sparsity"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear; close all;
%%%%%% The whole random running number, which is 100 in our paper.
RunningNum=100; %%%%   you can change it

H=48; W=48; n=H*W; P=3; 
load x;  %%%% the original sparse color image
m_num=160*[2.25  2.5, 2.75, 3, 3.25, 3.5, 3.75, 4];

figure; imshow(x, []); 
K=length(find(x(:)>0))/P;

for j=1:length(m_num),
     j
for iter=1:RunningNum,
    m=m_num(j);   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% step 2: Creating random projection matrix and measurements
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x=reshape(x, [n, P]);
    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(:,i)+0.01*e;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% step 3: Sparse recovery
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    %%%%  Method 1: MCS    
    tic;
    a = 1e2/0.01; b = 1;
    [xest, ml] = mt_CS(AA,y,a,b,1e-8);   
    xest1=reshape(xest, [H, W, P]);
    time1=toc;    
    
    %%%%  Method 2: SP
    tic;
    for i=1:P,
        Rec = DGS_CS(K, AA{i}, y{i}, 1, 0.0*ones(n,1), H, W, 'None');
        xest=Rec.x_hat; xest2(:,:,i)=reshape(xest, [H, W]);
    end   
    time2=toc; 
% 
    %%%%  Method 3: DGS
    tic;
    Rec = DGS_CS(K, AA, y, 1, 0.5*ones(n,1), H, W, 'Grid');
    time3=toc; xest=Rec.x_hat; xest3=reshape(xest, [H, W, P]);
%     
    %%%% Normalized Recovery Error
    diff1=norm(xest1(:)-x(:),2)/norm(x(:),2);
    diff2=norm(xest2(:)-x(:),2)/norm(x(:),2);
    diff3=norm(xest3(:)-x(:),2)/norm(x(:),2);

    result(:,iter)=[diff1,diff2,diff3,time1 time2,time3];
    x=reshape(x, [H, W, P]);
%     figure; subplot(1,4,1);imshow(x, []); title('(a)'); 
%     subplot(1,4,2); imshow(xest1, []); title('(b)');
%     subplot(1,4,3); imshow(xest2, []); title('(c)');
%     subplot(1,4,4); imshow(xest3, []); title('(d)');
end

CompareResults(:,j)=mean(result, 2);
CompareResultsStd(:,j)=std(result,0, 2);
end 

save CompareResults CompareResults;
save CompareResultsStd CompareResultsStd;

figure; hold on;
errorbar(m_num, CompareResults(1,:), CompareResultsStd(1,:), 'bv-.', 'linewidth', 2); 
errorbar(m_num, CompareResults(2,:), CompareResultsStd(2,:), 'gd--', 'linewidth', 2);
errorbar(m_num, CompareResults(3,:), CompareResultsStd(3,:), 'rs-', 'linewidth', 2);
ylabel('Recovery Error'); xlabel('Measurement Number'); box on;
legend('MCS','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)-20 m_num(end)+20 -0.1 1.3]);

figure; hold on;
errorbar(m_num, CompareResults(4,:), CompareResultsStd(4,:), 'bv-.', 'linewidth', 2);
errorbar(m_num, CompareResults(5,:), CompareResultsStd(5,:), 'gd--', 'linewidth', 2);
errorbar(m_num, CompareResults(6,:), CompareResultsStd(6,:), 'rs-', 'linewidth', 2);
ylabel('Running Time (Second)'); xlabel('Mesurment Number'); box on;
legend('MCS', '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); 