"get_number_of_optimal_solutions()" sometimes returns too high value

[Hermann-SW]

Issue Details

I am using "get_number_of_optimal_solutions()" on a 96 points convex set. The set consists of all integer x,y,z with x^2+y^2+z^2==65 (blue) … … linear transformed with 3x3 matrix G. See the green convex hull: https://github.com/Hermann-SW/cgal4gp?tab=readme-ov-file#ternary-quadratic-form-width

Calling Width_3 with that pointset and then "get_number_of_optimal_solutions()" returns 8 instead of 2.

Source Code

You can just sync my cgal4gp repo where I make 7 pointset width related CGAL functions available for PARI/GP: https://github.com/Hermann-SW/cgal4gp

This is minimal example on PARI/GP side, can be done with same result in CGAL C++ as well:

~/cgal4gp$ gp -q tqf_width.gp
? S=SquaresRepresentations3(65);
? #S
96
? G=[-67, -189, -200; 6, 17, 18; 1, 3, 3];
? CGAL_Width([G*y~ | y<-S])
? CGAL_get_all_build_directions(dir)
? dir
[[81920, 917504, -40960], [204800, 2293760, -102400], [122880, 1376256, -61440], [81920, 917504, -40960], [-122880, -1376256, 61440], [-40960, -458752, 20480], [-163840, -1835008, 81920], [-163840, -1835008, 81920]]
? CGAL_get_number_of_optimal_solutions()
8
? 

Simplifying the 8 vectors shows that there are exactly two different directions and not 8 as reported by "CGAL_get_number_of_optimal_solutions()" :

? [d/gcd(d) | d<-dir]
[[10, 112, -5], [10, 112, -5], [10, 112, -5], [10, 112, -5], [-10, -112, 5], [-10, -112, 5], [-10, -112, 5], [-10, -112, 5]]
? 

Environment

• Ubuntu 22.04 64bit on Intel and Amd, 64bit Raspberry Pi OS on ARM
• Compiler: g++
• Release or debug mode: release
• Specific flags used (if any): see Makefile
• CGAL version: 5.4
• Boost version: 1.74.0
• Other libraries versions if used (Eigen, TBB, etc.): --

[Hermann-SW]

Wait, [10, 112, -5] and [-10, -112, 5] are "the same", only opposite. Wrt the two parallel width planes both are "the same". If that is true, CGAL_get_number_of_optimal_solutions() should return 1 in above example.

[Hermann-SW]

The PARI/GP way to identify unique vectors:

? dir=[[81920, 917504, -40960], [204800, 2293760, -102400], [122880, 1376256, -61440], [81920, 917504, -40960], [-122880, -1376256, 61440], [-40960, -458752, 20480], [-163840, -1835008, 81920], [-163840, -1835008, 81920]];
? Set([d/gcd(d)*if(d[1]<0,-1,1) | d<-dir])
[[10, 112, -5]]
? 

Normalize by d/gcd(d) and discard opposite direction vectors.

[Hermann-SW]

Here is a much simpler example where "CGAL_get_number_of_optimal_solutions()" returns twice (8) as many as correct (4):

