function [retX,retSIdx,retRes] = gPOCREScreen(inY,inX,inFamily,inLink,inEIdx,inOption)
%function [retX,retSIdx,retRes] = gPOCREScreen(inY,inX,inFamily,inLink,inEIdx,inOption)
% 
% Using the generalized Penalized Orthogonal-Component Regression (gPOCRE)
% toscreen the variables for fitting a generalized linear model.
%
% Model: g(E[Y]) = mu + X*beta = mu + \sum_j (X*varpi_j)*vartheta_j'
%                              = mu + \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: Y response values of original data. It includes observed survival
%        time (col 1) & censoring status (col 2) for survival data.
%   inX: X;
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   inEIdx: Indices of exceptional variables which are always included in
%           the model (but not counted into maxvar);
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%       bctype -- type of bias correction: 0 ~ no bias correction,
%                 1 ~ approximate (default), type 2 -- optimal;
%       maxcmp -- maximum number of components;
%       maxvar -- maximum number of selected variables;
%       W -- prespecified weight vector for the weight matrix;
%
% Output:
%   retX: inX[,[retEIdx,retSIdx]];
%   retSIdx: indices of selected variables;
%   retRes: nCmp, mu, beta, deviance, W, vartheta, omega, varpi, xomega,
%           ([X_1*omega_1, ...] = X*varpi), theta.
%
% Notes on suggested Family & Link combinations,
%   'normal': identity (default), log, -1
%   'binomial': logit (default), probit, comploglog
%   'multibc' (multinomial baseline category model): bclogit (default) 
%   'multipo' (multinomial proportional odds model): pologit (default)
%   'multicr' (multinomial continuation ratio model): crlogit (default)
%   'poisson': log (default), identity, 0.5
%   'gamma': -1 (default), loglog, identity, log
%   'inverse gaussian': -2 (default), -1, identity, log
%   'cox': Cox proportinal hazards model for survival data (no link)
%   'quasi': logit (default), probit, comploglog, identity, log, loglog,
%            -1,-2,0.5, customized function
%
% Dabao Zhang    July 26, 2016
%
% Revised by:
%       Dabao Zhang, November 7, 2016 -- Default value of inOption.maxvar.
%
% Copyright (c) 2016 by Dabao Zhang.

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

% Special for multinomial models & Cox proportional odds model
if( sum(strcmp(inFamily,{'multibc','multipo','multicr'})) )
    if( sum(abs(sum(inY,2)-1))<eps )
        inY = inY(:,2:end);
    end
elseif( strcmp(inFamily,'cox') )
    inOption.bctype = 0; % No bias correction!

    % Prepares for fast calculation: inY to [ObsTime CStatus RIdx CIdx]
    % Sort the data
    if( size(inY,2)<=2 )
        inCens = inY(:,2);
        inRC = inY(:,3:4);
        inY = inY(:,1);
    else
        % Sort the data
        inCens = inY(:,2);
        [inY,tmpIdx] = sort(inY(:,1));
        inX = inX(tmpIdx,:);
        inCens = inCens(tmpIdx);

        % Calculate the rank of failture times for fast calculation
        [~,idxR] = ismember(inY,flipud(inY),'legacy');
        idxR = length(inY) + 1 - idxR;
        [~,idxC] = ismember(inY,flipud(inY));
        idxC = length(inY) + 1 - idxC;
        inRC = [idxR idxC];
    end
end

[n,p] = size(inX);    % X with n objects p variables;
k = size(inY,2);    % Y with n objects k variables;
if( n~=size(inY,1) )
    error('    gPOCRE: X and Y should have the same number of rows!');
end

if( nargin<6 ) inOption = []; end
if( ~isfield(inOption,'tol') ) inOption.tol = 1e-6; end
if( ~isfield(inOption,'maxit') ) inOption.maxit = 100; end
if( ~isfield(inOption,'bctype') ) inOption.bctype = 1; end
if( ~isfield(inOption,'maxcmp') ) inOption.maxcmp = 50; end
if( ~isfield(inOption,'maxvar') ) inOption.maxvar = min(n/log(n),p/2); end
inOption.maxvar = round(inOption.maxvar);
if( ~isfield(inOption,'W') ) inOption.W = []; end
if( ~isfield(inOption,'ptype') ) inOption.ptype = 'ebtz'; end
if( ~isfield(inOption,'covidx') ) inOption.covidx = []; end

