bugGNU Octave - Bugs: bug #36372, Improved ranks.m included:

 
 

bug #36372: Improved ranks.m included:

Submitted by:  None
Submitted on:  Wed 02 May 2012 04:50:36 PM UTC  
 
Category:  Octave Function Severity:  2 - Minor
Priority:  4 Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Dave Goel Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  dev
Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Tue 04 Dec 2018 05:25:29 AM UTC, comment #37: 

For point #1, Octave seems to be moving to mainting history and meta-information within the Mercurial version control system rather than within the files themselves.  So, we don't need explicit Author fields anymore.

For point #3, the current code uses a case statement to decode either a string or numeric argument at the same time.

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: unknown RTYPE '%s'", rtype);
endswitch

Is there more input validation you were thinking of?

Rik <rik5>
Project Administrator
Tue 04 Dec 2018 03:35:46 AM UTC, comment #36: 

Hi Rik,

(1) The description para. at the bottom seems to have gone missing.  Perhaps this needs to be added to the file:

## Author: Dave Goel <deego3@gmail.com>
## Description: Compute ranks using choice of methods when resolving ties.
## Thanks: Lachlan Andrew <lachanbis@gmail.com>, KH <Kurt.Hornik@wu-wien.ac.at>

or whatever else you guys see fit. (the code and algo. is completely new, but I maintained KH's original format.)

(2) I still have the revordinal method in my personal copy  (and my original submission included it). If you guys ever want to add it to octave's ranks, I will happily supply the code again.

(3) Perhaps ranktype could validaestring from   roptions = {"fractional", "competition", "modified", "ordinal", "revordinal"}. That's what I have done in my personal file.

Dave

Dave Goel <deego>
Thu 15 Nov 2018 08:13:28 PM UTC, comment #35: 

Hi Rik,

Much appreciated. Pleased to see it included. And, thanks again to all of you guys for all the hard work behind the scenes.  It's a thankless task merging other people's changes.

Thanks to all of you for Octave.

Dave

Dave Goel <deego>
Thu 15 Nov 2018 07:46:43 PM UTC, comment #34: 

@Dave: I apologize for this taking so long.  Even for Octave, 6 years is a bit much.  This bug report just slipped through the cracks.  I updated the lastest version of the patch to apply cleanly to the current development branch.  I also changed some of the code to follow Octave coding conventions.  You have been added to the list of contributors to Octave which shows up at the front of the manual.  I also announced the changes to ranks in the NEWS file (updated with every release) because I think this is an important development.  I checked everything in here https://hg.savannah.gnu.org/hgweb/octave/rev/31b443b5a6c1.

I put two FIXME notes in the code.  Note #1 is

    ## The algorithm works only on dim = 1, so permute if necessary.
    ## FIXME: Most all functions now accept a dim argument.
    ##        Would it be faster not to permute and use the dim argument
    ##        to sort, find, cumsum, etc.?

The second note is

    ## Offsets to map indices into each column to indices into the linear array.
    ## FIXME: Would sub2ind be faster here?
    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;

The 25X improvement on my machine is really wonderful, but perhaps there is still more that can be done.

Marking bug as fixed and closing report.

Rik <rik5>
Project Administrator
Fri 08 Sep 2017 04:31:43 PM UTC, comment #33: 

Minor comment. The changeset describes this as:
"Improve speed of ranks.m when dim!=2, and give choice of tie-breaks (bug #36372)"
I agree with the "give choice of tie-breaks" part.

But, regarding speed:

In fact, I think it improves the speed - and exponentially so - irrespective of dim, when there are (lots of) ties.

Dave

Anonymous
Fri 08 Sep 2017 03:41:20 PM UTC, comment #32: 

Avinoam,
You're right, I've got an updated one that applies clenly.
Attached.

(file #41761)

Philip Nienhuis <philipnienhuis>
Project Member
Tue 05 Sep 2017 08:46:54 PM UTC, comment #31: 

I've imported the patch just 5 days ago w/o issues ?
But I'll have a look

Philip Nienhuis <philipnienhuis>
Project Member
Tue 05 Sep 2017 08:12:09 PM UTC, comment #30: 

@Philip: Can you make a new patch?
I can't apply the current patch.
Thanks.

Avinoam Kalma <avinoam>
Project Member
Fri 19 Aug 2016 04:08:08 PM UTC, comment #29: 

can the latest cset be pushed?
I have this patch in my cross-builds for a long time.

Philip Nienhuis <philipnienhuis>
Project Member
Tue 03 May 2016 12:38:49 PM UTC, comment #28: 

I tried the latest patch to the dev tree.
The example cited in comment #8 and comment #9 run in 1.9 s instead of ~40 s w/o patch. Not bad.

>> test ranks
PASSES 13 out of 13 tests

Input validation works AFAICS. "ranks (1, 2, 3)" is accepted (although DIMS = 2) and returns 1, but I don't know if this is fine.

And the patch applies cleanly.

OS -> Any

Philip Nienhuis <philipnienhuis>
Project Member
Mon 08 Feb 2016 07:02:30 PM UTC, comment #27: 

That sneaky anonymous was me :-) (You can see who edited what in the wiki in of of the tabs of a page. AFAIK I was logged in when I commented).

I'll have another try.

Philip Nienhuis <philipnienhuis>
Project Member
Mon 08 Feb 2016 04:30:20 AM UTC, comment #26: 

Thanks to the anonymous person who tested this and wrote on the http://wiki.octave.org/Status_of_bugs page that the patch fails on all five hunks.

Which patch did you test?  The latest patch file #35951:  bug_36372.cset  seems to work for me.

Lachlan Andrew <lachlan>
Sat 09 Jan 2016 06:09:29 AM UTC, comment #25: 

Hi Lachlan

(1) To answer your question: I had the scores (for example, batting averages) of each individual in matrix. Each column represents a group (say, a team). The population is initially sorted by a parameter H (say, height.). I wanted to rank them, within each team, using this criterion: "Rank by the performance. In case of ties, prefer the tallest player."

(2) Re skipping revordinal, that's just fine; and again, really appreciate that you took the effort to create a patch for this. 

Dave

Dave Goel <deego>
Sat 09 Jan 2016 05:13:13 AM UTC, comment #24: 

It's interesting that it came up in your application.  Out of interest, what was that application?

I'm not sure what "Octave" will decide.  I don't get a say in that, but want to encourage the decision makers to actually include your patch.  Making it small will probably help.  If it is OK with you, I'll put up a changeset that omits reverse ordinal, and hope that they notice :)

Lachlan Andrew <lachlan>
Fri 08 Jan 2016 02:42:28 PM UTC, comment #23: 

Hi Lachlan,

Thanks for these hints. I have modified my code here to comply with these as well.

Regarding why we need reverse ordinal, I understand your points about leanness, and am fine one way or the other if Octave decides to not include it. In fact, I suspect the user can do some careful matrix reverses to get reverse ordinal out of ordinal.

My own reasons for inclusion would be:

(1) Completeness. It is a natural counterpart to ordinal.
(2) I actually need this for a project of mine. In fact, of the four ranking schemes, this is the one most naturally suited for that case. Implementing this is how I got started on this whole ranks project in the first place.

Many thanks again,
Dave

Dave Goel <deego>
Fri 08 Jan 2016 06:10:17 AM UTC, comment #22: 

Happy New Year again, Dave.

Thanks for the update.  I understand reverse ordinal ranking better now, but when would it be useful?  We want to keep Octave "lean", and only introduce options that will actually be used.  The others all seem useful, and the speedup is worthwhile.

I've modified your latest to comply with Octave's coding standards, and provided a patch relative the the current development version.

In future, could you try to use Octave's coding standards?  They require:
- space before () in a function call
- no space before () when indexing a matrix
- space after ',', and around '=', '+' etc
- Lines no longer than 80 characters

(They're all trivial whitespace issues, but uniformity makes the code easier to read.)

(file #35943)

Lachlan Andrew <lachlan>
Wed 06 Jan 2016 06:19:52 AM UTC, comment #21: 

Hi,

I am attaching a new iteration of ranks.m.

It uses a better example that clarifies reverse-ordinal ranking better.

Thanks,
Dave
--

(file #35930)

Dave Goel <deego>
Tue 29 Dec 2015 04:19:02 AM UTC, comment #20: 

Thanks Dave, but there's no need to list me as an author.  I'm just going through all these bug fixes in the hope that those with commit authority will apply them.

Happy New Year,
Lachlan

Lachlan Andrew <lachlan>
Sun 27 Dec 2015 06:34:56 PM UTC, comment #19: 

Hi Lachlan,

Thanks for your contributions to ranks, especially the algorithm to avoid a for loop when constructing the indexes over multiple dimensions.

I'm still new to hg changesets, but in my local copy, have adjusted the authors field thusly:

## Authors: Dave Goel <deego3@gmail.com>, Lachlan Andrew, KH <Kurt.Hornik@wu-wien.ac.at>

Dave

Dave Goel <deego>
Sat 05 Dec 2015 02:51:30 AM UTC, comment #18: 

Hi Lachlan,

(replying through the web interface this time.)

Thanks for your kind words, for going through the algos and for
the work you put in. Yeah, it was fun devising them :)

I browsed through your patch, and agree with your changes. Thanks. I had initially thought of not introducing a third argument in order to keep compatibility with any future matlab changes, but come to think of it, that way is indeed cleaner.

Thanks again,
Dave

Dave Goel <deego>
Thu 03 Dec 2015 11:40:29 PM UTC, comment #17: 

Dave, I like your new algorithms.  They're very elegant.  Even ignoring the performance improvement, I think that the reduced complexity is worthwhile.  The algorithm is also better with memory locality; I ran it on a machine with 512M, and the original ranks.m thrashes and takes 100 seconds to run even when dim=1, whereas the new one runs in a under 4 seconds.

I have:
1) tidied the code a little (in particular, it now takes the tie-break rule "rtype" as a third parameter, rather than as an element of a structure);
2) put it in Octave's coding style (as I understand it); and
3) made it a Mercurial change-set relative to the current development version.

