function [retRes,retStat] = gPOCRE(inY,inX,inFamily,inLink,inTP,inOption)
%function [retRes,retStat] = gGOCRE(inY,inX,inFamily,inLink,inTP,inOption)
% 
% Penalized Orthogonal-Component Regression to fit generalized linear
% models (gPOCRE) with high dimensional data. The deviance will be
% reported. The number of non-zero coefficients is limited (under maxvar
% =n/2 by default). Please refer gPOCRE to
%
% 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. 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(.);
%   inTP: tuning parameter;
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%       bctype -- type of bias correction: 0 ~ no bias correction,
%                 1 ~ approximate, type 2 -- optimal (default);
%       maxcmp -- maximum number of components;
%       maxvar -- maximum number of selected variables;
%       W -- prespecified weight vector for the weight matrix;
%       ptype -- type of penalization,
%                'ebt' ~ EBT w/o z-transformation,
%                'ebtz' ~ EBT w/ z-transformation (default),
%                '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);
%       covidx -- indices of covariates of no interest from components.
%
% Output:
%   retRes: nCmp, mu, beta, deviance, W, bSparse(is beta sparse), vartheta,
%           omega, varpi, xomega ([X_1*omega_1, ...] = X*varpi), theta,
%           poc(j).vartheta (vartheta of j-th component), poc(j).mu,
%           poc(j).beta.
%   retStat: some statistical/inference results based on the model with all
%            non-zero beta, e.g., pvalue, seqpv, indpv, loglik, effp
%            (effective number of predictors), separated.
%
% 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    May 28, 2013
%
% Revised by:
%   Dabao Zhang, 9/16/2013 -- Add inOption.covidx, retRes.pvalue,
%                             retRes.seqpv, retRes.indpv;
%   Dabao Zhang, 9/23/2013 -- Add retStat and calculate effp.
%   Dabao Zhang, 10/29/2013 -- Change the initialization of W.
%   Dabao Zhang, 10/31/2013 -- Add poc(j).beta.
%   Dabao Zhang, 1/19/2015 -- Multinomial data.
%   Dabao Zhang, 7/27/2016 -- Change function names.
%   Dabao Zhang, 12/19/2017 -- default maxvar.
%   Dabao Zhang, 12/19/2017 -- Construct all components for the first
%                              non-sparse tuning parameter.
%   Dabaqo Zhang, 7/2/2018 -- Error in reporting "perfect separation"?
%
% Copyright (c) 2013-2017 by Dabao Zhang.

