bugGNU Octave - Bugs: bug #63962, perms - Performance - Usage of...

 
 

bug #63962: perms - Performance - Usage of native C++ algorithm helpful

Submitter:  None
Submitted:  Sun 26 Mar 2023 05:25:02 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  koerhen Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * GNU/Linux Fixed Release:  9.1.0
Planned Release:  9.1.0
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 30 Mar 2023 10:14:05 PM UTC, comment #22: 

No complaints from CI in 48 hours. Closing as fixed.

@Hendrik: Thank you for the patch, and I'll be happy to support you in more patches in future.

@jwe, I'll start a Discourse thread about moving functions (not just perms) from libinterp to liboctave.

Arun Giridhar <arungiridhar>
Group Member
Tue 28 Mar 2023 05:46:28 PM UTC, comment #21: 

@jwe: What you say makes sense. I had been looking for a way to avoid the long list of string comparisons to call the functions of the appropriate type, but only finding, say, isinteger() and not isint32(), I didn't press it too much. It didn't occur to me to augment the octave_value class and Array class themselves, which would provide the required polymorphism for free.

I'd already pushed the patch before your comment, so let's treat the move from libinterp/corefcn to liboctave/array as a new project, along with the other DEFUNs you have in mind.

Btw, what's the reason for wanting to minimize the heavy lifting within DEFUNs?

Arun Giridhar <arungiridhar>
Group Member
Tue 28 Mar 2023 04:36:39 PM UTC, comment #20: 

Instead of using a series of if/else if statements or switch on data types, I'd prefer to see a perms function added to the octave_value class, and then dispatch using a virtual function in the octave_base_value class hierarchy.  If I understand correctly, most of those would just forward to the template function (currently GetPerms) but it looks like that could be moved to liboctave as part of the Array class or as a separate template function defined outside of class itself.

This arrangement would have the advantage of allowing perms to work for octave_value objects without having to call a DEFUN function, and to also be available for Array<T> objects directly, not just for octave_value objects.

As much as possible, I'd like to see DEFUN functions to do as little as possible.  Usually that means argument processing and and then calling octave_value methods.  Or, in some cases when the function only works on a limited number of types, processing the values and calling liboctave or standard library functions to do the real work.

It's not absolutely necessary that we make these changes before committing this patch, but it is something that I'd like to make happen, not just for this function, but for all the DEFUN functions in Octave.

John W. Eaton <jwe>
Group administrator
Tue 28 Mar 2023 04:13:57 PM UTC, comment #19: 