If anyone wants to check the style, or test the code, feel free.  However, I think the patch is ready to apply.

(file #35638)

Lachlan Andrew <lachlan>
Tue 10 Feb 2015 09:45:02 PM UTC, comment #16: 

looks good now, thanks

Nir Krakauer <nir_krakauer>
Tue 10 Feb 2015 09:20:08 PM UTC, comment #15: 

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.
##
##@end 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 <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>
## 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)

Dave Goel <deego>
Tue 10 Feb 2015 08:56:34 PM UTC, comment #14: 

Dave-- I get the error I reported with the version you have as long as I change the filename to ranks.m before I call it (or else change the calls in the file to ranksmy). Possibly the reason you don't get the error is that your filename (ranksmy) doesn't match the function name in the file (ranks), so Octave is still using the built-in version for the tests. --Nir

Nir Krakauer <nir_krakauer>
Tue 10 Feb 2015 08:45:22 PM UTC, comment #13: 

Nir, That's strange. It passes all 13 tests for me. I wonder if we have different versions of ranksmy.m.  I am including the one I have below.

- Dave Goel (through the web interface.)

## 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}, @var{dim})
##
## Return the ranks 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.
##
##@end deftypefn

## This function returns a result identical to GNU Octave's built-in
## ranks.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 tmpranksoctave.m in this
## directory.]]

## Authors: KH <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>
## Description: Compute ranks

function y = ranks (x, dim)

  if (nargin < 1);
    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 ("ranks: 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 (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)'')
%!assert (ranks (rand (1, 5), 1), ones (1, 5))

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

Dave Goel <deego>
Tue 10 Feb 2015 04:09:23 PM UTC, comment #12: 

With 3.6.2, I get the same error Mike reported:

> test ranksmy

ans =  1
  ***** error ranksmy (1, 2, 3)
!!!!! expected <.> but got no error

Nir Krakauer <nir_krakauer>
Tue 22 Oct 2013 03:49:10 AM UTC, comment #11: 

Hi Mike,

I am now trying to reply via the interface. Regarding your comment about it failing a test, I just ran "test ranksmy". Unfortunately, on my 3.6.2, it passes all tests, so I am at a loss regarding that part. :( 

(I realize I need to get with the times, and build the octave source some day, but last time I tried on my debian stable, back in 2012, I was unable to even check it out, because it depended on some package not in stable.)

Thanks.
Dave
--

Dave Goel <deego>
Tue 22 Oct 2013 02:54:53 AM UTC, comment #10: 

Also confirmed that the attached version does indeed improve performance on the case in question. Input validation is broken, the test suite fails when nargin > 2, but the rest of the unit tests pass for me. Function needs to be cleaned up to conform to Octave coding style.

Mike Miller <mtmiller>
Project Administrator
Tue 22 Oct 2013 02:48:20 AM UTC, comment #9: 

Pasting email reply and attachment from Dave below. Dave, please try to reply on the bug tracker web interface, it does not respond at all to emails sent to it.

> Hi Mike
>
> Been a while since I submitted this, so I am a bit hazy.
>
> But, it looks like in the current version, the slowness comes out when
> sorting on dim 2, as below:
>
> octave:4> version
> ans = 3.6.2
> octave:5> a=round(10*rand(100123,100));
> octave:6> tic; ranksmy(a,2); toc
> Elapsed time is 1.8450469971 seconds.
>
>     (attaching the latest ranksmy.m I have, which is an adaptation of
>     KH's original ranks.m, and .
>
>     (1) Jordigh: The comments about "ugliness" were not mine, but
>     KH's. I have removed them, nonetheless.
>
>     Note (2) Secondly, this ranksmy.m, while equivalent to the
>     built-in ranks.m by default, also provides several ranking
>     variants in case of collisions: competition, fractional, modified
>     competition, ordinal, reverse ordinal, and dense ranking. Again,
>     unlike the original case, each is efficient, and none uses for
>     loops.
>
>
> Now, let's try, for comparison, the built-in function:
>
> octave:7> tic; ranks(a,2); toc
> Elapsed time is 48.160676956 seconds.
>
> That's a factor of 25.
>
> Plus, the difference increases greatly as you increase the size, because
> one vectorizes and the other loops.


Ok, confirmed, I also get about 50 seconds when dim = 2.

(file #29438)

Mike Miller <mtmiller>
Project Administrator
Tue 22 Oct 2013 02:20:08 AM UTC, comment #8: 

This bug has been idle for some time.

As I understand it, the original bug report here says that the performance of the ranks function is poor when operating on some input like the following:

a = round (10 * rand (100123, 100));
x = ranks (a);

However, on my system with the current development version, the call to ranks returns in 2-3 seconds. So I'm not sure if I'm not understanding what the OP was reporting or if the problem is simply no longer present.

Also it would be real nice if the changes could be attached as a file, even better as a patch to the current function, and better still as a mercurial changeset. The bug tracker seems to flatten all indenting so the text pasted here is less useful.

Mike Miller <mtmiller>
Project Administrator
Thu 10 May 2012 04:39:22 AM UTC, comment #7: 

Summary: New ranks.m included:

(1) Support for dense added for real this time.
(2) New support for reverse-ordinal sorting.


Greetings:

(A) The previous ranks.m I sent never supported one of the new advertized methods: "dense ranking." I have now implemented it, and included as (A) below.

(B) Incidentally, I found a much simpler way to implement (the original, fractional) ranks with Kurt's original logic of using sort twice. The difference is that it avoids any loops this time  But, note that this simpler logic comes at this expense: (a) this is still about twice more expensive than (A) above, because it uses sort twice; (b) It can't be used to implement competition and modified ranking. So, if we want to stick with Kurt's original logic, and forgo adding any new methods, we could simply switch to this method.). I am including this as "sampleranksmy.m" below.

Either way (A) or (B) is still orders of magnitude faster than the current version.

-----

(A): ranks.m with support for dense and reverse ordinal methods.

## Copyright (C) 1995-2012 Kurt Hornik; Copyright (C) 2012 Dave Goel
##
## This file is 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}, @var{dim})
##
## Return the ranks 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.
##
##@end deftypefn

## This function returns a result identical to GNU Octave's built-in
## ranks.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 tmpranksoctave.m in this
## directory.]]

## Authors: KH <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>
## Description: Compute ranks

function y = ranks (x, dim)

  if (nargin < 1);
    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 ("ranks: 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;
      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

function linnew=_revordinal (lin,sx,sz)
  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

%!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)'')
%!assert (ranks (rand (1, 5), 1), ones (1, 5))

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


(B) sampleranksmy.m - Kurt's original logic, but avoids any for loops:
## Copyright (C) 1995-2012 Kurt Hornik Copyright (C) 2012 Dave Goel
##
## This file is 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} {} ranks2my (@var{x}, @var{dim})
## Return the ranks 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 for "reverse ordinal ranking"  and 5 or "dense"
## for dense ranking.
## @end deftypefn

## Authors: KH <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>
## Description: Compute ranks

## This function returns a result identical to GNU Octave's built-in
## ranks.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 tmpranksoctave.m in this
## directory.]]

function y = sampleranksmy(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 ("ranks: 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
    [xsorted, xidsa] = sort (x);
    [ignore,  xidsb] = sort (x,'descend');
    xidsb=flipdim(xidsb,1);
    [ig, y1]=sort(xidsa);
    [ig, y2]=sort(xidsb);
    ## Now, y1 is ordinal sorting. y2 is reverse ordinal sorting. We
    ## simply need the mean of the two.
    y=(y1+y2)/2;
    if (dim != 1)
      y = permute (y, perm);
    endif
  endif
endfunction

Anonymous
Wed 09 May 2012 03:27:42 AM UTC, comment #6: 

True. My constraints are a bit weird. :(

> It is unusual to think that a chroot is easier or safer than backports.


My bad, You are right. I had meant to say fakechroot.

Anonymous
Wed 09 May 2012 03:20:06 AM UTC, comment #5: 

As I recall, you want to do all that because you don't want to enable the Debian backports repository and use hg 1.8 or newer?

Backports are really very safe and good. The backports repository is only going to install a newer version of Mercurial. Nothing else. It won't touch your system libraries, and its pin priority is lower than any other repository, so you will never install backports by default, only upon request. It is unusual to think that a chroot is easier or safer than backports.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Wed 09 May 2012 03:04:51 AM UTC, comment #4: 

Hi Jordi, No progress yet, unfortunately.

I would love to find time and get my hands dirty with chroot, debootstrap, etc. But, unfortunately, I have been called upon to meet some deadlines and that has kept me very busy. :(

Dave

Anonymous
Wed 09 May 2012 02:39:51 AM UTC, comment #3: 

Ping. Any progress on this? If not, I can attempt to distill a proper patch from your submissions so far.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Wed 02 May 2012 06:08:14 PM UTC, comment #2: 

Oops, I meant to refer you to this part of the Octave manual:

http://www.gnu.org/software/octave/doc/interpreter/Basics-of-Generating-a-Changeset.html

Also, some instructions on how to write the commit message:

http://www.octave.org/wiki/index.php?title=Commit_message_guidelines

In your case, I would greatly appreciate if you could explain a little in your commit message about what changed. Also, you may remove the part of the patch that says "this may not become a part of GNU Octave" and associated text. Part of that text should become a unit test (look at the bottom of the file for the syntax for writing tests).

Your comments about the code being ugly are fine, although I have a distaste for multiple exclamation marsk!!!

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Wed 02 May 2012 06:00:01 PM UTC, comment #1: 

Hi, thanks for the patch.

I know we've already made you jump through a bunch of hoops to give us this patch. I'm going to try to ask you to jump through a few more. If you don't want to, that's ok. However, bear in mind that if you don't do this further work, then someone else (e.g. me) will have to do it, and we probably don't care as much about your problem as you do, so it will take longer for us to do it.

That being said, if you can produce an hg changeset according to these instructions, it would be hugely beneficial:

http://www.gnu.org/software/octave/doc/interpreter/How-to-Contribute.html
http://jordi.inversethought.com/blog/how-to-write-a-patch-for-octave/

If those instructions are not clear, please let me know.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Wed 02 May 2012 04:50:36 PM UTC, original submission:  

Hi.

I believe that ranks.m (Though I originally tested from version 3.2.4,
the dev. version's algo is basically unchanged) is very inefficient
except when the elements of the input matrix are distinct.

For example,  for this input, a=round(10*rand(100123,100));

my alternative finished computation in 1.4s while the current function
hasn't yet returned an answer in the last 15-20 minutes.

For smaller inputs, I checked carefully, and the outputs do agree.

---
Additionally, the alternative also implements multiple ranking schemes:
fractional (default), standard competition ranking, modified competition
ranking, and finally, ordinal ranking; examples of all of which are seen
here: http://en.wikipedia.org/wiki/Ranking


Though I used a completely new algorithm, I managed to pen the
alternative as a derivative of the current function. Thus, attaching
(a) diff; (b) as the complete file.

The "diff -u" is included below. Below that follows the entire file. Thanks very much.


1a2

> ## Copyright (C) 2012 Dave Goel

6,7c7,8
< ## under the terms of the GNU General Public License as published by
< ## the Free Software Foundation; either version 3 of the License, or (at
---

> ## under the terms of the GNU General Public License as published by the
> ## Free Software Foundation; either version 3 of the License, or (at

10,13c11,14
< ## 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.
---

> ## 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.

18c19
<
---

> ##

20,27c21,40
< ## @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.
< ## @seealso{spearman, kendall}
< ## @end deftypefn
<
< ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
---

> ##
> ## @deftypefn {Function File} {} ranks (@var{x}, @var{dim}) Return the
> ## ranks of @var{x} along the first non-singleton dimension adjust for
> ## ties. If the optional argument @var{dim} is given, operate along this
> ## dimension. Third argument @var{rtype} 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 and 4 or
> ## "dense" for dense ranking. When using the third argument, you can
> ## supply us '' for @var{dim} to have it auto-computed. @end deftypefn
> ##
> ## This function returns a result identical to GNU Octave's built-in
> ## ranks.m but is 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 octave's mailing list and may (or may not)
> ## become part of GNU Octave.
> ##
> ##
> ## Authors: KH <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>

30,32d42
< ## 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..
34c44
< function y = ranks (x, dim)
---

> function y = ranks (x, dim, rtype)

36c46
<   if (nargin != 1 && nargin != 2)
---

>   if (nargin < 1);

40,41c50,51
<   if (! (isnumeric (x) || islogical (x)))
<     error ("ranks: X must be a numeric vector or matrix");
---

>   if nargin<3||isempty(rtype);
>     rtype=0;

43c53
<
---

>  

46c56
<   if (nargin != 2)
---

>   if (nargin<2)||isempty(dim);

48c58,64
<     (dim = find (sz > 1, 1)) || (dim = 1);
---

>     dim  = 1;
>     while (dim < nd + 1 && sz(dim) == 1)
>       dim = dim + 1;
>     endwhile
>     if (dim > nd)
>       dim = 1;
>     endif

50,52c66,69
<     if (!(isscalar (dim) && dim == fix (dim))
<         || !(1 <= dim && dim <= nd))
<       error ("ranks: DIM must be an integer and a valid dimension");
---

>     if (! (isscalar (dim) && dim == round (dim))
> && dim > 0
> && dim < (nd + 1))
>       error ("ranks: dim must be an integer and valid dimension");

66,80c83,105
<     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
---

>     sz=size(x);
>     [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"};
>         ## lin=lin;
>       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;
>

84a110

> endfunction

85a112,123

> 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;

88a127,148

> 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);
> 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
>
>

104d163
<

====================================================

## Copyright (C) 1995-2012 Kurt Hornik
## Copyright (C) 2012 Dave Goel
##
## This file is 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}, @var{dim}) Return the
## ranks of @var{x} along the first non-singleton dimension adjust for
## ties. If the optional argument @var{dim} is given, operate along this
## dimension. Third argument @var{rtype} 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 and 4 or
## "dense" for dense ranking. When using the third argument, you can
## supply us '' for @var{dim} to have it auto-computed. @end deftypefn
##
## This function returns a result identical to GNU Octave's built-in
## ranks.m but is 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 octave's mailing list and may (or may not)
## become part of GNU Octave.
##
##
## Authors: KH <Kurt.Hornik@wu-wien.ac.at>, Dave Goel <deego3@gmail.com>
## Description: Compute ranks

function y = ranks (x, dim, rtype)

  if (nargin < 1);
    print_usage ();
  endif

  if nargin<3||isempty(rtype);
    rtype=0;
  endif

 

  nd = ndims (x);
  sz = size (x);
  if (nargin<2)||isempty(dim);
    ## 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

  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);
    [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"};
        ## lin=lin;
      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=_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);
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 (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))
%!assert (ranks (rand (1, 5), 1), ones (1, 5))

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

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #41761:  bug_36372_ranks_v2.cset added by philipnienhuis (6KiB - application/octet-stream)
file #35951:  bug_36372.cset added by lachlan (6KiB - application/octet-stream)
file #35943:  bug_36372_patch.txt added by lachlan (6KiB - text/plain)
file #35930:  ranks.m added by deego (5KiB - text/x-objcsrc)
file #29438:  ranksmy.m added by mtmiller (6KiB - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by deego (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-11-15 mtmiller Carbon-CopyRemoved 80942 => -
    2018-11-15 rik5 StatusPatch Submitted => Fixed
        Open/ClosedOpen => Closed
    2017-09-08 philipnienhuis Attached File- => Added bug_36372_ranks_v2.cset, #41761
    2016-05-03 philipnienhuis Operating SystemGNU/Linux => Any
    2016-02-08 lachlan Severity3 - Normal => 2 - Minor
        Priority5 - Normal => 4
    2016-01-09 lachlan Attached File- => Added bug_36372.cset, #35951
    2016-01-08 lachlan Attached File- => Added bug_36372_patch.txt, #35943
    2016-01-06 deego Attached File- => Added ranks.m, #35930
    2015-12-03 lachlan Attached File- => Added bug_36372-ranks-patch.txt, #35638
    2013-10-22 mtmiller Attached File- => Added ranksmy.m, #29438
        StatusNeed Info => Patch Submitted
    2013-10-22 mtmiller CategoryNone => Octave Function
        StatusPatch Submitted => Need Info
    2012-05-09 jordigh StatusNone => Patch Submitted

    Back to the top


    Powered by Savane 3.5