function retRes = gPOCRECox1POC(inY,inCens,inRC,inX,inPreRes,inTP,inZeta,inOption)
%function retRes = gPOCRECox1POC(inY,inCens,inRC,inX,inPreRes,inTP,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).
%   inTP: tuning parameter (lambda);
%   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: omega, [vartheta_1,...], X, varpi, theta;
%
% Dabao Zhang    May 23, 2013
%
% Revised by:
%       Dabao Zhang, July 26, 2016 -- Change the name.
%
% Copyright (c) 2013 by Dabao Zhang.

error(nargchk(8, 8, 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);

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 = CalcSLEV(Z,inX,inTP,inOption);
    
    %if( norm(retRes.omega)<inOption.tol )
    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) < inOption.tol )
        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

%if( norm(retRes.omega)>=inOption.tol )
if( norm(retRes.omega)>eps )
    % 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
else
    retRes.mu = inPreRes.mu;
end
