
function varargout = IterGS(yt, gram, gramBias, allCG, sp, uu, nm, q0, scale,maxIter, stdThresh, betaI)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% itGS.m: Run the iterative Graphlet Screening procedure.
%%%
%%% INPUT:
%%% yt   X'Y
%%% gram: the thresholded gram matrix
%%% gramBias: the trimed part in the threholding.
%%% allCG:  all the connected connected subgraphs of gram with size no more than nm.
%%% sp: a user-friendly version of the tuning paramters v(vartheta) of the
%%% graphlet screening. sp is the expected sparsity level.
%%% uu: a user-friendly version of the tuning paramters v(vartheta) of the
%%% graphlet screening. uu can be interpreted as the minimal detected signal
%%% strength.
%%% nm: the maximal size of the connected subgraphs considered in the
%%% screening step.
%%% q0: the minimal threshold in the screening step when the signial is very weak.
%%% scale: Optional parameter of the screening step. The default is 1.
%%% maxIter: Optional, maximal iterations of the refinement. The default is 3.
%%% stdThresh: Optional, the threshold on the std(w) of stopping the loop.
%%% The default is 1.05.
%%% betaI: the initial estimate of the regression coefficient beta, the
%%% default is uu.*sign(yt).*(abs(yt) > uu).
%%%
%%% OUTPUT ARGUMENTS:
%%% betaGS:  The GS estimate for beta. 
%%% nIteration: Optional output, the number of iterations before it stops.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Input defaults.
if nargin < 7
    error('Insufficient number of input arguments')   
elseif nargin == 7
    q0 = 0.1;
    scale = 1;
    maxIter = 3;
    stdThresh = 1.05;
    betaI = uu.*sign(yt).*(abs(yt) > uu);
elseif nargin == 8
    scale = 1;
    maxIter = 3;
    stdThresh = 1.05;
    betaI = uu.*sign(yt).*(abs(yt) > uu);
elseif nargin == 9
    maxIter = 3;
    stdThresh = 1.05;
    betaI = uu.*sign(yt).*(abs(yt) > uu);
elseif nargin == 10
    stdThresh = 1.05;
    betaI = uu.*sign(yt).*(abs(yt) > uu);
elseif nargin == 11
    betaI = uu.*sign(yt).*(abs(yt) > uu);
end

%%
[p0, p] = size(gram);
r = uu^2/2/log(p);
v = 1-log(sp)/log(p);
lambda = sqrt(2*v*log(p));

%% The initial values
betaGS = betaI;
w = yt;

%% Iterative refinement
for it = 1:maxIter
    %%% prepare
    last_w_std = std(w);
    last_beta = betaGS;
    last_nonzero = (last_beta ~= 0);
    %%% update w
    w = yt -gramBias(:,last_nonzero) * last_beta(last_nonzero);
    if (std(w) > stdThresh * last_w_std)
        %%% abort the current iteratoin if std(w) is not improved much
        nIteration = it - 1;
        break
    end
    %%% update betaGS
    survivor = ScreeningStep(w, gram, allCG, nm, v, r, scale);
    betaGS = CleaningStep(survivor, w, gram, lambda, uu);
    nIteration = it; % count finished iterations
end % for it

%% define the output with various length
varargout = cell(nargout,1);
if nargout == 1
    varargout = {betaGS};
elseif nargout == 2
    varargout = {betaGS, nIteration};
end
end

