## Copyright (C) 1995-2012 Kurt Hornik; Copyright (C) 2012 Dave Goel
##
## This file is (likely to be) part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
## for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
## -*- texinfo -*-
## @deftypefn {Function File} {} ranks (@var{x})
## @deftypefnx {Function File} {} ranks (@var{x}, @var{dim})
## @deftypefnx {Function File} {} ranks (@var{x}, @var{dim}, @var{rtype})
##
## Return the ranks (in the sense of order statistics) of @var{x} along
## the first non-singleton dimension
## adjusted for ties. If the argument @var{dim} is given and non-empty,
## operate along this dimension instead.
##
## The parameter @var{rtype} determines the way ties are handled:
## 0 or "fractional" (the default) for fractional, which gives ranks
## like (1, 2.5, 2.5, 4);
## 1 or "competition" for competition ranking (Example: [1, 2, 2, 4], @var{x} is [11 12 12 14];
## 2 or "modified" for modified competition ranking (1, 3, 3, 4);
## 3 or "ordinal" for ordinal ranking (1, 2, 3, 4);
## 4 or "revordinal" for reverse ordinal (1, 3, 2, 4)), and
## 5 or "dense" for dense ranking (1, 2, 2, 3).
## @seealso{spearman, kendall}
## @end deftypefn
## Authors: Dave Goel , KH ,
## Description: Compute ranks using choice of methods when resolving ties.
## Thanks: Lachlan Andrew
function y = ranks (x, dim, rtype)
if (nargin < 1) || (nargin > 3);
print_usage ();
endif
if (nargin < 3) || ( isempty(rtype));
rtype = 0;
endif
nd = ndims (x);
sz = size (x);
if (nargin < 2) || isempty(dim);
(dim = find (sz > 1, 1)) || (dim = 1);
endif
if (sz(dim) == 1)
y = ones(sz); ## x only has element along the dimension. Thus, all ranks equal 1.
else
## The algorithm works only on dim = 1, so permute if necessary.
if (dim != 1)
perm = [1 : nd];
perm(1) = dim;
perm(dim) = 1;
x = permute (x, perm);
sz=size(x);
endif
switch rtype;
case {4,"revordinal"};
[sx ids]=sort(x,'descend');
ids=flipdim(ids,1);
otherwise
[sx ids]=sort(x); ## sx is sorted x.
endswitch
lin=cumsum(ones(size(x)),1); ## A linearly increasing array.
switch rtype;
case {0,"fractional"};
lin=(_competition(lin,sx,sz)+_modified(lin,sx,sz))/2;
case {1,"competition"};
lin=_competition(lin,sx,sz);
case {2,"modified"};
lin=_modified(lin,sx,sz);
case {3,"ordinal"};
## no processing needed here.
case {4,"revordinal"};
## no processing needed here.
case {5,"dense"};
lin=_dense(lin,sx,sz);
otherwise
if (!ischar (rtype))
rtype = num2str (rtype);
end
error ("ranks: Illegal RTYPE '%s'", rtype);
endswitch
y=NaN(size(lin));
## Offsets to map indices into each column to indices into the linear array. (Thanks, Lachlan Andrew.)
idf = zeros (sz);
idf(1,:) = 0:sz(1): (numel (ids)-1);
idf(:,:) = repmat (idf (1,:), [sz(1), ones(1,length(sz)-1)]);
y(ids + idf) = lin;
if (dim != 1)
y = permute (y, perm);
endif
endif
endfunction
function linnew=_dense (lin,sx,sz)
infvec = -Inf * ones ([1, sz(2 : end)]);
fnewp= logical(diff([infvec;sx]));
linnew=cumsum(fnewp,1);
endfunction
function linnew=_competition (lin,sx,sz)
## Stop increasing lin when sx does not increase. Same as before
## otherwise.
infvec = -Inf * ones ([1, sz(2 : end)]);
fnewp= find(diff([infvec;sx]));
linnew=zeros(size(lin));
linnew(fnewp)=lin(fnewp);
linnew=cummax(linnew,1);
endfunction
function linnew=_modified (lin,sx,sz)
## Traverse lin backwards. Stop decreasing it when sx doesn't
## decrease.
infvec = Inf * ones ([1, sz(2 : end)]);
fnewp= find(diff([sx;infvec]));
linnew=Inf(size(lin));
linnew(fnewp)=lin(fnewp);
linnew=flipdim(cummin(flipdim(linnew,1)),1);
endfunction
%!test assertmy (ranks (1:2:10), [1:5]'');
%!test assertmy (ranks (10:-2:1), [5:-1:1]'');
%!test assertmy (ranks ([2, 1, 2, 4]), [2.5, 1, 2.5, 4]);
%!test assertmy (ranks (ones (1, 5)), 3*ones (1, 5)'');
%!test assertmy (ranks (1e6*ones (1, 5)), 3*ones (1, 5)'');
%!test assertmy (ranks (rand (1, 5), 1), ones (1, 5));
%% Test input validation
%!error ranks ()
%!error ranks (1, 2, 3, 4)
%!error ranks ({1, 2})
%!error ranks (['A'; 'B'])
%!error ranks (1, 1.5)
%!error ranks (1, 0)
%!error ranks (1, 3)