pi@raspberrypi5:~/cgal4gp $ gp -q cgal4gp.gp 
? CGAL_Width([[1, 0, 0], [2, -1, 0], [2, 0, -1], [3, 0, 0], [2, 1, 0], [2, 0, 1]])
? CGAL_get_all_build_directions(dir)
? dir
[[4096, 4096, 4096], [4096, -4096, 4096], [-4096, -4096, 4096], [-4096, 4096, 4096], [-4096, 4096, -4096], [4096, 4096, -4096], [4096, -4096, -4096], [-4096, -4096, -4096]]
? CGAL_get_number_of_optimal_solutions()
8
? Set([d/gcd(d)*if(d[1]<0,-1,1) | d<-dir])
[[1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]
? 

[afabri]

Can you please provide a minimal example that only uses CGAL on gist.github.com

[Hermann-SW]

@afabri Sure, with demo execution in comment, matching to discussion of two cases above in this issue: https://gist.github.com/Hermann-SW/3fbb577ab0ef2abb0e89fb06965c3256

[afabri]

Hello Hermann,
When you add a -DSIMPLIFY as written in the section Large Numbers of the reference manual page, you will see for the second data set in the self contained example, that there are also only two vectors . So let me confirm that your found a bug. At least the width computation itself seems not to be compromised.

A quick fix in the code would probably to generate the planes, and run std::unique() with an equal that also returns true for opposite planes.

[Hermann-SW]

Hello Andreas,

thank you for "add a -DSIMPLIFY" (I am relatively new to CGAL). The results are better, but still wrong.

While the 8 entries listed in first case are all different, "1 1 1 1" and "-1 -1 -1 1" are the same width (width does not have a direction). So there are 4 distinct pairs, and get_number_of_optimal_solutions() should return 4. (see screenshot from July 16)

In second case there are two groups of 4 identical entries. And since width has no direction, all are "the same" (multiply with "-1"). So get_number_of_optimal_solutions() should return 1.

Do you agree?

pi@raspberrypi5:~/cgal4gp $ g++ -DSIMPLIFY 8358.cpp -lgmp 2>err
pi@raspberrypi5:~/cgal4gp $ ./a.out 
opt: 8 (should be 4)
1 1 1 1
1 -1 1 1
-1 -1 1 1
-1 1 1 1
-1 1 -1 1
1 1 -1 1
1 -1 -1 1
-1 -1 -1 1

opt: 8 (should be 1)
10 112 -5 1
10 112 -5 1
10 112 -5 1
10 112 -5 1
-10 -112 5 1
-10 -112 5 1
-10 -112 5 1
-10 -112 5 1

pi@raspberrypi5:~/cgal4gp $ 

[afabri]

I agree and that's why I wrote "a quick fix would be".

[afabri]

The SIMPLIFY only simplifies the quotients.

[Hermann-SW]

With latest revision of the gist https://gist.github.com/Hermann-SW/3fbb577ab0ef2abb0e89fb06965c3256 the workaround you described is implemented:

• compile with -DSIMPLIFY
• "normalize" away the direction
• std::sort()
• std::unique()

hermann@7950x:~$ g++ -DSIMPLIFY 8358.cpp -lgmp
hermann@7950x:~$ ./a.out 
opt: 8 (should be 4)
1 -1 -1 1
1 -1 1 1
1 1 -1 1
1 1 1 1

opt: 8 (should be 1)
10 112 -5 1

hermann@7950x:~$ 

This only works with -DSIMPLIFY, but it would be good if CGAL source could be fixed along the workaround so that .get_number_of_optimal_solutions() returns the correct number as well if compiled with that flag.

[Hermann-SW]

I moved PARI/GP workaround to C++, and that can be removed when this issue will be fixed: https://github.com/Hermann-SW/cgal4gp/commit/aa4083605426af308cd25b79168b97a3a55f3bf4

pi@raspberrypi5:~/cgal4gp $ gp -q < simple.gp 
- Width([[1, 0, 0], [2, -1, 0], [2, 0, -1], [3, 0, 0], [2, 1, 0], [2, 0, 1]])
- get_squared_width(num,denom)
4/3
4/3
- get_width_planes(e1,e2)
[1, 1, 1, -3] [1, 1, 1, -1]
- get_width_coefficients(A,B,C,D,K)
1 1 1 -3 -1
width-plane e1(e2) is given by the equation Ax+By+Cz+D(K)=0
- get_build_direction()
[1, 1, 1]
- get_all_build_directions(dir)
[[1, -1, -1], [1, -1, 1], [1, 1, -1], [1, 1, 1]]
- get_number_of_optimal_solutions()
4
pi@raspberrypi5:~/cgal4gp $ 

[afabri]

I looked into it and came to the conclusion that we only make it explicit in the documentation, that there can be opposite build direction vectors being returned, as well as vectors which are different but having the same direction. If we wanted to get rid of such multiples, we would have to do the simplification of numbers all the time. As it is documented that the number type must provide the modulo operation for the simplification, the removal of multiples would be conditionally available as well.

Is this acceptable for you? What you observe happens only for very regular polyhedra, right?

[Hermann-SW]

... If we wanted to get rid of such multiples, we would have to do the simplification of numbers all the time. As it is documented that the number type must provide the modulo operation for the simplification, the removal of multiples would be conditionally available as well.

Is this acceptable for you?

Perhaps the fix can be framed with

#ifdef SIMPLIFY
...
#endif 

because in that case the simplifications will happen anyway and are not caused by the fix. The only "overhead" of the fix in this case is normalization, std::sort() and std::unique() on a small list of candidates.

What you observe happens only for very regular polyhedra, right?

Yes.

From initial posting, in case of -DSIMPILFY with sort() and unique() the first example, for the green polyhedron (linearly transformed from the blue)

? [vecmin([(G*y~)[1] | y<-S]),vecmax([(G*y~)[1]|y<-S])]
[-2279, 2279]
? [vecmin([(G*y~)[2] | y<-S]),vecmax([(G*y~)[2]|y<-S])]
[-205, 205]
? [vecmin([(G*y~)[3] | y<-S]),vecmax([(G*y~)[3]|y<-S])]
[-35, 35]
? 

two opposite vectors will remain, with additional normalization exactly one.