function [retRes,retStat] = gPOCRE4Cox(inY,inDelta,inX,inTP,inOption)
%function [retRes,retStat] = gGOCRE4Cox(inY,inDelta,inX,inTP,inOption)
% 
% Penalized Orthogonal-Component Regression to fit generalized linear
% models (gPOCRE) with high dimensional survival data, following Cox's
% proportional hazard model. The number of non-zero coefficients is limited
% (under maxvar =n/2 by default).
%
% Model: pl(eta), with eta = X*beta = \sum_j (X*varpi_j)*vartheta_j'
%                              = \sum_j (X_j*omega_j)*vartheta_j';
%
% Input:
%   inY: Y response values of original data;
%   inDelta: observation status, 0 if censored, 1 if not censored;
%   inX: X;
%   inTP: tuning parameter;
%   inOption:
%       tol -- tolerance;
%       maxit -- maximum number of iterations;
%       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).
%
% Dabao Zhang    January 6, 2015
%
% Revised by:
%
% Copyright (c) 2015 by Dabao Zhang.

error(nargchk(3, 5, 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('    gPOCRE4Cox: X and Y should have the same number of rows!');
end

if( nargin<5 ) 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 = n/log(n); end
if( ~isfield(inOption,'W') ) inOption.W = []; end
if( ~isfield(inOption,'ptype') ) inOption.ptype = 'ebtz'; end
if( ~isfield(inOption,'covidx') ) inOption.covidx = []; end

% Only implemented for single response!
if( k>1 )
    error('    gPOCRE4Cox: only implemented for single time!');
end

% Initial values
retRes.nCmp = 0;
%retRes.mu = glnk(mean(inY),inLink);
retRes.beta = zeros(p,k);
retRes.vartheta = zeros(k,0);
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.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) )
    idxCmp=idxCmp+1;

    % Construct the component
    tmpRes = Build1SOC4Cox(inY,tmpX,inLink,retRes,inTP,inOption);
    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
if( nargout>1 )
    % Need to implement ...

end