Thanks. I made a couple minor changes (removed some extra parentheses and went back to "cater for" -- I had no idea that it had the same meaning as "cater to" because somehow I'd always used "cater for" as in provide food for an event. Learned something new today) and pushed to https://hg.savannah.gnu.org/hgweb/octave/rev/7782d1ead0a0

Marking as Ready For Test.

Arun Giridhar <arungiridhar>
Group Member
Tue 28 Mar 2023 02:59:24 PM UTC, comment #18: 

@Arun-

Thanks. Looks good and I am fine with the changes including the warning and then giving back an empty octave value.

For the commit log you can use my full name "Hendrik Koerner <koerhen@web.de>"

i)
I did tow last content changes: Changing the AND order of the check for "unique" for cells allowing for a faster execution time as a simple int comparison is faster than an octave class "is_equal" method call comparison - note that it will not matter much as it is already very fast in C++.
((Ar[i].is_equal (Ar[j]) && (myvidx[j] > myvidx[i])) -->

((myvidx[j] > myvidx[i]) && (Ar[i].is_equal (Ar[j]))

Also checking that the number of input arguments is not more than 2 (and not just less than one). I had this in the original code but lost it during the last iteration.


ii)
@ "Setting do_sort to false for the case of booleans?"
This is by design following the existing functionality and assert testings: for historical reasons logicals are not sorted but value but by order of input (in contrast to numericals) which makes some sense when you think about it.
E.g.
%!assert (perms (logical ([1 0])), logical ([0 1;, 1 0]))
%!assert (perms (logical ([0 1])), logical ([1 0; 0 1]))

So for historical compatibility I kept this but I now give some more comments
    // Do not sort logicals by value as the other numerical data types.
    // This reflects past design design decisions using the order positioning
    // for logicals which makes more sense than ordering logicals by value.


iii)
Note that there is a reopening and re-discussion of the perms ordering: https://savannah.gnu.org/bugs/?50426
I guess that eventually octave may support the Matlab convention to have reverse ordering by vector position for everything and not just for logicals (and never by value like now)
The change in perms.cc will be minimal:
a) Change the default value of bool do_sort = true to false (keep the functionality even though not used) in GetPerms
b) Remove providing the do_sort flag and the comment in the (clname == "logical")
b) Remove or change the assert (perms (1:5), perms ([2 5 4 1 3]'))

One could already change this now / wait until a decision is formed and solve 50426 together with 63962 referencing each other in the bugs.
Or solve it separately: so migrate now and make later a second migration for 50426 once agreed.
Your call and not much difference although in general a 1:1 relation between bug and migration from the latest versioned state to the solved state may be easier to trace.


PS: The switch from British to American English code comments is fine as well: In British English, to cater FOR something means to take it into account. In American English, you say you cater TO something. :-)


(file #54538)

Hendrik K <koerhen>
Tue 28 Mar 2023 01:21:31 PM UTC, comment #17: 

@Hendrik: Thanks for pulling in the BISTs from bug #63965.

Nearly there. I'm attaching a version for your review. I stripped some trailing spaces from the ends of lines, capitalized whole-line comments, etc. I also made the if-else ladder fit on one line where possible, shortening one variable name (unique_v --> uniq_v) to make things fit.

I see that you're returning 0 for the case that the input class was not recognized. I changed that to returning an empty octave_value with a warning that the class was not recognized. Is that OK with you?

Btw, why are you setting do_sort to false for the case of booleans? There is a comment that you're doing so consciously by design, but the reason isn't clear. Leaving it true makes it more consistent with "reverse lexicographic order", which the existing perms.m doesn't do properly for e.g. `perms ([false false true])`. Is it OK to change it to true? If not, a longer comment would be helpful.

Also, what is your preferred name and email for the commit log? Right now I have "Hendrik K <koerhen@web.de>".


(file #54537)

Arun Giridhar <arungiridhar>
Group Member
Tue 28 Mar 2023 06:17:57 AM UTC, comment #16: 

Attached is the patch created in my octave build directory. It works and passes all tests for perms when executing make check.
I hope everything is ok: it is my first patch using mercurial and for octave....
 

The patch already caters for identical unique and non-unique reverse lexicographic ordering (https://hg.savannah.gnu.org/hgweb/octave/rev/29e0d557a3be) including the documentation changes.

Also it supports "unique" for cell arrays in perms as discussed in this post and in https://savannah.gnu.org/bugs/?63965.
"unique" works for cell arrays containing all basic numerical types and strings albeit not for nested data containers where existing octave equality is not supported throwing an error like "error: binary operator '==' not implemented for 'cell' by 'cell' operations".
 
In theory the method "is_equal" for the generic data container octave_value could be extended to cover recursively nested data container and then the perms.cc code would work out of the box. But Octave is a math language, so supporting such a construct seems for me of a little value as it is not in scope of the normal use cases and a normal user.

Note that this "unique" support for cells is not a regression compared to current perms, so keeping it makes sense I think. But it can be deactivated by simply removing the parameter unique_v when calling GetPermsNoSort in the c++ code as well as removing the asserts covering this if decided as such.

In any case the implementation is fast enough requiring n * (n-1) / 2 mutual comparisons. I tested it with 170 identical values (it would be pure madness in reality giving so many values to perms) and it takes less than 5ms on my machine.  
As we are speaking about permutations, n=171 factorial(171) would blow the system even with doubles...


https://savannah.gnu.org/bugs/?63965 continues to track the cell array discussion.



  

@Arun-
thanks for your great support.

(file #54536)

Hendrik K <koerhen>
Mon 27 Mar 2023 10:33:22 PM UTC, comment #15: 

@Hendrik: I've attached a patch version of your function. The changes I made were these:

1. Copy the new file perms.cc to libinterp/corefcn
2. hg add libinterp/corefcn/permc.cc (add to repository)
3. Add perms.cc to libinterp/corefcn/module.mk (add to build system)
4. Delete scripts/specfun/perms.m
5. hg forget scripts/specfun/perms.m  (remove from repository) and .
6. Remove perms.m from scripts/specfun/module.mk  (remove from build system)
7. Edit perms.cc and change DLD to DEFUN. Use one of the DEFUN functions in data.cc as an example.
8. Copy include files at the top, using _isprimelarge_.cc as an example.
9. Build and fix compilation errors: remove oct.h, and fix the "template required" errors for all the "const Array foo = ...".
10. Add a print_usage at the start of the function.

Now it builds and passes "make check". You can compare the DLD format to the DEFUN format.

It still needs formatting work, which I will leave to you. The Octave code guidelines are here, for your reference: https://wiki.octave.org/C%2B%2B_style_guide

For benchmarking, this code:

tic; for n = 1:11, perms(1:n); end; toc

takes 6.5 seconds on the unpatched version and 1.2 seconds on the patched version, while this code:

tic; for n = 1:15, v = round((1:n) / sqrt(n)); perms(v, "unique"); end; toc

takes 22.4 seconds on unpatched and 2.05 seconds on patched. Well done!


(file #54535)

Arun Giridhar <arungiridhar>
Group Member
Mon 27 Mar 2023 05:20:41 PM UTC, comment #14: 

unique should work for cellstrings like {"foo", "b", "b"}, but definitely won't work for other cells.

Nicholas Jankowski <nrjank>
Group Member
Mon 27 Mar 2023 05:13:07 PM UTC, comment #13: 

Re handling cells etc, I think calling `isequal (A, B)` or its equivalent in C++ would solve some problems. If it fails that sort of test, then the cell array is too complicated for perms too. The isequal function is slow though, so any C++ equivalent would be preferred.

Right now the problem with perms.m being unable to handle cells with "unique" is that the unique() function does not accept cell input. I can do this for small n by comparing rows, but that won't scale for bigger inputs. Even sorting and comparing only consecutive rows won't scale. https://savannah.gnu.org/bugs/?63965  But it does work properly for say `perms ({"foo", "b", "b"}, "unique")`


Arun Giridhar <arungiridhar>
Group Member
Mon 27 Mar 2023 05:07:17 PM UTC, comment #12: 

Pushed this change: https://hg.savannah.gnu.org/hgweb/octave/rev/29e0d557a3be

Now both unique and non-unique return reverse lexicographic order consistently.

Arun Giridhar <arungiridhar>
Group Member
Mon 27 Mar 2023 03:48:52 PM UTC, comment #11: 

to close off my topic derailment -

a decade ago it was agreed that matlab didn't dictate ordering and octave likewise shouldn't try to chase undoc'd behavior. (bug #30114)

at some point in the pre-2015 timeframe, matlab added a declaration that ordering was reverse lexicographic. but there was matlab bug report such that 2013-15 matlab didn't follow it's own specification.  Apparently that was fixed after 2015, acknowledged by bug #50426 and octave implemented the equivalent of sortrows(perms(x)) to follow suit.

However, as discussion around that matlab bug in https://www.mathworks.com/matlabcentral/answers/248750-the-working-of-perms#comment_316425
it is explained that "Reverse lexicographical order" does not refer to the value of the inputs. it refers to the ordering of the inputs.  so (checking with v2022b) an input of perms([1 2 3]) as expected produces:


>> perms([1 2 3])
ans =
     3     2     1
     3     1     2
     2     3     1
     2     1     3
     1     3     2
     1     2     3


while perms ([2 1 3]) appears as:


>> perms([2 1 3])
ans =
     3     1     2
     3     2     1
     1     3     2
     1     2     3
     2     3     1
     2     1     3


if it was value-based reverse lexicographical sorting, both would be the same and match sortrows(perms([2 1 3])), as octave currently does.

However, Matlab is doing order-based sorting.  (two examples above have the same positional sorting if you just swap with the 2 and the 1).

in conclusion - matlab's perms is now sorting consistent with it's stated documentation, if you understand it's wording to not simply mean what it says.  Octave, doesn't currently do this.

Bringing it back home:  i think unique should match similar to the base sorting case. I'm not sure how to do that with low overhead, but assume the order passed to unique could be retained with the I J outputs. separately, i'll comment about the octave-matlab mismatch over on the other report.

Nicholas Jankowski <nrjank>
Group Member
Mon 27 Mar 2023 03:21:22 PM UTC, comment #10: 

(now seeing some older discussions about perms ordering, and the wisdom of trying to follow undocumented matlab...) Older versions of the Matlab doc (<2005) didn't make any claims about ordering. Sometime later while the output from matlab never changed, the statement i quoted below about "in reverse lexicographic order" appeared in the public docs.  looking at the results of perms with different input ordering, i'm not quite sure how they define 'reverse lexicographical ordering.:


>> perms(1:4)

ans =

     4     3     2     1
     4     3     1     2
     4     2     3     1
     4     2     1     3
     4     1     3     2
     4     1     2     3
     3     4     2     1
     3     4     1     2
     3     2     4     1
     3     2     1     4
     3     1     4     2
     3     1     2     4
     2     4     3     1
     2     4     1     3
     2     3     4     1
     2     3     1     4
     2     1     4     3
     2     1     3     4
     1     4     3     2
     1     4     2     3
     1     3     4     2
     1     3     2     4
     1     2     4     3
     1     2     3     4

>> perms([2 3 1 4])

ans =

     4     1     3     2
     4     1     2     3
     4     3     1     2
     4     3     2     1
     4     2     1     3
     4     2     3     1
     1     4     3     2
     1     4     2     3
     1     3     4     2
     1     3     2     4
     1     2     4     3
     1     2     3     4
     3     4     1     2
     3     4     2     1
     3     1     4     2
     3     1     2     4
     3     2     4     1
     3     2     1     4
     2     4     1     3
     2     4     3     1
     2     1     4     3
     2     1     3     4
     2     3     4     1
     2     3     1     4

>> perms([4 3 2 1])

ans =

     1     2     3     4
     1     2     4     3
     1     3     2     4
     1     3     4     2
     1     4     2     3
     1     4     3     2
     2     1     3     4
     2     1     4     3
     2     3     1     4
     2     3     4     1
     2     4     1     3
     2     4     3     1
     3     1     2     4
     3     1     4     2
     3     2     1     4
     3     2     4     1
     3     4     1     2
     3     4     2     1
     4     1     2     3
     4     1     3     2
     4     2     1     3
     4     2     3     1
     4     3     1     2
     4     3     2     1


other than 'the inputs vector is always last, and you can watch the permutation flips progress upward starting from the rightmost value', I think the only real 'rule' is "if your inputs starts sorted lowest to highest, your output will be in reverse lexicographical order"

Nicholas Jankowski <nrjank>
Group Member
Mon 27 Mar 2023 03:06:20 PM UTC, comment #9: 

are you suggesting changing the order of the regular or 'unique' results?  matlab docs for perms specifies that they return in " in reverse lexicographic order"  which seems to be matched by the default perms behavior. If the overhead isn't too painful I'd support having the (octave only) 'unique' option do the same.

Nicholas Jankowski <nrjank>
Group Member
Mon 27 Mar 2023 02:01:17 PM UTC, comment #8: 

Another thing: Currently we change the permutation sort order when using "unique" for numerical values. Having the same sort order with or without "unique" seems more intuitive.

Any special reason for this ? Can we change this ?



perms ([1, 2, 3])
          ans =

             3   2   1
             3   1   2
             2   3   1
             2   1   3
             1   3   2
             1   2   3

perms ([1, 2, 3],"unique")
          ans =

             1   2   3
             1   3   2
             2   1   3
             2   3   1
             3   1   2
             3   2   1

Note: I already implemented the same existing reversing logic for cell arrays.
perms ({0.1, "foo", "foo"})
          ans =
          {
            [1,1] = foo
            [2,1] = foo
            [3,1] = foo
            [4,1] = foo
            [5,1] = 0.1000
            [6,1] = 0.1000
            [1,2] = foo
            [2,2] = 0.1000
            [3,2] = foo
            [4,2] = 0.1000
            [5,2] = foo
            [6,2] = foo
            [1,3] = 0.1000
            [2,3] = foo
            [3,3] = 0.1000
            [4,3] = foo
            [5,3] = foo
            [6,3] = foo
          }

perms ({0.1, "foo", "foo"},"unique")
          ans =
          {
            [1,1] = 0.1000
            [2,1] = foo
            [3,1] = foo
            [1,2] = foo
            [2,2] = 0.1000
            [3,2] = foo
            [1,3] = foo
            [2,3] = foo
            [3,3] = 0.1000
          }

Hendrik K <koerhen>
Mon 27 Mar 2023 10:33:58 AM UTC, comment #7: 

I would suggest dropping the class comparison, as for numerical types one would expect the typical octave results:

int32(1) == int8(1) or double(1) == single(1)
1

Or:
class(int32(1)) = int32 is not equal to
class(int8(1)) == int8

and also when a cell contains a cell, the class is identical but the comparison will also fail with an error message
a = {0.1, "foo"};
b = {0.1, "foo1"};
c = {0.1, "foo", a};
d = {0.1, "foo", b};
c{3} == d{3}
error: binary operator '==' not implemented for 'cell' by 'cell' operations

Hendrik K <koerhen>
Mon 27 Mar 2023 10:02:56 AM UTC, comment #6: 

Fair point. The function unique() does not handle heterogenous cell arrays either.

One way to extend equality for this context would be something like "if (classname(foo(i)) == classname(foo(j)) && all (foo(i) == foo(j))" in pseudocode, but anything more complicated wouldn't be necessary.

Arun Giridhar <arungiridhar>
Group Member
Mon 27 Mar 2023 09:43:04 AM UTC, comment #5: 

Everything what can be done in Octave can be relatively easily be done in OCT C++, so yes it can be done and yes I am happy to support.

Note that the problem is tough though, as Cells and structs can contain anything including other cells and struct. Not even Octave can do a comparison on these things. Example:

a = {0.1, "foo"};
b = {0.1, "foo", a};
b{1} == b{3}
error: binary operator '==' not implemented for 'scalar' by 'cell' operations

But if one limits the input to perms (or throws error messages) of cells which contain only fundamental numerical data types and strings (using 'all' allows comparing strings) it can be done.
I suggest an n:n comparison to detect the same values and create the index vector, so no sort (which does not work). Some overhead will occur but the number of items for perms must be small anyway....

Thoughts ?

Hendrik K <koerhen>
Mon 27 Mar 2023 08:49:08 AM UTC, comment #4: 

In the course of reviewing your patch, I realized that the existing perms.m does not do unique properly for cell inputs. The workaround is to construct a numerical index vector and pass that to perms, then use it to index into the cell array. That should ideally be done by perms though, so that's a bug. I'll file a separate report.

Would your code be able to accommodate unique for cells by doing that sort of numerical index translation? If not, an alternative is to do that index conversion in perms.m, pass it to the new perms.cc, and then translate the result back in perms.m. Thoughts?

Arun Giridhar <arungiridhar>
Group Member
Mon 27 Mar 2023 03:12:49 AM UTC, comment #3: 

That is great news.
Some support would be very useful as I am not familiar with mercurial to create patches nor with the C++ software functions as well as the core development architecture of octave.
 
So e.g. it appears that native C++ octave functions use the template DEFFUN using 'doc' for documentation and other <includes> instead of just DEFUN_DLD <octave/oct.h>.
Also makefiles and documentation driver files would need to be adapted to replace perms.m with perms.cc, the directory location for perms.cc needs to be chosen etc etc.

 
Some guidance/support of what additional files (see above some points coming to my mind) need to be changed and how as well as what needs to change within the current perms.cc when switching from DEFFUN_DLD to DEFFUN would be highly appreciated.


I attach the current fully working perms.cc DEFFUN_DLD version.



(file #54531)

Hendrik K <koerhen>
Sun 26 Mar 2023 01:37:01 PM UTC, comment #2: 

Right, I see that std::next_permutation does not generate repetitions, so it won't spike memory usage: https://en.cppreference.com/w/cpp/algorithm/next_permutation

There's a note there that its behavior changes in C++20 (becomes a constexpr) but that should not be a concern.

Post your code when ready.

Arun Giridhar <arungiridhar>
Group Member
Sun 26 Mar 2023 12:58:19 PM UTC, comment #1: 

You can always post a Mercurial patch if you're familiar with them. I'm happy to review it. If not, post the C++ code you wrote and I will work with you to make it a patch.

Btw, your message isn't clear on whether your C++ replacement code also does "unique" or not. In particular, does it generate all repetitions and then discard duplicates (high memory consumption intermediate) or does it generate only the unique permutations?

Arun Giridhar <arungiridhar>
Group Member
Sun 26 Mar 2023 05:25:02 AM UTC, original submission:  

The function perms have been recently reworked (to cater for "unique") and enhanced for performance but it is not really fast and the algorithm seems to be pretty complex.

Or C++ offers a standard permutation algorithm std::next_permutation which is standard and is very performant.

As perms was for me a bottleneck, I implemented a C++ routine which is straightforward (the key algorithm is only 5 lines):

    octave_idx_type i = 0;
    do {
      for (octave_idx_type j = 0; j < m; j++)
         Res[ i + j*n ] = Myarray[myvidx[j]];
      i++;
    } while (std::next_permutation(myvidx, myvidx + m,std::greater<int>()));

Thanks to the C++ template approach this code snippet supports all octave datatypes including the non documented perms features for cell arrays and even structs.
The overall  C++ code is about 200 lines, very readable, easy to change/adapt and materially faster than the octave language implementation.

It is a 1:1 replacement of the current perms function keeping even the sorting order identical (even though for permutations not required). All (pretty extensive) assert test functions of perms are passed.

Performance wise it improved by factor 4 to 50 with the sweat spot of permutation between 4 to 7:
N       Performance Factor
1 4
2 5
3 7
4 50
5 60
6 45
7 14
8 6
9 6
10 4
11 4
 
I am happy to share the code and get it incorporated in the mainline for the benefits of the community if there is interest.


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #54538:  patch_63962_v4 added by koerhen (21KiB - application/octet-stream)
file #54537:  perms_v3.patch added by arungiridhar (21KiB - text/x-patch)
file #54536:  patch_63962 added by koerhen (21KiB - application/octet-stream)
file #54535:  perms.patch added by arungiridhar (20KiB - text/x-patch)
file #54531:  perms.cc added by koerhen (11KiB - text/x-c++src)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by koerhen (Updated the item)
  • -email is unavailable- added by arungiridhar (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-03-30 arungiridhar Fixed ReleaseNone 9.1.0
        Planned ReleaseNone 9.1.0
    2023-03-30 arungiridhar StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2023-03-28 arungiridhar StatusNone Ready For Test
    2023-03-28 koerhen Attached File- Added patch_63962_v4, #54538
    2023-03-28 arungiridhar Attached File- Added perms_v3.patch, #54537
    2023-03-28 koerhen Attached File- Added patch_63962, #54536
    2023-03-27 arungiridhar Attached File- Added perms.patch, #54535
    2023-03-27 koerhen Attached File- Added perms.cc, #54531

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code