--- a/scripts/statistics/base/ranks.m	Fri Apr 30 04:12:37 2004 +0000
+++ b/scripts/statistics/base/ranks.m	Fri Apr 30 04:21:33 2004 +0000
@@ -18,37 +18,78 @@
 ## 02111-1307, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} ranks (@var{x})
+## @deftypefn {Function File} {} ranks (@var{x}, @var{dim})
 ## If @var{x} is a vector, return the (column) vector of ranks of
 ## @var{x} adjusted for ties.
 ##
-## If @var{x} is a matrix, do the above for each column of @var{x}.
+## If @var{x} is a matrix, do the above for along the first 
+## non-singleton dimension. If the optional argument @var{dim} is
+## given, operate along this dimension.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
 ## Description: Compute ranks
 
-## This code is rather ugly, but is there an easy way to get the ranks
-## adjusted for ties from sort?
+## 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)
+
+  if (nargin != 1 && nargin != 2)
+    usage ("ranks (x, dim)");
+  endif
 
-function y = ranks (x)
-
-  if (nargin != 1)
-    usage ("ranks (x)");
+  nd = ndims (x);
+  sz = size (x);
+  if (nargin != 2)
+    %% Find the first non-singleton dimension
+    dim  = 1;
+    while (dim < nd + 1 && sz (dim) == 1)
+      dim = dim + 1;
+    endwhile
+    if (dim > nd)
+      dim = 1;
+    endif
+  else
+    if (! (isscalar (dim) && dim == round (dim)) && dim > 0 && 
+	dim < (nd + 1))
+      error ("ranks: dim must be an integer and valid dimension");
+    endif
   endif
 
-  y = [];
-
-  [r, c] = size (x);
-  if ((r == 1) && (c > 0))
-    p = x' * ones (1, c);
-    y = sum (p < p') + (sum (p == p') + 1) / 2;
-  elseif (r > 1)
-    o = ones (1, r);
-    for i = 1 : c;
-      p = x (:, i) * o;
-      y = [y, (sum (p < p') + (sum (p == p') + 1) / 2)'];
-    endfor
+  if (sz (dim) == 1)
+    y = ones(sz);
+  else
+    ## The algorithm works only on dim=1, so permute if necesary
+    if (dim != 1)
+      perm = [1 : nd];
+      perm (1) = dim;
+      perm (dim) = 1;
+      x = permute (x, perm);
+    endif
+    sz  = size (x);
+    infvec = -Inf * ones ([1, sz(2 : end)]);
+    [xs, y] = sort (x);
+    eq_el = find (diff ([xs; infvec]) == 0);
+    if (isempty (eq_el))
+      [eq_el, y] = sort (y);  
+    else
+      runs = complement (eq_el+1, eq_el);
+      runs = reshape (y (runs), size (runs)) + 
+             floor (runs ./ sz (1)) * sz(1);
+      len = diff (find (diff ([Inf; eq_el; -Inf]) != 1)) + 1;
+      [eq_el, y] = sort (y);  
+      for i = 1 : length(runs)
+	p = y (runs (i)) + (len (i) - 1) / 2;
+	for j = 0 : len (i) - 1
+	  y (runs (i) + j) = p;
+	endfor
+      endfor
+    endif  
+    if (dim != 1)
+      y = permute (y, perm);
+    endif
   endif
 
 endfunction