if( isempty(inLink) )
    switch inFamily
        case 'normal'
            inLink = 'identity'; 
        case 'binomial'
            inLink = 'logit';
        case 'multibc'
            inLink = 'bclogit';
        case 'multipo'
            inLink = 'pologit';
        case 'multicr'
            inLink = 'crlogit';zc
        case 'poisson'
            inLink = 'log';
        case 'gamma'
            inLink = -1;
        case 'inverse gaussian'
            inLink = -2;
        case 'quasi' % Incomplete
            inLink = 'identity';
        case 'cox'
            inLink = 'cox';
        otherwise
            error('  gPOCRE: The family is not supported...');
    end
end

% Initial values
retRes.nCmp = 0;
if( strcmp(inFamily,'cox') )
    retRes.mu = 0;
else
    retRes.mu = glnk(mean(inY,1),inLink);
end
retRes.beta = zeros(p,k);
if( sum(strcmp(inFamily,{'multipo','multicr'})) )
    retRes.vartheta = zeros(1,0);
else
    retRes.vartheta = zeros(k,0);
end
if( isempty(inOption.W) )
    retRes.W = ones(n,1)/n;
else
    retRes.W = inOption.W/sum(inOption.W);
end
retRes.varpi = zeros(p,0);
retRes.xomega = zeros(n,0);
retRes.omega = zeros(p,0);
retRes.theta = zeros(p,0);
if( size(retRes.W,2)~=1 ) retRes.W = retRes.W'; end

% Define indicator vectors for inEIdx (idcE) & retSIdx (idcS)
idcE = zeros(p,1);
idcE(inEIdx) = 1;
idcS = zeros(p,1);

% Centralize X appropriately according to W:
mX = mean(bsxfun(@times,inX, retRes.W),1);
tmpX = bsxfun(@minus,inX,mX);

% Construct all components
idxCmp = 0;     omeganrm = 1;
nSele = round(inOption.maxvar/inOption.maxcmp);
fprintf(' Screening variables ');
%zeta = CalcZeta4BC(tmpX,retRes.W,inOption.bctype);
zeta = 1-retRes.W;    % sum(inW)=1.
%while( (idxCmp<inOption.maxcmp) && retRes.bSparse && (omeganrm>eps) )
while( (sum(idcS)<inOption.maxvar)&&(idxCmp<inOption.maxcmp) )
    idxCmp=idxCmp+1;

    if( idxCmp==inOption.maxcmp )
        nSele = inOption.maxvar-nSele*(inOption.maxcmp-1);
    end
    
    % Construct the component
    if( strcmp(inFamily,'cox') )
        %tmpRes = gPOCREScreenCox1POC(inY,inCens,inRC,tmpX,retRes,inTP,zeta,inOption);
        tmpRes = gPOCREScreenCox1POC(inY,inCens,inRC,tmpX,retRes,idcE+idcS,nSele,zeta,inOption);
    else
        %tmpRes = gPOCREScreen1POC(inY,tmpX,inLink,retRes,inTP,zeta,inOption);
        tmpRes = gPOCREScreen1POC(inY,tmpX,inFamily,inLink,retRes,idcE+idcS,nSele,zeta,inOption);
    end
    omeganrm = norm(tmpRes.omega);
    idcS = idcS+tmpRes.idcS;

    retRes.mu = tmpRes.mu;
    retRes.vartheta = tmpRes.vartheta;
    retRes.omega(:,idxCmp) = tmpRes.omega;
    retRes.varpi(:,idxCmp) = tmpRes.varpi;
    retRes.xomega(:,idxCmp) = tmpX*tmpRes.omega;
    retRes.theta(:,idxCmp) = tmpRes.theta;
    retRes.nCmp = idxCmp;

    % Update tmpX
    tmpX = tmpRes.X;
    fprintf('.');
end
fprintf('\n');

% Calculate beta
retRes.beta = retRes.varpi*retRes.vartheta';

% Update mu
retRes.mu = retRes.mu-mX*retRes.beta;

retSIdx = find(idcS);
retX = inX(:,[inEIdx,retSIdx]);



