# HG changeset patch
# User Dave Goel
# Date 1466542264 -7200
# Tue Jun 21 22:51:04 2016 +0200
# Node ID 3725f4a01810c0464911828a9dd8ddb3faf0e3ac
# Parent da218a61ce4a2a8fe30b262bee871d5c1cf87afc
Improve speed of ranks.m when dim!=2, and give choice of tie-breaks (bug #36372)
* scripts/statistics/base/rank.m: Replace algorithm.
diff -r da218a61ce4a -r 3725f4a01810 scripts/statistics/base/ranks.m
--- a/scripts/statistics/base/ranks.m Tue Jun 21 13:08:25 2016 -0700
+++ b/scripts/statistics/base/ranks.m Tue Jun 21 22:51:04 2016 +0200
@@ -1,4 +1,5 @@
## Copyright (C) 1995-2015 Kurt Hornik
+## Copyright (C) 2012 Dave Goel
##
## This file is part of Octave.
##
@@ -17,24 +18,31 @@
## .
## -*- texinfo -*-
-## @deftypefn {} {} ranks (@var{x}, @var{dim})
-## Return the ranks of @var{x} along the first non-singleton dimension
-## adjusted for ties.
+## @deftypefn {} {} ranks (@var{x})
+## @deftypefnx {} {} ranks (@var{x}, @var{dim})
+## @deftypefnx {} {} ranks (@var{x}, @var{dim}, @var{rtype})
##
-## If the optional argument @var{dim} is given, operate along this dimension.
+## 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 (1, 2, 2, 4);
+## 2 or "modified" for modified competition ranking (1, 3, 3, 4);
+## 3 or "ordinal" for ordinal ranking (1, 2, 3, 4);
+## 4 or "dense" for dense ranking (1, 2, 2, 3).
## @seealso{spearman, kendall}
## @end deftypefn
-## Author: KH
-## Description: Compute ranks
+## Authors: Dave Goel , KH ,
+## Description: Compute ranks using a choice of methods when resolving ties.
-## This code was rather ugly, since it didn't use sort due to the
-## fact of how to deal with ties. Now it does use sort and its
-## even uglier! At least it handles NDArrays.
+function y = ranks (x, dim, rtype)
-function y = ranks (x, dim)
-
- if (nargin != 1 && nargin != 2)
+ if (nargin < 1) || (nargin > 3);
print_usage ();
endif
@@ -42,9 +50,14 @@
error ("ranks: X must be a numeric vector or matrix");
endif
+ if (nargin < 3) || ( isempty (rtype));
+ rtype = 0;
+ endif
+
nd = ndims (x);
sz = size (x);
- if (nargin != 2)
+
+ if (nargin < 2) || isempty (dim);
## Find the first non-singleton dimension.
(dim = find (sz > 1, 1)) || (dim = 1);
else
@@ -55,35 +68,79 @@
endif
if (sz(dim) == 1)
- y = ones (sz);
+ y = ones (sz); ## Dimension dim is singleton, so all are ranked first.
else
- ## The algorithm works only on dim = 1, so permute if necesary.
+ ## 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
- sz = size (x);
- infvec = -Inf ([1, sz(2 : end)]);
- [xs, xi] = sort (x);
- eq_el = find (diff ([xs; infvec]) == 0);
- if (isempty (eq_el))
- [eq_el, y] = sort (xi);
- else
- runs = setdiff (eq_el, eq_el+1);
- len = diff (find (diff ([Inf; eq_el; -Inf]) != 1)) + 1;
- [eq_el, y] = sort (xi);
- for i = 1 : length (runs)
- y (xi (runs (i) + [0:(len(i)-1)]) + floor (runs (i) ./ sz(1))
- * sz(1)) = eq_el(runs(i)) + (len(i) - 1) / 2;
- endfor
- endif
+
+ [sx ids] = sort (x); ## sx is sorted x.
+ 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, "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.
+ 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
@@ -96,7 +153,7 @@
## Test input validation
%!error ranks ()
-%!error ranks (1, 2, 3)
+%!error ranks (1, 2, 3, 4)
%!error ranks ({1, 2})
%!error ranks (['A'; 'B'])
%!error ranks (1, 1.5)