%%% ScreeningStep of the function
function indicator = ScreeningStep(yt, gram, allCG, nm, v, r, q0, scale)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% This function conduct the GS-step of the graphet screening
%%%
%%% INPUT:
%%% yt: X'Y where X and Y are the predictor matrix and the observed response
%%% gram: the regularized gram matrix
%%% allCG: a structure with one field cg, cg contains all the connected
%%% subgraphs of gram with no more than nm nodes.
%%% nm: the maximal subgraph size to be considered in the screening step
%%% v: vartheta, tuning parameter tied to the sparse level
%%% r: tuning parameter tied to the signal strength
%%% q0: the minimal threshold when the signial is very weak.
%%% scale: q(D,F)=q^{max}(D,F)*scale, default is scale=1
%%% 
%%% The outcomes are
%%% indicator: a p by 1 indicator vector, where 1 means retains as a
%%% protential signal.
%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 6
    error('Too few input arguments')
elseif nargin == 6
    q0 = 0.1;
    scale = 1;
elseif nargin == 7
    scale = 1;
end

if length(allCG)<nm
    error('The provided connected subgraphs are not enough')
end

%% options for quadprog
if  verLessThan('matlab', '7.12.0')
    options = optimset('algorithm', 'interior-point', 'LargeScale','off','display', 'off'); 
elseif  verLessThan('matlab', '7.14.0')
    options = optimset('Algorithm', 'interior-point-convex', 'LargeScale','off','display', 'off');
else
    options = optimset('Algorithm','interior-point-convex','Display','off');
end


%% univeriate screening, identical to the screening step of UPS
p = length(yt);

q1 = (v+r)^2/r^2/4;
tau = sqrt(2*r*log(p));
tau1 = sqrt(q1)*tau*sqrt(scale);
indicator = (abs(yt)>tau1);
indices = find(indicator);

%% multivariate screening
if nm > 1
    for ii = 2:nm
        %%by our theory, we do not need to consider the clusters longer than 1/v
        if ii > 1/v
            break;
        end
        %%extract all the CG with ii nodes in the graph
        iCG = allCG(ii).cg;
        %%if there is no such CG, break
        if length(iCG) == 0
            break;
        else
            indicatorii = zeros(size(indicator));
            [nii, nk] = size(iCG);
            for jj = 1:nii
                %%find the indices to be tested.
                ds = setdiff(iCG(jj, :), indices);
                d = length(ds);
                if d > 0
                    %%define the test statistic.
                    es = intersect(iCG(jj, :), indices);
                    if length(es) > 0
                        gramD = gram(ds, ds)-gram(ds, es)*inv(gram(es, es))*gram(es, ds);
                    else
                        gramD = gram(ds, ds);
                    end
                    if d == 1
                        ww = gramD;
                    elseif d == 2
                        ww = d*min(eig(gramD));
                    else
                        wwV = zeros(1, 2^(d-1));
                        %%% for a length d vector, let the sign of the first coordinate to be 1, there are 2^(d-1) combinations.
                        %%% by symmetry of this quadratic form, the the above assumption is ok.
                        for kk = 0:(2^(d-1)-1)
                            idx = VectorizeBase(kk, 2, d-1); % express k on base 2 as a row vector
                            signs = [1, 2.*idx - 1]'; % convert to -1, 1 to indicate the signs
                            [bt, Lt] = quadprog(full(gramD), zeros(d, 1), -diag(signs), -ones(d, 1),[],[],[],[],1.1*signs, options);
                            %[bt, Lt] = quadprog(full(gramD), zeros(d, 1), -diag(signs), -ones(d, 1));
                            wwV(kk+1) = 2*Lt;
                        end
                        ww = min(wwV);
                    end
                    if ww > (v/r)*(2*d-(d-sqrt(d^2-1))*((d/2) ~= floor(d/2)))
                        core = ww*r;
                        delta = v^2/core/4-sqrt((core^2-2*core*d*v+v^2)*(core^2-2*core*d*v+v^2>0))+sqrt((core^2-2*core*d*v)*(core^2-2*core*d*v>0));
                        threshold = scale*2*log(p)*(core*5/4-v*d/2-sqrt((core^2-2*core*d*v)*(core^2-2*core*d*v>0))+delta*(d/2~=floor(d/2)));
                    else
                        threshold = q0;
                    end
                    test = yt(iCG(jj,:))'*inv(gram(iCG(jj,:),iCG(jj,:)))*yt(iCG(jj,:))-yt(es)'*inv(gram(es,es))*yt(es);
                    if test > threshold
                        %%% as long as a node pass in one test, it will be included in
                        %%% the survival indices.
                        indicatorii(ds) = 1;
                    end
                end                
            end
           indicator = max(indicator, indicatorii);
           indices = find(indicator);
        end
    end
end
%% the output must be a logical vector
indicator = logical(indicator);
end


%% a private function, transform a non-negative integer to a vector on a given base.
function vector = VectorizeBase(i, base, vectorLength)

if (i < 0)||(base < 2)||(vectorLength < 1)
    error('Unreasonable arguments')
end

if(i > base^vectorLength -1)
    warning('ScreeningStep:vectorizeBase:IncompatibleArguments',' the vector is not long enough')
end

vector = zeros(1, vectorLength);
for j = vectorLength:-1:1
    vector(j) = i - base *floor(i/base); % remainder
    i = (i - vector(j))/base;
end

end
% end private function vectorizeBase