%--------------------------------------------------------------------------
%
function retRes = gPOCREScreenCox1POC(inY,inCens,inRC,inX,inPreRes,inIdc,inNX,inZeta,inOption)
%function retRes = gPOCREScreenCox1POC(inY,inCens,inRC,inX,inPreRes,inIdc,inNX,inZeta,inOption)
% 
% Construct one orthogonal component for Cox Proportional Hazards Model in gPOCRE.
%
% Model: log(h(t)/h0(t)) = X*beta = \sum_j (X*varpi_j)*vartheta_j'
%                        = \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: observed survival data;
%   inCens: 0 if inY is not censored, 1 if censored;
%   inRC: ranks of observed failture times with ties the same rank;
%   inX: X;
%   inPreRes: mu,W,omega,varpi,poc(j).vartheta(varthea of j-th component).
%   inIdc: a vector of 0 & 1 indicating whether the variables have been
%          selected in prevously constructed components;
%   inNX: number of variables to be selected for current component;
%   inZeta: zeta for bias correction;
%   inOption: 
%       tol -- tolerance,
%       maxit -- maximum number of iterations
%       maxvar: maximum number of selected variables.
%
% Output:
%   retRes: omega, [vartheta_1,...], X, varpi, theta, idcS.
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(9, 9, nargin));
[n,p] = size(inX);    % X with n objects p variables;
k = size(inY,2);    % Y with n objects k variables;
if( k>1 )
    error(' Multiple responses are not supported!');
end
if( n~=size(inY,1) )
    error('gPOCRE: X and Y should have the same number of rows!');
end
idxCmp = inPreRes.nCmp+1; % Index of current component

% Initial Values:
retRes.mu = inPreRes.mu;
retRes.vartheta = [inPreRes.vartheta zeros(k,1)];
retRes.omega = zeros(p,1);
retRes.varpi = zeros(p,1);
retRes.idcS = zeros(p,1);

tmpDen = ones(idxCmp,1);
for j=1:(idxCmp-1)
    tmpDen(j) = inPreRes.xomega(:,j)'...
                *bsxfun(@times,inPreRes.xomega(:,j),inPreRes.W);
end

omegaPre = retRes.omega;    omegaDiff = 1;  idIter = 0;
while( (omegaDiff>inOption.tol) && (idIter<inOption.maxit) )
    idIter = idIter+1;
    omegaPre = retRes.omega;

    % Update eta
    %eta = retRes.mu+[inPreRes.xomega inX*retRes.omega]*retRes.vartheta';
    eta = [inPreRes.xomega inX*retRes.omega]*retRes.vartheta';
    
    % Update Z: Z = eta + [Score Function]/[Empirical Fisher Information];
    Z = eta + CalcCoxRatio(eta,inCens,inRC)*(sqrt(5)-1)/2;

    % Robust values of Z
    Z = robustize(Z);

    % Weighted Z
    Z = bsxfun(@times,Z,inPreRes.W);

    % Update mu
    retRes.mu = sum(Z,1);   % sum(inPreRew.W) = 1;
    disp([ 'mu = ' num2str(retRes.mu)]);
    
    % Update omega
    retRes.omega = gPOCREScreenSLEV(Z,inX,inIdc,inNX,inOption);
    
    if( norm(retRes.omega)<eps )
        omegaDiff = 0;
        retRes.vartheta(:,idxCmp) = zeros(k,1);
    else
        % Update all vartheta
        for j=1:(idxCmp-1)
            retRes.vartheta(:,j) = Z'*inPreRes.xomega(:,j)/tmpDen(j);
        end
        tmpDen(idxCmp) = retRes.omega'*inX'...
                         *bsxfun(@times,inX,inPreRes.W)*retRes.omega;
        if( tmpDen(idxCmp) < eps )
            retRes.vartheta(:,idxCmp) = zeros(k,1);
        else
            retRes.vartheta(:,idxCmp) = Z'*inX*retRes.omega/tmpDen(idxCmp);
        end
    
        omegaDiff = norm(omegaPre-retRes.omega);
    end
end

% Update X
retRes.theta = inX'*bsxfun(@times,inX,inPreRes.W)*retRes.omega...
               /tmpDen(idxCmp);     % px1
retRes.X = inX-inX*retRes.omega*retRes.theta';

