function [ x_rk ] = tiedrankrelative(x, idxRankable)
%TIEDRANKRELATIVE	Performs tied ranking for all elements, relative to a subset of the elements
%	[ x_rk ] = tiedrankrelative(x, idxRankable) ranks all of the elements
%	in the vector x relative to the elements indexed by idxRankable. Thus,
%	elements not indexed by idxRankable are given a rank, but do not
%	influence the ranking of points. This is useful if you would, for
%	example, like to rank points from an unused test set relative to a
%	training set.
%   
%   If idxRankable is not supplied, then this function uses all elements to
%   calculate their ranking, averaging the ranks of tied elements.
%
%   Note: this function may not work for vectors with Inf elements.
%   
%   If x is a matrix, the function operates column-wise.
%
%	Inputs:
%		x           - The input data elements.
%       idxRankable - Elements to be used in the ranking (i.e.,
%           idxTraining). Must be a logical vector of indices.
%		
%
%	Outputs:
%		x_rk        - Points in x ranked according to the elements indexed
%           by idxRankable.
%		
%
%	Example
%       load PhysionetDataSetA;
%       xtrain=data(1:500,:);
%       ytrain=outcome(1:500);
%       xtest=data(501:1000,:);
%       ytest=outcome(501:1000);
%       
%       idxTraining = [true(size(xtrain));...
%           false(size(xtest));
%		[ x_rk ] = tiedrankrelative([xtrain;xtest], idxTraining);
%	
%	See also TIEDRANK

%	Copyright 2012 Alistair Johnson

%	$LastChangedBy: alistair $
%	$LastChangedDate: 2012-05-15 12:26:18 +0100 (Tue, 15 May 2012) $
%	$Revision: 1 $
%	Originally written on GLNXA64 by Alistair Johnson, 30-Apr-2012 10:30:47
%	Contact: alistairewj@gmail.com


[N,P] = size(x);

if nargin<2
    idxTest = false(size(x));
else
    idxTest = ~idxRankable;
end

%=== Sort data
[x_sorted, idxRk] = sort(x,1,'ascend');

%=== Create default rankings 1:N before ties
x_nan = isnan(x_sorted);
x_rk = repmat((1:N)',1,P);
x_rk(x_nan) = NaN;

%=== Get logical indices of ties
idxTied = x_sorted;
idxTied = ((x_sorted(2:end,:)-idxTied(1:end-1,:)) == 0);


%=== Sort test indices and create adjustment vectors
idxTestRk = idxTest(idxRk);
adj1 = double(idxTestRk); % adjustment for -0.5 (observations == test)
adj2 = double(idxTestRk); % adjustment for -1 (observations > test)


for c=1:P
    %=== Get column specific numerical indices of ties
    idxTiedCol = [find(idxTied(:,c)==1);N+2];
    N_Tied = numel(idxTiedCol);
    i=1;
    while i < N_Tied
        idxTieStart = idxTiedCol(i);
        nties = 2;
        %=== Count number of ties
        while (idxTieStart+nties-1 == idxTiedCol(i+1))
            nties = nties+1;
            i=i+1;
        end
        
        %=== Take and impute average
        idxTieEnd = idxTiedCol(i)+1;
        idxTieAdjust = idxTieStart:idxTieEnd;
        x_rk(idxTieAdjust,c) = sum(x_rk(idxTieAdjust,c)) / nties;
        
        ntesttied = sum(idxTestRk(idxTieAdjust,c));
        idxTrainTestTied = ~idxTestRk(idxTieAdjust,c);
        if any(idxTrainTestTied) && ntesttied>0
            %=== Add -0.5 for training points tied with test points
            adj1(idxTieAdjust(idxTrainTestTied),c) = adj1(idxTieAdjust(idxTrainTestTied),c) + 1;
        end
        
        %=== If multiple tied points are from test set
        if ntesttied>0
            adj2(idxTieAdjust,c) = 0;
            adj2(idxTieAdjust(end),c) = ntesttied; % number of tied test points
            
            if ntesttied>1
                % must add more -0.5s to adj1
                adj1(idxTieAdjust,c) = adj1(idxTieAdjust,c) + ntesttied-1;
            end
        end
        i=i+1;
    end
end

%=== Shift adj2 down and do cumulative sum
adj2 = cumsum([zeros(1,size(adj2,2));adj2(1:end-1,:)],1);

%=== Create linear indices
idxLinearRk = ones(N,P); idxLinearRk = cumsum(idxLinearRk,2);
idxLinearRk = sub2ind([N,P],idxRk,idxLinearRk);
idxLinearRk = idxLinearRk(:);

x_rk(idxLinearRk) = x_rk - adj1*0.5 - adj2;
end