bugGNU Octave - Bugs: bug #65244, Simplify code complexity of...

 
 

bug #65244: Simplify code complexity of perms.cc

Submitter:  Hendrik K <koerhen>
Submitted:  Sat 03 Feb 2024 03:10:21 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Other
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  10.1.0 (current default)
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 14 Feb 2024 12:16:36 AM UTC, comment #21: 

It's been a couple of days without comment.  Marking bug as fixed and closing report.

Rik <rik5>
Group administrator
Sun 11 Feb 2024 06:55:34 AM UTC, comment #20: 

Thanks for this well-thought-out changeset.  I made a few changes and checked it in here https://hg.savannah.gnu.org/hgweb/octave/rev/349b4adf686a

Changes:

1) I was not comfortable with "long double" since this is not used frequently elsewhere in the Octave code base.  And, as you noted, there isn't a PC available that has flintmax (2^53) bytes of memory anyways so using "long double" is not necessary.  I reverted Factorial to using just "double".

2) The Octave coding convention is to avoid CamelCase so I renamed "IsEqual" to "isequal".

3) In general I try to make the minimum changes necessary and to mimic the existing coding conventions in the file being edited.  There was no functional need to change postfix increment (i++) to prefix increment (++i) so I reverted that change.

4) I expanded the commit message to mention each function changed and what was changed in that function.  Octave inherits this convention from the GNU coding standard.

Marking as Ready for Test.

Rik <rik5>
Group administrator
Sun 11 Feb 2024 02:18:56 AM UTC, comment #19: 

A new change set is attached.
The c++11 compatible approach depicted by Markus using template specialisation is now used instead of "if constexpr" requiring c++17.


The comments are amended to reflect the additional analysis and to include a note, that the code can be simplified once allowing c++17.

Thanks to everyone for their valuable contribution.

(file #55683)

Hendrik K <koerhen>
Sat 10 Feb 2024 05:23:48 AM UTC, comment #18: 

The discussion is interesting. To shed a little light on the practical implications, one needs to keep in mind that perms deals with factorials and factorials explode faster than anybody usually imagines and that time complexity is only an indication for CPU execution time. In a nutshell:

a) The code simplification using mutual comparisons is in nearly all cases faster(!) than using the (previous) sort approach.
b) Either removing constexpr or taking the C++11 compatible constexpr solution depicted by Markus is both fine.
c) For the "factorial" function in perms, one can use both doubles or long doubles.


Personally (suprise suprise - this reflects my coding past decisions...), I prefer using constexpr via Markus's solution and long doubles for the factorial function used by perms, but I am fine either way as in terms of practical implications, the different solutions for b) and (c) are not different...



More detailed analysis
======================

Whilest sorting has a lower time complexity (typical n*log n - worst case is also n*n by the way), ones needs to post-process the sort result in order to determine the number of unique values and how often they occur increasing both code complexity and time complexity. So one has in reality a complexity of n*logn + (n-1).
Also we do not need the sorted (interim) result, so we overfill the requirements and the sort overhead to reserve memory, fill with values, rearrange values etc. is "wasted".
 
 
Mutual comparison always requires n*(n-1)/2 complexity and does not have equivalent (unnecessary) overhead.

for i=1:20
 printf("n = %2d, Quicksort approach complexity %8.1f versus mutual comparison approach complexity %6d \n", i, i*log(i) + i - 1, i * (i-1)./2);
endfor
n =  1, Quicksort approach complexity      0.0 versus mutual comparison approach complexity      0
n =  2, Quicksort approach complexity      2.4 versus mutual comparison approach complexity      1
n =  3, Quicksort approach complexity      5.3 versus mutual comparison approach complexity      3
n =  4, Quicksort approach complexity      8.5 versus mutual comparison approach complexity      6
n =  5, Quicksort approach complexity     12.0 versus mutual comparison approach complexity     10
n =  6, Quicksort approach complexity     15.8 versus mutual comparison approach complexity     15