% Calculate varpi
retRes.varpi = retRes.omega;
for j=(idxCmp-1):-1:1
    retRes.varpi = retRes.varpi-inPreRes.omega(:,j)...
                   *(inPreRes.theta(:,j)'*retRes.varpi);
end

% Update retRes.idcS
retRes.idcS = max((abs(retRes.omega)>eps)-inIdc,0);



%--------------------------------------------------------------------------
%
function retRes = gPOCREScreen1POC(inY,inX,inFamily,inLink,inPreRes,inIdc,inNX,inZeta,inOption)
%function retRes = gPOCRE1POC(inY,inX,inFamily,inLink,inPreRes,inIdc,inNX,inZeta,inOption)
% 
% Construct one orthogonal component in gPOCRE. Please refer details in
%
% Lin, Y., Zhang, M., and Zhang, D. (2011) Penalized orthogonal-component
% regression for fitting high-dimensional generalized linear models. Purdue
% Statistics Technical Report ???.
%
% Model: g(E[Y]) = mu + X*beta = mu + \sum_j (X*varpi_j)*vartheta_j'
%                              = mu + \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: Y response values of original data;
%   inX: X;
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   inPreRes: mu,W,omega,varpi,poc(j).vartheta(varthea of j-th component).
%   inIdc: a vector of 0 & 1 indicating whether the variables have been
%          selected in prevously constructed components;
%   inNX: number of variables to be selected for current component;
%   inZeta: zeta for bias correction;
%   inOption: 
%       tol -- tolerance,
%       maxit -- maximum number of iterations
%       maxvar: maximum number of selected variables;
%       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:
%   retRes: mu, omega, [vartheta_1,...], X, varpi, theta, idcS.
%
% Dabao Zhang    July 27, 2016
%
% Revised by:
%       Dabao Zhang, Oct. 3, 2018 -- Add inFamily.
%       Yan Li, Feb. 18, 2025 -- the number of input: error(nargchk(8, 8,
%       nargin))--error(nargchk(8, 9, nargin))
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(8, 9, nargin));
[n,p] = size(inX);    % X with n objects p variables;
k = size(inY,2);    % Y with n objects k variables;
if( n~=size(inY,1) )
    error('gPOCRE: X and Y should have the same number of rows!');
end
idxCmp = inPreRes.nCmp+1; % Index of current component

bLatent = 0;    % Latent Variable Model: multiple responses share identical X*beta.
if( sum(strcmp(inLink,{'pologit','crlogit'})) )
    bLatent = 1;
end

% Initial Values:
retRes.mu = inPreRes.mu;
if( bLatent )
    retRes.vartheta = [inPreRes.vartheta zeros(1,1)];
else
    retRes.vartheta = [inPreRes.vartheta zeros(k,1)];
end
retRes.omega = zeros(p,1);
retRes.varpi = zeros(p,1);
retRes.idcS = zeros(p,1);

tmpDen = ones(idxCmp,1);
for j=1:(idxCmp-1)
    tmpDen(j) = inPreRes.xomega(:,j)'...
                *bsxfun(@times,inPreRes.xomega(:,j),inPreRes.W);
end

