function retRatio = CalcCoxRatio(inEta,inCens,inRC)
%function retRatior = CalcCoxRatio(inEta,inCens,inRC)
% 
% Calculate U*inv(I), with U as the score and I as the empirical Fisher
% information matrix, for the partial likelihood function following Cox
% proportional hazards model.
%
% Input:
%   inEta: eta;
%   inCens: 0 if Y is not censored, 1 if censored;
%   inRC: ranks of observed failture times with ties the same rank;
%
% Output:
%   retRatio: Score*inv(Empirical Fisher Information].
%
% Dabao Zhang    March 19, 2015
%
% Revised by:
%
% Copyright (c) 2015 by Dabao Zhang.

error(nargchk(3, 3, nargin));

% Compute the score: dL/dEta
bUnC = ~inCens;
expEta = exp(inEta);
%riskSum = flipud(cumsum(flipud(bUnC.*expEta)));
riskSum = flipud(cumsum(flipud(expEta)));
riskSum = riskSum(inRC(:,1));
invSum = cumsum(bUnC./riskSum);
invSum = invSum(inRC(:,2));
tmpSc = bUnC-expEta.*invSum;

% Compute the empirical Fisher Information: related to -d2L/dEta2
invSum2 = cumsum(bUnC./(riskSum.^2));
invSum2 = invSum2(inRC(:,2));
tmpI = expEta.*invSum-(expEta.^2).*invSum2; 

%idxMin = find(tmpI==min(tmpI)); idxMin = idxMin(1);
%idxMax = find(tmpI==max(tmpI)); idxMax = idxMax(1);

%disp([' min(Hessian): ' num2str(tmpI(idxMin))]);
%disp([' Score: ' num2str(tmpSc(idxMin))]);

%disp([' max(Hessian): ' num2str(tmpI(idxMax))]);
%disp([' Score: ' num2str(tmpSc(idxMax))]);

retRatio = tmpSc./tmpI;

%disp([' min(Ratio): ' num2str(min(retRatio))]);
%disp([' max(Ratio): ' num2str(max(retRatio))]);

retRatio(retRatio>10) = 10;
retRatio(retRatio<-10) = -10;