n =  7, Quicksort approach complexity     19.6 versus mutual comparison approach complexity     21
n =  8, Quicksort approach complexity     23.6 versus mutual comparison approach complexity     28
n =  9, Quicksort approach complexity     27.8 versus mutual comparison approach complexity     36
n = 10, Quicksort approach complexity     32.0 versus mutual comparison approach complexity     45
n = 11, Quicksort approach complexity     36.4 versus mutual comparison approach complexity     55
n = 12, Quicksort approach complexity     40.8 versus mutual comparison approach complexity     66
n = 13, Quicksort approach complexity     45.3 versus mutual comparison approach complexity     78
n = 14, Quicksort approach complexity     49.9 versus mutual comparison approach complexity     91
n = 15, Quicksort approach complexity     54.6 versus mutual comparison approach complexity    105
n = 16, Quicksort approach complexity     59.4 versus mutual comparison approach complexity    120
n = 17, Quicksort approach complexity     64.2 versus mutual comparison approach complexity    136
n = 18, Quicksort approach complexity     69.0 versus mutual comparison approach complexity    153
n = 19, Quicksort approach complexity     73.9 versus mutual comparison approach complexity    171
n = 20, Quicksort approach complexity     78.9 versus mutual comparison approach complexity    190


For up to n=6 the mutual comparison time complexity is actually better than using quicksort, so using mutual comparison will be faster.

I took n=7 and measured the processing time (using std::chrono::high_resolution_clock in the code of GetPerms mutual comparison):
n = 7; a = perms(floor(rand(1,n)*n),"unique");

The CPU proportion spent to determine the unique elements and their frequency was typically less than 1% for n=7. Factorials grow materially faster than n*log(n) so higher n will mean much lower proportion spent for this task.

Then I measured the GetPerms using the sort approach, and actually the (not used) sort overhead is pretty material leading to MORE time using sort compared to using the mutual comparisons for n = 7,8,9,10 

Note that the overhead to check parameters, extract them etc of perms itself before calling GetPerms was NOT included and which is proportionally material for smaller n.

Conclusions:
In practice mutual comparison is (nearly always) faster and the time spent in the code part to perform "unique" is only a tiny proportion of the total execution time, so it does not matter whether one uses constexpr or not.




Perms produces a matrix of the size (n!, n). The (minimal assuming int8) memory requirements are

s = {  "bytes", "kilobytes", "megabytes", "gigabytes", "terabytes", "pentabytes", "exabytes", "zettabytes" };
for i=1:20
  N = factorial(i) * i;
  s_i = 1;
  r_i = N;
  st = "";
  if N > flintmax(0.0)
    st = [ st " flintmax reached! "];
  endif
  if N >= intmax(uint64(0));
    st = [ st " 64 bit limit reached! "];
  endif
  while (r_i >= 1024)
    s_i++;
    r_i ./= 1024;
  endwhile
  printf("n = %2d, Required memory size: ~ %7.1f %10s %s\n", i, r_i, s{s_i},st);
endfor

n =  1, Required memory size: ~     1.0      bytes
n =  2, Required memory size: ~     4.0      bytes
n =  3, Required memory size: ~    18.0      bytes
n =  4, Required memory size: ~    96.0      bytes
n =  5, Required memory size: ~   600.0      bytes
n =  6, Required memory size: ~     4.2  kilobytes
n =  7, Required memory size: ~    34.5  kilobytes
n =  8, Required memory size: ~   315.0  kilobytes
n =  9, Required memory size: ~     3.1  megabytes
n = 10, Required memory size: ~    34.6  megabytes
n = 11, Required memory size: ~   418.7  megabytes
n = 12, Required memory size: ~     5.4  gigabytes
n = 13, Required memory size: ~    75.4  gigabytes
n = 14, Required memory size: ~     1.1  terabytes
n = 15, Required memory size: ~    17.8  terabytes
n = 16, Required memory size: ~   304.5  terabytes
n = 17, Required memory size: ~     5.4 pentabytes
n = 18, Required memory size: ~   102.4 pentabytes  flintmax reached!
n = 19, Required memory size: ~     2.0   exabytes  flintmax reached!
n = 20, Required memory size: ~    42.2   exabytes  flintmax reached!  64 bit limit reached!