omegaPre = retRes.omega;    omegaDiff = 1;  idIter = 0;
while( (omegaDiff>inOption.tol) && (idIter<inOption.maxit) )
    idIter = idIter+1;
    omegaPre = retRes.omega;

    % Update eta
    %eta = retRes.mu+[inPreRes.xomega inX*retRes.omega]*retRes.vartheta';
    eta = bsxfun(@plus,[inPreRes.xomega inX*retRes.omega]...
                 *retRes.vartheta',retRes.mu);
    
    % Update Z
    if( k==1 )  
        Z = eta + (inY+inZeta/2-(1+inZeta).*ginv(eta,inLink))...
            ./((1+inZeta).*dginv(eta,inLink));

        % Robust values of Z
        Z = robustize(Z);
    %elseif( bLatent ) % How to correct the bias???
        %Z = arrayfun(@(i)(eta(i,:)+(inY(i,:)-ginv(eta(i,:),inLink))...
        %            /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        
        %Z = cell2mat(Z);
    else
        %Z = arrayfun(@(i)(eta(i,:)+(inY(i,:)+inZeta(i)/2 ...
        %           -(.5*(k+1)*inZeta(i)+1).*ginv(eta(i,:),inLink)) ...
        %           /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        Z = arrayfun(@(i)(eta(i,:)+((inY(i,:)+inZeta(i)/2)./ ...
                   (1+.5*(k+1)*inZeta(i))-ginv(eta(i,:),inLink)) ...
                   /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        %Z = arrayfun(@(i)(eta(i,:)+((inY(i,:)+inZeta(i)/(k+1))./ ...
        %           (1+inZeta(i))-ginv(eta(i,:),inLink)) ...
        %           /dginv(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        Z = cell2mat(Z);
        
        % Robust values of Z
        Z = arrayfun(@(j)(robustize(Z(:,j))),1:k,'UniformOutput',false);
        Z = cell2mat(Z);
    end

    % Weighted Z
    Z = bsxfun(@times,Z,inPreRes.W);

    % Update mu
    retRes.mu = sum(Z,1);   % sum(inPreRew.W) = 1;
    %disp([ 'mu = ' num2str(retRes.mu)]);
    
    % Update omega
    if( bLatent )
        %tmpv = arrayfun(@(i)(CalcCatVar(eta(i,:),inLink)),[1:n]','UniformOutput',false);
        %optw = sum(cell2mat(tmpv),1);
        %optw = optw/sum(optw);
        %%disp([' category-wise weights: ' num2str(optw)]);
        
        % Any better way to handle here?
        optw = 1./sum(CalcIRLSResponseVar(eta,inFamily,inLink),1);
        optw = optw/sum(optw);

        optZ = Z*optw';
        colZ = 1;
    else
        optZ = Z;
        colZ = k;
    end
    retRes.omega = gPOCREScreenSLEV(optZ,inX,inIdc,inNX,inOption);
    
    %if( norm(retRes.omega)<inOption.tol )
    if( norm(retRes.omega)<eps )
        omegaDiff = 0;
        %retRes.vartheta(:,idxCmp) = zeros(k,1);
        retRes.vartheta(:,idxCmp) = zeros(colZ,1);
    else
        % Update all vartheta
        for j=1:(idxCmp-1)
            %retRes.vartheta(:,j) = Z'*inPreRes.xomega(:,j)/tmpDen(j);
            retRes.vartheta(:,j) = optZ'*inPreRes.xomega(:,j)/tmpDen(j);
        end
        tmpDen(idxCmp) = retRes.omega'*inX'...
                         *bsxfun(@times,inX,inPreRes.W)*retRes.omega;
        %if( tmpDen(idxCmp) < inOption.tol )
        if( tmpDen(idxCmp) < eps )
            %retRes.vartheta(:,idxCmp) = zeros(k,1);
            retRes.vartheta(:,idxCmp) = zeros(colZ,1);
        else
            %retRes.vartheta(:,idxCmp) = Z'*inX*retRes.omega/tmpDen(idxCmp);
            retRes.vartheta(:,idxCmp) = optZ'*inX*retRes.omega/tmpDen(idxCmp);
        end
    
        omegaDiff = norm(omegaPre-retRes.omega);
    end
end

% Update X
retRes.theta = inX'*bsxfun(@times,inX,inPreRes.W)*retRes.omega...
               /tmpDen(idxCmp);     % px1
retRes.X = inX-inX*retRes.omega*retRes.theta';

% Calculate varpi
retRes.varpi = retRes.omega;
for j=(idxCmp-1):-1:1
    retRes.varpi = retRes.varpi-inPreRes.omega(:,j)...
                   *(inPreRes.theta(:,j)'*retRes.varpi);
end

% Update retRes.idcS
retRes.idcS = max((abs(retRes.omega)>eps)-inIdc,0);



%-------------------------------------------------------------------------
%
function retLEV = gPOCREScreenSLEV(inZ,inX,inIdc,inNX,inOption)
%function retLEV = gPOCREScreenSLEV(inZ,inX,inIdc,inOption)
% 
% Calculate the leading sparse eigenvecor of X'ZZ'X.
%
% Input:
%   inZ: Z, which can be multivariate
%   inX: X
%   inIdc: a vector of 0 & 1 indicating whether the variables have been
%          selected in prevously constructed components;
%   inNX: number of variables to be selected for current component;
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%
% Output:
%   retLEV: the leading sparse eigenvector of X'ZZ'X.
%
% Dabao Zhang     July 27, 2016
%
% Revised by:
%
% Copyright (c) 2016 by Dabao Zhang.

error(nargchk(5, 5, 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);
    
    % Q: Hard Threshold or Soft Threshold?
    tmpCutoff = kthlargestvalue(abs(omega(~inIdc)),inNX+1);
    omega = max(abs(omega)-tmpCutoff,0).*sign(omega);
            
    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
