# HG changeset patch # User Lachlan Andrew # Date 1449185073 -39600 # Node ID 1b59c82e0332f6eb1ff347b2d35515652c2a1073 # Parent c06cbec3a88328a3bf531f431d830567f3d4dcf2 Fix bug #36372: Faster ranking of large matrices, and new tie-break options * scripts/statistics/base/ranks.m: New algorithms using cummax, not for loop diff -r c06cbec3a883 -r 1b59c82e0332 scripts/statistics/base/ranks.m --- a/scripts/statistics/base/ranks.m Sun Nov 29 19:58:52 2015 +1100 +++ b/scripts/statistics/base/ranks.m Fri Dec 04 10:24:33 2015 +1100 @@ -1,4 +1,5 @@ ## Copyright (C) 1995-2015 Kurt Hornik +## Copyright (C) 2012 Dave Goel ## ## This file is part of Octave. ## @@ -17,34 +18,44 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {} ranks (@var{x}, @var{dim}) -## Return the ranks of @var{x} along the first non-singleton dimension -## adjusted for ties. ## -## If the optional argument @var{dim} is given, operate along this dimension. +## @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. +## +## 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 "revordinal" for reverse ordinal (4, 3, 2, 1), and +## 5 or "dense" for dense ranking (1, 2, 2, 3). +## ## @seealso{spearman, kendall} ## @end deftypefn -## Author: KH +## Authors: KH , Dave Goel ## Description: Compute ranks -## 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 - if (! (isnumeric (x) || islogical (x))) - error ("ranks: X must be a numeric vector or matrix"); + if (nargin < 3 || isempty (rtype)) + rtype = 0; endif nd = ndims (x); sz = size (x); - if (nargin != 2) + if (nargin == 1 || isempty (dim)) ## Find the first non-singleton dimension. (dim = find (sz > 1, 1)) || (dim = 1); else @@ -57,38 +68,84 @@ function y = ranks (x, dim) if (sz(dim) == 1) y = ones (sz); 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); + + if (rtype == 4 || strcmp (rtype, "revordinal")) + [sx ids] = sort (x, "descend"); + ids = flipud (ids); 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 + [sx ids] = sort (x); ## sx is sorted x. endif + 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 + 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 + zeros ([1, sz(2 : end)]); + fnewp = logical (diff ([infvec;sx])); + linnew = cumsum (fnewp, 1); +endfunction -%!assert (ranks (1:2:10), 1:5) -%!assert (ranks (10:-2:1), 5:-1:1) +function linnew = _competition (lin, sx, sz) + ## Like _dense but stop increasing lin when sx does not increase. + infvec = -Inf + zeros ([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 + zeros ([1, sz(2 : end)]); + fnewp= find (diff ([sx;infvec])); + linnew = Inf (size (lin)); + linnew(fnewp) = lin(fnewp); + linnew = flipud (cummin (flipud (linnew))); +endfunction + +%!assert (ranks (1:2:10), [1:5]) +%!assert (ranks (10:-2:1), [5:-1:1]) %!assert (ranks ([2, 1, 2, 4]), [2.5, 1, 2.5, 4]) %!assert (ranks (ones (1, 5)), 3*ones (1, 5)) %!assert (ranks (1e6*ones (1, 5)), 3*ones (1, 5)) @@ -96,10 +153,9 @@ endfunction ## 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) %!error ranks (1, 0) %!error ranks (1, 3) -