[Octave-bug-tracker] [bug #36372] Improved ranks.m included:

[Dave Goel]

Follow-up Comment #15, bug #36372 (project octave):

Nir,

Ah, thanks for spotting that. I have now fixed the problem. The test expected
ranks(my).m to fail when nargin was greater than 2. It now does.

I have named it ranksmy currently since I maintain both ranks and ranksmy
separately, pending acceptance of this file.

It now passes all 13 tests with calls to its own "ranksmy."

Attaching the new file. This web interface wouldn't let me attach it in a
follow-up comment, so just pasting it here:

## 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} {} ranksmy (@var{x}, @var{dim})
##
## Return the ranksmy of @var{x} along the first non-singleton dimension
## adjust for ties. If the optional argument @var{dim} is given, operate
## along this dimension. @var{x} can optionally be a structure. In that
## case, its field named x is interpreted as the matrix, and its field
## named rtype is interpreted as the rank-type. The field rtype then
## donates the type of ranking: 0 or "fractional" for fractional, which
## is the default, 1 or "competition" for competiton ranking, 2 or
## "modified" for modified competition ranking, 3 or "ordinal" for
## ordinal ranking, 4 or "revordinal" for reverse ordinal,  and 5 or
## "dense" for dense ranking.
##
address@hidden deftypefn


## This function returns a result identical to GNU Octave's built-in
## ranksmy.m but should be much faster (For example, 1.4s vs. "no result
## in 33mins" when the input was a=round(10*rand(100123,100));
## Additionally, we also handle several ranking schemes. This function
## has been submitted to GNU Octave's mailing list and may become part
## of GNU Octave. [[LOCAL NOTES 20120503 : To compare this file to the
## latest octave version, see the file tmpranksmyoctave.m in this
## directory.]]


## Authors: KH <address@hidden>, Dave Goel <address@hidden>
## Description: Compute ranksmy


function y = ranksmy (x, dim)

  if (nargin < 1) || (nargin > 2);
    print_usage ();
  endif

  if isstruct(x);
    rtype=0;
    if isfield(x,"rtype");
      rtype=x.rtype;
      if isempty(rtype);
        rtype=0;
      endif
    endif
    x=x.x;
  else
    rtype=0;
  endif


  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 ("ranksmy: dim must be an integer and valid dimension");
    endif
  endif

  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);
    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;
        ## We could also have taken the avg. of _ordinal and _revordinal
        ## here, except that the current setup didn't separate out the
        ## logic of defining them. That would also hvae involved two
        ## calls to sort, one for ordinal, and another for revordinal.
      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
        rtype
        error("Illegal value of rtype specified.");
    endswitch
    y=NaN(size(lin));
    ## If input was a vector, we could have simply done this:
    ## y(ids)=lin;
    y(_ids(ids,sz))=lin;
    if (dim != 1)
      y = permute (y, perm);
    endif
  endif
endfunction

function idf=_ids(ids,sz);
  oo=ones(sz);
  allids=[{ids-1}];
  nn=numel(sz);
  for ii=2:nn;
    allids=[allids,{cumsum(oo,ii)-1}];
  endfor
  idf=allids{end};
  for jj=(nn-1):-1:1;
    idf=idf*sz(jj)+allids{jj};
  endfor
  idf+=1;
endfunction



function linnew=_dense (lin,sx,sz)
  infvec = -Inf * ones ([1, sz(2 : end)]);
  fnewp= logical(diff([infvec;sx]));
  linnew=cumsum(fnewp,1);
  linnew
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 



%!assert (ranksmy (1:2:10), [1:5]'')
%!assert (ranksmy (10:-2:1), [5:-1:1]'')
%!assert (ranksmy ([2, 1, 2, 4]), [2.5, 1, 2.5, 4])
%!assert (ranksmy (ones (1, 5)), 3*ones (1, 5)'')
%!assert (ranksmy (1e6*ones (1, 5)), 3*ones (1, 5)'')
%!assert (ranksmy (rand (1, 5), 1), ones (1, 5))

%% Test input validation
%!error ranksmy ()
%!error ranksmy (1, 2, 3)
%!error ranksmy ({1, 2})
%!error ranksmy (['A'; 'B'])
%!error ranksmy (1, 1.5)
%!error ranksmy (1, 0)
%!error ranksmy (1, 3)






    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?36372>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/