function retLEV = CalcSLEV(inZ,inX,inTP,inOption)
%function retLEV = CalcSLEV(inZ,inX,inTP,inOption)
% 
% Calculate the leading sparse eigenvecor of X'ZZ'X.
%
% Input:
%   inZ: Z, which can be multivariate
%   inX: X
%   inTP: tuning parameter (lambda)
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%       ptype -- type of penalization:
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation,
%                'l1' ~ L1,
%                'scad' ~ SCAD (gamma=3.7 by default),
%                'mcp' ~ MCP (gamma=3.7 by default);
%       gamma -- parameter for SCAD or MCP (=3.7 by default);
%
% Output:
%   retLEV: the leading sparse eigenvector of X'ZZ'X.
%
% Dabao Zhang     May 30, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(4, 4, nargin));

% Initialize retRes
p = size(inX,2);
omega = zeros(p,0);
tmpMat = inZ'*inX;

% Calculate the leading eigenvector of X_j'YY'X_j
[tmpU,tmpS,omega] = svds(tmpMat,1);

% Calculate the regularized component vector
idIter = 0;
omegapre = omega;   omegadf = 1;  omeganrm = norm(omega);
alpha = tmpMat'*(tmpMat*omega); alpha = alpha/norm(alpha);
while((omeganrm>eps)&&(omegadf>inOption.tol)&&(idIter<=inOption.maxit))
    idIter = idIter + 1;
    omegapre = omega;
    omega = tmpMat'*(tmpMat*alpha);
    
    switch( inOption.ptype )
        case 'ebt'  % POCRE0 via Empirical Bayes Thresholding
            tmpSD = median(abs(omega))*inTP/norminv(0.75,0,1);
            tmpRes = GEBCauchy(omega,tmpSD,[1,1,1]);
            omega = tmpRes.mu;
            
        case 'ebtz'  % POCRE via Empirical Bayes Thresholding
            omega = omega/norm(omega);
            omega = .5*log((1+omega)./(1-omega));
            tmpSD = inTP/sqrt(length(omega)-3);
            tmpRes = GEBCauchy(omega,tmpSD,[1,1,1]);
            omega = 1-2./(exp(2*tmpRes.mu)+1);
            
        case 'l1'
            tmpSD = median(abs(omega))*inTP/norminv(0.75,0,1);
            omega = max(abs(omega)-tmpSD*sqrt(2*log(p)),0)...
                    .*sign(omega);
            
        case 'scad'
            if( ~isfield(inOption,'gamma') ) inOption.gamma = 3.7; end
            tmpSD = median(abs(omega))*inTP/norminv(0.75,0,1);
            tmpT = tmpSD*sqrt(2*log(length(omega))); % Threshold of L1

            % |z|<=2*lambda
            tIdx = find(abs(omega)<=2*tmpT);
            omega(tIdx) = max(abs(omega(tIdx))-tmpT,0).*sign(omega(tIdx));

            % 2*lambda<|z|<=gamma*lambda
            tIdx = setdiff(find(abs(omega)<=inOption.gamma*tmpT),tIdx);
            tmpSO = max(abs(omega(tIdx))-tmpT/(1-1/inOption.gamma),0);
            omega(tIdx) = tmpSO.*sign(omega(tIdx))/(1-1/(inOption.gamma-1));
            
        case 'mcp'
            if( ~isfield(inOption,'gamma') ) inOption.gamma = 3.7; end
            tmpSD = median(abs(omega))*inTP/norminv(0.75,0,1);
            tmpT = tmpSD*sqrt(2*log(length(omega))); % Threshold of L1
            tIdx = find(abs(omega)<=inOption.gamma*tmpT);
            tmpSO = max(abs(omega(tIdx))-tmpT,0);
            omega(tIdx) = tmpSO.*sign(omega(tIdx))/(1-1/inOption.gamma);
            
        otherwise
            error('gPOCRE: The penalization method is not supported!');
    end

    alpha = tmpMat'*(tmpMat*omega);
    omeganrm = norm(omega);
    if( omeganrm>eps )
        alpha = alpha/norm(alpha);  % norm(omega)>eps ==> norm(alpha)>0?
        omegadf = norm(omegapre-omega)/omeganrm;
    end
end

if( omeganrm>eps )
    retLEV = omega/omeganrm;
else
    retLEV = omega;
end