error(nargchk(5, 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 subjects p variables;
k = size(inY,2);    % Y with n subjects 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 = 2; end
if( ~isfield(inOption,'maxcmp') ) inOption.maxcmp = 10; end
if( ~isfield(inOption,'maxvar') ) inOption.maxvar = min(n,p)/2; end
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';
        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.bSparse = 1;
retRes.n = n;
retRes.p = p;
retRes.id = 'gpocre';
retRes.poc = [];
retRes.omega = zeros(p,0);
retRes.theta = zeros(p,0);
if( size(retRes.W,2)~=1 ) retRes.W = retRes.W'; end

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

% Construct all components
idxCmp = 0;     idxNZBeta = [];
omeganrm = 1;
fprintf('\b  ');
zeta = CalcZeta4BC(tmpX,retRes.W,inOption.bctype);
%while( (idxCmp<inOption.maxcmp) && retRes.bSparse && (omeganrm>eps) )
while( (idxCmp<inOption.maxcmp) && (omeganrm>eps) )
    idxCmp=idxCmp+1;

    % Construct the component
    if( strcmp(inFamily,'cox') )
        tmpRes = gPOCRECox1POC(inY,inCens,inRC,tmpX,retRes,inTP,zeta,inOption);
    else
        tmpRes = gPOCRE1POC(inY,tmpX,inFamily,inLink,retRes,inTP,zeta,inOption);
    end
    omeganrm = norm(tmpRes.omega);
    idxNZBeta = union(idxNZBeta,find(abs(tmpRes.varpi)>eps));
    if( length(idxNZBeta)>inOption.maxvar )
        retRes.bSparse = 0;
    elseif( omeganrm>eps )        % Save the information on constructed component
        retRes.mu = tmpRes.mu;
        retRes.vartheta = tmpRes.vartheta;
        retRes.poc(idxCmp).mu = tmpRes.mu;
        retRes.poc(idxCmp).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
end
fprintf('\n');

% Calculate beta
retRes.beta = retRes.varpi*retRes.vartheta';
for j=1:retRes.nCmp
    retRes.poc(j).beta = retRes.varpi(:,1:j)*retRes.vartheta(:,1:j)';
end

% Update mu
retRes.mu = retRes.mu-mX*retRes.beta;
for j=1:retRes.nCmp
    retRes.poc(j).mu = retRes.poc(j).mu-mX*retRes.poc(j).beta;
end

% Calculate p-values and log-likelihood value if requested
% Currently only implemented for univariate response
% See mnrfit() to implement likelihood ratio tests
if( nargout>1 )
    % Special for multinomial models
    if( sum(strcmp(inFamily,{'multibc','multipo','multicr'})) )
        inY = [1-sum(inY,2) inY];
    end
    
    tmpX = bsxfun(@minus,inX,mX);
    if( strcmp(inFamily,'cox') )
        retStat = EvalCoxCompsSig(inY,inCens,tmpX,retRes,inOption);
    else
        retStat = EvalCompsSig(inY,tmpX,inFamily,inLink,retRes,inOption);
    end
    
    % Calculate the deviance for SelectModel() & Report 'PerfectSeparation'
    origstat = warning('off','all');
    warning('on','stats:glmfit:PerfectSeparation');
    msgid = '';         lastwarn('');

    %tmpX = inX(:,abs(retRes.beta)>eps);
    tmpX = inX(:,sum(abs(retRes.beta),2)>eps);
    [tmpR,tmpJ] = rref(tmpX);
    retStat.effp = length(tmpJ)+1;
    tmpX = tmpX(:,tmpJ);

    %[msgstr,msgid] = lastwarn;
    switch inLink
        case 'bclogit'
            [tb,retStat.deviance] = mnrfit(tmpX,inY,'model','nominal');
            %retStat.separated = 0;

        case 'pologit'
            [tb,retStat.deviance] = mnrfit(tmpX,inY,'model','ordinal','link','logit');
            %retStat.separated = 0;

        case 'crlogit'
            [tb,retStat.deviance] = mnrfit(tmpX,inY,'model','hierarchical',...
                                'interactions','off');
            %retStat.separated = 0;

        case 'cox' % -2*logL here!
            [tb,retStat.deviance,tH,ts ] = coxphfit(tmpX,inY,'censoring',inCens);
            retStat.deviance = -2*retStat.deviance;
            %retStat.separated = 0;  % Need to work out!
            
        otherwise  % univariate response
            [tb,retStat.deviance] = glmfit(tmpX,inY,inFamily,'link',inLink);
            %retStat.separated = strcmp(msgid,'stats:glmfit:PerfectSeparation');
    end

    [msgstr,msgid] = lastwarn;
    retStat.separated = strcmp(msgid,'stats:glmfit:PerfectSeparation');
    warning(origstat);
else
    retStat = [];
end




%-------------------------------------------------------------------------
%
function retZeta = CalcZeta4BC(inX,inW,inType)
%function retZeta = CalcZeta4BC(inX,inW,inType)
% 
% Calculate the diagonal elements of W^{1/2}X(X'WX)^{-1}X'W^{1/2} for bias
% correction.
%
% Input:
%   inType: 0 -- no bias correction, i.e., zeta=0_{nxn};
%           1 -- approximate values, i.e., zeta=(1-1/n)I_{nxn};
%           2 -- exact values
%
% Output:
%   retZeta: a vector.
%
% Dabao Zhang    May 23, 2013
%
% Revised by:
% add line 323 by Yan Li --Feb 18,2025
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(3, 3, nargin));
[n,p] = size(inX);
rowX = n
retZeta = zeros(n,1);
switch( inType )
    %case 0
        %retZeta = retZeta + 0.25;
    case 1
        retZeta = 1-inW;    % sum(inW)=1.
    case 2
        % X=USV': nxp; X(X'X)^{-1}X'=UU'
        tmpM = bsxfun(@times,[ones(rowX,1) inX],sqrt(inW));
        [tmpU,noU1,noU2] = svd(tmpM,'econ');
        retZeta = diag(tmpU*tmpU');
end



%--------------------------------------------------------------------------
%
function retStat = EvalCompsSig(inY,inX,inFamily,inLink,inRes,inOption)
%function retStat = EvalCompsSig(inY,inX,inFamily,inLink,inRes,inOption)
%
% Evaluate the significance of constructed components.
%
% Note: This function is implemented only for univariate Y.
%
% Input:
%   inY: Y --- allow for multiple response variables
%   inX: X
%   inFamily: the distribution family;
%   inLink: the link function g(.);
%   inRes: first output of gPOCRE;
%   inOption:
%       covidx -- indices of covariates of no interest from components.
%
% Dabao Zhang     November 6, 2013
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

switch inFamily
    case {'normal','binomial','poisson','gamma','inverse gaussian'}
        bMultiN = 0;
    case {'multibc','multipo','multicr'}
        k = size(inY,2)-1;
        bMultiN = 1;
    case 'quasi'
        error('gPOCRE: quasi is not supported...');
    otherwise
        error('gPOCRE: The family is not supported...');
end
%tmpX = bsxfun(@minus,inX,mX);

% Covariates of no interest
covx = [];
for j = 1:inRes.nCmp
    if( sum(abs(inRes.varpi(inOption.covidx,j)))>eps )
        covx = [covx ...
          inX(:,inOption.covidx)*inRes.varpi(inOption.covidx,j)];
    end
end
ncovx = size(covx,2);

% Calculate the p-values for each component:
tidx = setdiff(1:size(inX,2),inOption.covidx);
tmpX = [covx inX(:,tidx)*inRes.varpi(tidx,:)];
switch inLink
    case 'bclogit'
        [tb,td,ts] = mnrfit(tmpX,inY,'model','nominal');
        retStat.pvalue = ts.p((ncovx+2):end,:);

        retStat.seqpv = retStat.pvalue;
        for j=1:(inRes.nCmp-1)
            [tb,td,ts] = mnrfit(tmpX(:,1:(ncovx+j)),inY);
            retStat.seqpv(j,:) = ts.p(j+ncovx+1,:);
        end

        retStat.indpv = retStat.seqpv;
        for j=2:inRes.nCmp
            [tb,td,ts] = mnrfit(tmpX(:,[1:ncovx j+ncovx]),inY);
            retStat.indpv(j,:) = ts.p(ncovx+2,:);
        end
    
    case 'pologit'
        [tb,td,ts] = mnrfit(tmpX,inY,'model','ordinal','link','logit');
        retStat.pvalue = ts.p((ncovx+k+1):end);

        retStat.seqpv = retStat.pvalue;
        for j=1:(inRes.nCmp-1)
            [tb,td,ts] = mnrfit(tmpX(:,1:(ncovx+j)),inY,'model','ordinal');
            retStat.seqpv(j) = ts.p(j+ncovx+k);
        end

        retStat.indpv = retStat.seqpv;
        for j=2:inRes.nCmp
            [tb,td,ts] = mnrfit(tmpX(:,[1:ncovx j+ncovx]),inY,'model','ordinal');
            retStat.indpv(j) = ts.p(ncovx+k+1);
        end

    case 'crlogit'
         [tb,td,ts] = mnrfit(tmpX,inY,'model','hierarchical','interactions','off');
        retStat.pvalue = ts.p((ncovx+k+1):end);

        retStat.seqpv = retStat.pvalue;
        for j=1:(inRes.nCmp-1)
            [tb,td,ts] = mnrfit(tmpX(:,1:(ncovx+j)),inY,'model','hierarchical',...
                'interactions','off');
           retStat.seqpv(j) = ts.p(j+ncovx+k);
           %retStat.seqpv(j) = ts.p(j+ncovx+k-1);
        end

        retStat.indpv = retStat.seqpv;
        for j=2:inRes.nCmp
            [tb,td,ts] = mnrfit(tmpX(:,[1:ncovx j+ncovx]),inY,'model','hierarchical',...
                'interactions','off');
            retStat.indpv(j) = ts.p(ncovx+k+1);
            %retStat.indpv(j) = ts.p(ncovx+k);
        end
        
    otherwise  % univariate response
        [tb,td,ts] = glmfit(tmpX,inY,inFamily,'link',inLink);
        retStat.pvalue = ts.p((ncovx+2):end);

        retStat.seqpv = retStat.pvalue;
        for j=1:(inRes.nCmp-1)
            [tb,td,ts] = glmfit(tmpX(:,1:(ncovx+j)),inY,inFamily,'link',inLink);
            retStat.seqpv(j) = ts.p(j+ncovx+1);
        end

        retStat.indpv = retStat.seqpv;
        for j=2:inRes.nCmp
            [tb,td,ts] = glmfit(tmpX(:,[1:ncovx j+ncovx]),inY,inFamily,'link',inLink);
            retStat.indpv(j) = ts.p(ncovx+2);
        end
end



%--------------------------------------------------------------------------
%
function retStat = EvalCoxCompsSig(inY,inCens,inX,inRes,inOption)
%function retStat = EvalCoxCompsSig(inY,inCens,inX,inRes,inOption)
%
% Evaluate the significance of constructed components.
%
% Note: This function is implemented only for univariate Y.
%
% Input:
%   inY: Y --- allow for multiple response variables
%   inX: X
%   inCens: censoring status, 1 if censored, 0 otherwise;
%   inRes: first output of gPOCRE;
%   inOption:
%       covidx -- indices of covariates of no interest from components.
%
% Dabao Zhang
%
% Revised by:
%
% Copyright (c) 2013 by Dabao Zhang.

% Covariates of no interest
covx = [];
for j = 1:inRes.nCmp
    if( sum(abs(inRes.varpi(inOption.covidx,j)))>eps )
        covx = [covx ...
          inX(:,inOption.covidx)*inRes.varpi(inOption.covidx,j)];
    end
end
ncovx = size(covx,2);

% Calculate the p-values for each component:
tidx = setdiff(1:size(inX,2),inOption.covidx);
tmpX = [covx inX(:,tidx)*inRes.varpi(tidx,:)];

[tb,logL,tH,ts] = coxphfit(tmpX,inY,'censoring',inCens);
retStat.pvalue = ts.p((ncovx+1):end);

retStat.seqpv = retStat.pvalue;
for j=1:(inRes.nCmp-1)
    [tb,logL,tH,ts] = coxphfit(tmpX(:,1:(ncovx+j)),inY,'censoring',inCens);
    retStat.seqpv(j) = ts.p(j+ncovx);
end

retStat.indpv = retStat.seqpv;
for j=2:inRes.nCmp
    [tb,logL,tH,ts] = coxphfit(tmpX(:,[1:ncovx j+ncovx]),inY,'censoring',inCens);
    retStat.indpv(j) = ts.p(ncovx+1);
end
