
function [pval chi2]=manteltwo(tA, cA, tB, cB, alpha)
%--------------------------
% Mantel's LogRank Test (1966)  
% Represents a generalization of 
% Savage (1956) Two-sample Rank Test
% to censored samples. 
%
% [pval chi2]=manteltwo(tA, cA, tB, cB, alpha)
%  Output: pval, chi2 -  p-value logrank test and statistic
%  Input:tA,  cA, tB, cB -- obseved times for group  A, followed by
%        censooring information; the same for group B.
%  Code in cA, cB:   0 - observed, 1- censored
%
%
% Example of use:
% tA =[1 3 3 6 7   7 10 12 14 15 ...
%          18 19 22 26 28    29 34 40 48 49];
% cA=[0 0 0 0 0   0 0 0 0 0  ...
%          0 0 0 0 1   0 0 0 1 1];
% tB =[1 1 2 2 3    4 5 8 8 9   ...
%     11 12 14 16 18    21 27 31 38 44]
% cB=[0 0 0 0 0    0 0 0 0 0   ...
%          0 0 0 0 0   0 1  0 1 0];
%      [pval chi2]=manteltwo(tA,  cA, tB,  cB)
%--------------------
     if nargin < 5 
    alpha=0.05;
end
if nargin < 3
    cA = zeros(size(tA));
    cB = zeros(size(tB));
end
    tA=tA(:)'; tB=tB(:)'; cA=cA(:)'; cB=cB(:)'; 
     tAo = tA( cA == 0 );
     tBo = tB( cB == 0 );
     t=[tAo tBo];
     from=[1*ones(1,length(tAo))  2*ones(1,length(tBo))];
     [st in] = sort(t);
     froms=from(in);
     %
     tu = unique(t);
     r = length(tu);
     
     for i=1:r
         d(i)=0; s(i)=0;
         for j=1:length(t);
             if( t(j) == tu(i) ) 
                  d(i) = d(i)+1;
                  s(i)=s(i) + froms(j)-1;
             end
         end
             d2(i)=s(i);
             d1(i)=d(i)-d2(i);
     end       
         % 
        n1(1)=length(tA);
        n2(1)=length(tB);
         for i=2:r
             n1(i) =  n1(1) -  sum( tA < tu(i));
             n2(i) =  n2(1) -  sum( tB < tu(i)); 
         end
         n = n1 + n2;       
         % 
         U0 = sum( d2 - d.*n2./n);
         % I0 = sum( d .* n1 .* n2 ./ n.^2 ); %is usual, 
         % but more precise is:  
         I0  = sum( d .* n2 ./ n   .* (1 - n2./n) .*  (n - d)./max(1,n-1)  );
               
         chi2 = U0^2/I0;         
         pval = 1-chi2cdf(chi2, 1);
end