## 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
## <http://www.gnu.org/licenses/>.

## -*- 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 <deego3@gmail.com>, KH <Kurt.Hornik@wu-wien.ac.at>, 
## Description: Compute ranks using choice of methods when resolving ties. 
## Thanks: Lachlan Andrew <lachanbis@gmail.com>


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)