Conclusions:
There is no practical difference whether for the factorial function one uses doubles or long doubles, except if one happens to have a computer with 100+ pentabytes of memory and a near eternity of CPU time to create and fill so much memory....
Either way is fine.

Hendrik K <koerhen>
Fri 09 Feb 2024 08:38:54 PM UTC, comment #17: 

I personally in my projects use different functions for different data types and am not sure if the provided patch is more readable and understandable than the original implementation.
However I appreciate any contribution even those that to my opinion aren't OK.
Moreover it is also good for me if one of the core maintainers says that something is reasonable.

Anonymous
Fri 09 Feb 2024 05:29:11 PM UTC, comment #16: 

The changeset adds the following comment:

+      // Mutual Comparison is used to detect duplicate values.
+      // Using sort would be possible for numerical values and be of
+      // n log n complexity instead of n * (n / 2). But sort
+      // is not supported for the octave-value container (structs/cells).
+      // As the perms element size n must be very small, any potential
+      // gains would be minimal as nearly all CPU is spent to create the
+      // actual permutations.


I didn't compare the performance. But the argument looked reasonable to me. Is that not the case?

Markus Mützel <mmuetzel>
Group administrator
Fri 09 Feb 2024 04:54:50 PM UTC, comment #15: 

It looks good to me. But it may be OK to have code duplication for different data structures.
I haven't tested the patch but I think it imposes overhead on Octave numeric types.
Currently GetPerms has complexity of n*log(n) and GetPermsNoSort has complexity of n^2.
If I correctly understand the proposed patch makes the complexity of both function equal to n^2.

Anonymous
Fri 09 Feb 2024 03:50:47 PM UTC, comment #14: 

