function [X , S , r] = BLOOMP(A,b,s,band,neighbor)
%% BLOOMP ''Coherence Pattern-Guided Compressive Sensing with Unresolved
%% Grids''
%% input: A:        sensing matrix
%%        b:        measurement
%%        band:     radius of excluded band in Band Exclusion
%%        neighbor: radius of searching band in Local Optimization
%%        ep:       noise level Ax + e = b where ||e||_2 \le ep. By
%%                  default, ep = 1e-12.
%%        search:   search == 1: locally optimize the object one by one
%%                  in an increasing order; search == 2: locally optimize the object
%%                  in a random order. By default, search == 1.
%% output: X:       output vector
%%         S:       support
%%         r:       A*X-b


[m,N] = size(A);
itermax = s;
%%

r = b;

S = [];
Region = 1:N;

for iter=1:itermax
    % step1: pick the next index
    mu = abs(A(:,Region)'*r);
    [q,lambda] = max(mu);
    NewIndex = Region(lambda);
    
    % step2 : add the new index to the support
    S = [S , NewIndex]; 
    

    %local search
    [X , S , r] = LS(A,b,neighbor,iter,S);
    
    % update the complement of the excluded band
     Region = 1 : N;
     for k = 1 : iter
         OffRegion = S(k)-band:S(k)+band;
         if OffRegion(1)<0.5
             Ind = OffRegion<0.5;
             OffRegion = OffRegion.*(1-Ind)+(OffRegion+N).*Ind;
         elseif OffRegion(end)>(N+0.5)
             Ind = OffRegion>(N+0.5);
             OffRegion = OffRegion.*(1-Ind)+(OffRegion-N).*Ind;
         end
         Region    = setdiff(Region,OffRegion);
     end
     
%    if(res<ep) 
%        return;
%    end
end

S = S';



