function [x, e,norm_save] = DGS_BGColor(A, y, x0, e0, Paras)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Multiple Tasks/Groups sparsity
%%% Input:
%%%   A:   projection matrix. Cell structure, One cell for one task.
%%%   y:   CS measurements. Cell structure, One cell for one task.
%%%   Jb:  the group size for each task
%%%   K:   the sparsity  group number
%%%   beta   betaM=alfa*ones(nc*k, 5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Variable:
%%%  Jy =length(y), the number of tasks. 
%%%  [dJy, TJb]=size(A{Jy})
%%%  Jx=Jb*Jy, the total group size for all tasks
%%%  Jx*K: Total sparsity number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Output:
%%%   Rec.x_hat [TJb, Jy]: the reconstructed sparse variable
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  Written by Junzhou Huang, Rutgers University, Nov 6, 2008 ver1.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iscell(y)  
    Jy = length(y); 
else Jy = 1; 
    A = {A}; y = {y}; 
end
m=size(A{1}, 1); T=size(A{1}, 2);
Jb=length(Paras); TJb=T+m;  %%%% Jb must be 2
in = 1; KJb=0; K=Paras{1}.k; Ke=Paras{2}.k;


for jjb=1:Jb,
    KJb=KJb+Paras{jjb}.k;
    index_T{jjb}=[1:Paras{jjb}.T]';
    if jjb==1
        index_N{jjb}=Edge4Index(Paras{jjb}.H, Paras{jjb}.W);
    else
        Tau=Paras{jjb}.Tau;
%         index_N{jjb}=Edge8Index(Paras{jjb}.H, Paras{jjb}.W);
        index_N{jjb}=Edge16IndexSimple(Paras{jjb}.H, Paras{jjb}.W, Tau);
    end
end



if isempty(x0)
    y_r = y;

for jj=1:Jy,
    a= y_r{jj}'*A{jj}; 
    cv((jj-1)*T+1:jj*T, 1) =a';
    cve((jj-1)*m+1:jj*m, 1) = y_r{jj};
end

cv_index=GkqPruneRobust(cv, cve, index_T, index_N, Jy, Paras);

Index_save{in} = cv_index;
for jj=1:Jy,
    [xs, xe, er]=Lsq_Occlusion(A{jj}, y{jj}, cv_index);
    
    x_p((jj-1)*K+1:jj*K, 1)=xs;
    e_p((jj-1)*Ke+1:jj*Ke, 1)=xe;
%     [x_p((jj-1)*KJb+1:jj*KJb, 1), er]=Lsq_Occlusion(A{jj}, y{jj}, cv_index);
    y_r{jj}=er; norm_yr(jj) =norm(y_r{jj});
end


else  %%%% given x0
    jj=1;
    cv_index{1}=find(x0((jj-1)*T+1:jj*T, 1)~=0);
    cv_index{2}=find(e0((jj-1)*m+1:jj*m, 1)~=0);   
    for jj=1:Jy,
        y_r{jj}=y{jj}-A{jj}*x0((jj-1)*T+1:jj*T, 1)-e0((jj-1)*m+1:jj*m, 1);
        norm_yr(jj) =norm(y_r{jj});
        xtmp=x0((jj-1)*T+1:jj*T, 1);
        etmp=e0((jj-1)*m+1:jj*m, 1);
        x_p((jj-1)*length(cv_index{1})+1:jj*length(cv_index{1}), 1)=xtmp(cv_index{1});
        e_p((jj-1)*length(cv_index{2})+1:jj*length(cv_index{2}), 1)=etmp(cv_index{2}); 
    end
    Index_save{in} = cv_index;     
end

x_curr{1}=x_p; e_curr{1}=e_p;
norm_save(in) = norm(norm_yr);


while 1
    in = in+1;
    for jj=1:Jy,
        a= y_r{jj}'*A{jj}; 
        cv((jj-1)*T+1:jj*T, 1) =a';
        cve((jj-1)*m+1:jj*m, 1) = y_r{jj};
    end
    cv_index=GkqPruneRobust(cv, cve, index_T, index_N, Jy, Paras);
%     cv_index=GkqPruneRobust(cv', index_T, index_N, Jy, Paras);
    for jjb=1:Jb,
        a1=Index_save{in-1}{jjb}; a2=cv_index{jjb};
        cv_add_tmp=unique([a1;a2], 'rows'); 
        cv_add{jjb}=cv_add_tmp(:);
    end
   
    Kxp=length(cv_add{1});Kep=length(cv_add{2});
    x_p=zeros(Kxp*Jy, 1); e_p=zeros(Kep*Jy, 1);
    for jj=1:Jy,
        [xp, xep, er]=Lsq_Occlusion( A{jj}, y{jj}, cv_add);
%         [xp3, xep3, er3]=Lsq_Occlusion3(A{jj}, y{jj}, cv_add);
        x_p((jj-1)*Kxp+1:jj*Kxp, 1)=xp;
        e_p((jj-1)*Kep+1:jj*Kep, 1)=xep;
    end   
    
    cv_index=GkqPruneRobust(x_p, e_p, cv_add, index_N, Jy, Paras);
    
%     e=zeros(8064,1);
%     e(Index_save{in}{2})=e_curr;
%     im=reshape(e ,[96 84]);
%     figure; imshow(im, []);

    Index_save{in}=cv_index;
    
    x_p=zeros(Jy*K,1);e_p=zeros(Jy*Ke,1);
    for jj=1:Jy,
        [xs, xe, er] = Lsq_Occlusion( A{jj}, y{jj}, cv_index);
        x_p((jj-1)*K+1:jj*K, 1)=xs;
        e_p((jj-1)*Ke+1:jj*Ke, 1)=xe;
        y_r{jj}=er; norm_yr(jj) =norm(y_r{jj});
    end
        
    norm_save(in) = norm(norm_yr);
    if ( norm_save(in) == 0 ) ||(norm_save(in)>=norm_save(in-1) )
        break;
    end
    x_curr{1}=x_p; e_curr{1}=e_p;    
end
x_hat = zeros(TJb*Jy,1);

sup1=[]; sup2=[];
for jj=1:Jy,
    temp1=Index_save{in-1}{1}+T*(jj-1);
    temp2=Index_save{in-1}{2}+m*(jj-1);
    sup1=[sup1; temp1]; sup2=[sup2; temp2];
end

e=zeros(Paras{2}.T*Jy,1); 
x=zeros(Paras{1}.T*Jy,1); 
% for i=1:Jy,
    x(sup1)=x_curr{1};
    e(sup2)=e_curr{1};
% end

return

function gs_index=GkqPruneRobust(x, xe, xindex, nindex, Jy, Paras);
%%%% x=[X1; X2; ...; XJ], XJ-> T x 1   
%%%% x [TJx, 1]  xindex [T, 1]
%%%% K is sparse for T, thus total sparsity is K*Jx for T*Jx

TJy=length(x); mJy=length(xe);
Ts=TJy/Jy; ms=mJy/Jy; Jb=length(Paras);
XM{1}=reshape(x, [Ts, Jy]);
XM{2}=reshape(xe, [ms, Jy]);
T_curr=0; m_curr=0; 

for jjb=1:Jb,
    T=Paras{jjb}.T;
    betaM=Paras{jjb}.beta;
    K=Paras{jjb}.k;
    mm=length(xindex{jjb});
%     XMB=XM{jjb}(m_curr+1:m_curr+mm, :);
    XMB=XM{jjb};
    xtmp=zeros(T,Jy); 
    index_curr=xindex{jjb};
    xtmp(index_curr, :)=XMB;
    clear tmlp;
    for ne=2:size(nindex{jjb},2) 
        tmp=nindex{jjb}(index_curr,ne);
        for p=1:Jy,
            tmlp(:,p)=betaM(index_curr, ne-1).*xtmp(tmp,p);
        end
%         XMB=[XMB, betaM(index_curr, ne-1).*xtmp(tmp,:)];  %%% 03/04
        XMB=[XMB, tmlp];  
    end
    XMJnorm=sqrt(sum(XMB.*XMB, 2));
    [xmj_sort, xmj_index] = sort(XMJnorm,'descend');
    xmj_index = sort( xmj_index(1:K) ); 
    xmj_index2 = xindex{jjb}( xmj_index);
    xmj_index2 = sort(xmj_index2 );
    temp=xmj_index2;  %%%%%% modify later
    gs_index{jjb}=xmj_index2;
    T_curr=T_curr+T; m_curr=m_curr+mm;
end
return

function Edge4=Edge4Index(H, W);
N = H*W;
Edge4=zeros(N,5);
%% Self
Edge4(:,1)=[1:N;];
%% down
is1=Edge4(:,1)+1;
is1(H:H:N,1)=[H-1:H:N-1;];
Edge4(:,2)=is1;
%% up
is2=Edge4(:,1)-1;
is2(1:H:N-H+1,1)=[2:H:N-H+2;];
Edge4(:,3)=is2;
%% LEFT
is=Edge4(:,1)-H;
is(1:H,1)=[H+1:2*H;];
Edge4(:,4)=is;
if W==1
   Edge4(:,4)=is1;
end
%% RIGHT
is=Edge4(:,1)+H;
is(N-H+1:N,1)=[N-2*H+1:N-H;];
Edge4(:,5)=is;
if W==1
   Edge4(:,5)=is2;
end
return