Ah. Thank you for the example.
My previous tests were with an example that was too simple. (For more complex examples I didn't understand the assembly.)

Based on your example, would something like this work?

#include <vector>
#include <iostream>

struct octave_value
{
  octave_value ():val(0.0) {}
  explicit octave_value (double d) :val(d){}
  bool is_equal(const octave_value& b)
  {
    return val == b.val;
  }

  double val;
};

template <typename T>
bool is_equal_T (T a, T b)
{
  return a == b;
}

template <>
bool is_equal_T<octave_value> (octave_value a, octave_value b)
{
  return a.is_equal (b);
}

template <typename T>
bool perm ()
{
  std::vector <T> Ar;
  Ar.push_back(T(2.0) );
  Ar.push_back(T(3.0) );

  return is_equal_T<T>(Ar[0], Ar[1]);
}

int main()
{
  if (perm<double>())
    std::cout << "equal double" << std::endl;
  else
    std::cout << "unequal double" << std::endl;
  if (perm<octave_value>())
    std::cout << "equal octave_value" << std::endl;
  else
    std::cout << "unequal octave_value" << std::endl;
}


That might still allow to reduce part of the code duplication in the GetPerms and GetPermsNoSort templates.

If that works, we should still leave a note about the possible simplification using constexpr if once we allow C++17.

Markus Mützel <mmuetzel>
Group administrator
Fri 09 Feb 2024 01:45:42 PM UTC, comment #13: 

A similar code that doesn't compile on both gcc and clang:
 

#include <type_traits>
#include <vector>

struct octave_value
{
  octave_value ():val(0.0) {}
  explicit octave_value (double d) :val(d){}
  bool is_equal(const octave_value& b)
  {
    return val = b.val;
  }

  double val;
};

template <typename T>
bool perm ()
{
  std::vector <T> Ar;
  Ar.push_back(T(2.0) );
  Ar.push_back(T(3.0) );
  static constexpr bool is_same_type = std::is_same<T, octave_value>::value;

  if (is_same_type)
    return Ar[0].is_equal (Ar[1]);
  else
   return (Ar[0] == Ar[1]);
}

int main()
{
    perm<double>();
    perm<octave_value>();
}


Anonymous
Fri 09 Feb 2024 12:12:03 PM UTC, comment #12: 

I don't see the compiler error for the example in comment #9 (after fixing the typo at the end of the first line).

But looking at the generated assembly (without that intermediate constexpr value), at least GCC and Clang seem to produce the same code with "normal" if or with constexpr if even at -O0. (At least with a simple example where the assembly doesn't get too complicated.)
So, just omitting the constexpr is probably good enough to make the code compatible with C++11.

Markus Mützel <mmuetzel>
Group administrator
Fri 09 Feb 2024 11:30:37 AM UTC, comment #11: 
Anonymous
Fri 09 Feb 2024 11:25:07 AM UTC, comment #10: 

I think in C++11 the template tricks like enable_if and function overloading is used for such conditions.  "if (something constexpre)" is different than "if constexpr (somthing)".
I the example in comment #9 the compiler will throw error.

Anonymous
Fri 09 Feb 2024 11:02:44 AM UTC, comment #9: 


> Note that "the constexpr if statement" that is used in the patch is a feature of C++17.


Good point.

I wonder if something like the following would work with C++11 and be optimized to virtually the same code as the constexpr if:

static constexpr bool is_same_type = std::is_same<T, octave_value>::value:
if (is_same_type)
  IsEqual =  Ar[i].is_equal (Ar[j]);
else
  IsEqual =  (Ar[i] == Ar[j]);


(Not sure if `static` is necessary.)

Markus Mützel <mmuetzel>
Group administrator
Fri 09 Feb 2024 10:50:59 AM UTC, comment #8: 

Note that "the constexpr if statement" that is used in the patch is a feature of C++17.

Anonymous
Fri 09 Feb 2024 09:26:05 AM UTC, comment #7: 

The changes look good to me in general. I'm not so sure about the part concerning `long double` though. This would be the first time that that type would be used unconditionally in Octave.

They are used conditionally in `octave_int<T>`. But that is guarded by `OCTAVE_INT_USE_LONG_DOUBLE`.
ISTR that `long double` are a mess when it comes to cross-platform compatibility. Is it save to use them unconditionally here?

Markus Mützel <mmuetzel>
Group administrator
Fri 09 Feb 2024 02:49:11 AM UTC, comment #6: 

I have created a new patch incorporating the change of myidx to octave_idx_type.

In line with this change to use octave_idx_type for octave array indexing, I changed N_el (also related to indexing) from "int" to "octave_idx_type".
On top, I also changed the factorial from double to long double to prevent potential precision loss for the factorial result as it is used to size the resulting array.
And for code clarity, I modified the way on how the resulting array size is derived when using "unique" (build up a denominator and do one final division reflecting the analytical reasoning instead of successive individual divisions).

Note that in practice, the computer memory would be exhausted long beforehand before the any of these changes were to become relevant.
But the code soundness and clarity is improved, so worth the effort.

(file #55679)

Hendrik K <koerhen>
Thu 08 Feb 2024 05:16:15 PM UTC, comment #5: 

Thanks for the review.

I think that declaring myvidx as octave_idx_type would (slightly) improve the robustness of the code without much harm done in terms of used memory space:


In practical terms, the number of permutations explodes extremely quickly, so everything beyond myvidx of 12 will use all the available computer memory (e.g. myvidx = 15 would require hundreds of terabytes of memory)

Therefore I did not even think about anything beyond an "int" or even an "int8" for myvidx...

But through the parameter "unique", one could come up with (one) extreme case: a very long vector of identical values (size larger than an int) for which the perms solution with "unique" is the vector of identical values itself. This would not work with the current implementation.
Note that if there are two non-identical values in such a long vector beyond intmax values, the result would explode beforehand.


So changing myvidx to octave_idx_type would (slightly) increase the robustness but also especially the logical soundness/clarity of the code:

When one gets an Octave array, one should expect a maximum index fitting in an octave_idx_type variable and therefore use octave_idx_type for index variables for Octave input arrays. This is already done in the code.

But this should also be true for interim/dependent variables which are directly driven by the size of the input array (the case of myvidx we are looking at here) 

So good catch...

Hendrik K <koerhen>
Thu 08 Feb 2024 04:27:13 PM UTC, comment #4: 


> This problem seems like something we should examine and make a decision about for many Octave functions where I we use int to hold the value returned from numel.


This particular case is a bit different to the issue we touched on recently: The return value of numel is assigned to a variable of type `octave_idx_type`. Everything fine here.
But in a subsequent step, values up to numel are assigned to a buffer of type `int`. So, the potential integer overflow would happen a bit later (and a simple pattern match would probably not catch it).

Markus Mützel <mmuetzel>
Group administrator
Thu 08 Feb 2024 04:14:36 PM UTC, comment #3: 

"Alternatively, would it make sense to add a check that `ar_in.numel ()` doesn't exceed INT_MAX (or another reasonable limit)?"

This problem seems like something we should examine and make a decision about for many Octave functions where I we use int to hold the value returned from numel.

John W. Eaton <jwe>
Group administrator
Thu 08 Feb 2024 03:35:09 PM UTC, comment #2: 

This looks good to me. A great application of C++11 features to simplify the code. Also a good opportunity to remove dead code. 👍

Tangentially related: Would it make sense to change the type of `myvidx` to `octave_idx_type`? At least theoretically, `ar_in` could have more than INT_MAX elements. But I guess the problem would explode much earlier.

Alternatively, would it make sense to add a check that `ar_in.numel ()` doesn't exceed INT_MAX (or another reasonable limit)?

Markus Mützel <mmuetzel>
Group administrator
Sat 03 Feb 2024 03:12:55 AM UTC, comment #1: 

Changeset attached with change-set ID.

(file #55659)

Hendrik K <koerhen>
Sat 03 Feb 2024 03:10:21 AM UTC, original submission:  

Perms.cc had been originally created when the ordering of the perms output varied by class type and was (in retro perspective) not Matlab compatible.

With the clarification (see bug #50426 for the lengthy discussion about the reverse lexicographical ordering based on position and ignoring element value), the code complexity can be materially reduced (40% less LOC) by removing unused code and by using more CXX-11 features (std::iota and constexpr allowing to merge templates).

This simplification will reduce future maintenance efforts.




Hendrik K <koerhen>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #55683:  patch_65244_v2.cset added by koerhen (8KiB - application/octet-stream)
file #55679:  patch_65244_v1.cset added by koerhen (7KiB - application/octet-stream)
file #55659:  patch_65244.cset added by koerhen (6KiB - application/octet-stream)

 

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 jwe (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by koerhen (Submitted the item)
  •  

    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 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-02-14 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2024-02-11 rik5 StatusPatch Reviewed Ready For Test
    2024-02-11 koerhen Attached File- Added patch_65244_v2.cset, #55683
    2024-02-09 koerhen Attached File- Added patch_65244_v1.cset, #55679
    2024-02-08 mmuetzel StatusNone Patch Reviewed
        Operating SystemGNU/Linux Any
        Planned ReleaseNone 10.1.0 (current default)
    2024-02-03 koerhen Attached File- Added patch_65244.cset, #55659

    Back to the top

    Powered by Savane 3.13-3e34.
    Corresponding source code