default library runs forever (or segfaults) in simple examples

[parrilo]

On some simple examples the default library seems to run forever, eventually running out of memory or segfaulting, while CDDLib works OK (on both Linux and Mac).

The example below is in exact arithmetic, but I've seen the same thing happening in floating point.

using Polyhedra                                                                                                         
using CDDLib                                                                                                            

# Octahedron                                                                                                                 
V = [1 0 0 ; 0 1 0 ; 0 0 1 ; -1 0 0 ; 0 -1 0 ; 0 0 -1];                                                                 

# Using CDDLib (works OK)                                                                                               
P = polyhedron(vrep(V),CDDLib.Library());                                                                               
hrep(P+P)                                                                                                               

# Using default library (runs forever, or segfaults)                                                                    
P = polyhedron(vrep(V));                                                                                                
hrep(P+P)

[blegat]

It seems to be a performance issue. With the debugging added in https://github.com/JuliaPolyhedra/Polyhedra.jl/pull/151, running

using Polyhedra

# Octahedron
V = [ 1  0  0
      0  1  0
      0  0  1
     -1  0  0
      0 -1  0
      0  0 -1]

voct = vrep(V)
doubledescription(voct + voct, verbose=3)

gives

Intersecting halfspace 1/19
Removing duplicates: 1 points, 8 rays and 0 lines.
Removing redundancy: 1 points, 1 rays and 3 lines.
After intersection:  1 points, 1 rays and 3 lines.
Intersecting halfspace 2/19
Removing duplicates: 1 points, 9 rays and 0 lines.
Removing redundancy: 1 points, 3 rays and 2 lines.
After intersection:  1 points, 2 rays and 2 lines.
Intersecting halfspace 3/19
Removing duplicates: 1 points, 10 rays and 0 lines.
Removing redundancy: 1 points, 5 rays and 1 lines.
After intersection:  1 points, 4 rays and 1 lines.
Intersecting halfspace 4/19
Removing duplicates: 1 points, 9 rays and 0 lines.
Removing redundancy: 1 points, 8 rays and 0 lines.
After intersection:  1 points, 8 rays and 0 lines.
Intersecting halfspace 5/19
Removing duplicates: 1 points, 16 rays and 0 lines.
Removing redundancy: 1 points, 14 rays and 0 lines.
After intersection:  1 points, 11 rays and 0 lines.
Intersecting halfspace 6/19
Removing duplicates: 1 points, 35 rays and 0 lines.
Removing redundancy: 1 points, 30 rays and 0 lines.
After intersection:  1 points, 15 rays and 0 lines.
Intersecting halfspace 7/19
Removing duplicates: 1 points, 36 rays and 0 lines.
Removing redundancy: 1 points, 28 rays and 0 lines.
After intersection:  1 points, 22 rays and 0 lines.
Intersecting halfspace 8/19
Removing duplicates: 1 points, 64 rays and 0 lines.
Removing redundancy: 1 points, 50 rays and 0 lines.
After intersection:  1 points, 39 rays and 0 lines.
Intersecting halfspace 9/19
Removing duplicates: 1 points, 264 rays and 0 lines.
Removing redundancy: 1 points, 183 rays and 0 lines.
After intersection:  1 points, 77 rays and 0 lines.
Intersecting halfspace 10/19
Removing duplicates: 1 points, 617 rays and 0 lines.
Removing redundancy: 1 points, 250 rays and 0 lines.
After intersection:  1 points, 111 rays and 0 lines.
Intersecting halfspace 11/19
Removing duplicates: 1 points, 581 rays and 0 lines.
Removing redundancy: 1 points, 537 rays and 0 lines.
After intersection:  1 points, 202 rays and 0 lines.
Intersecting halfspace 12/19

and then it takes several minutes before printing the next line. After some time, it gets to

Removing duplicates: 1 points, 4955 rays and 0 lines.
Removing redundancy: 1 points, 4474 rays and 0 lines.
After intersection:  1 points, 352 rays and 0 lines.
Intersecting halfspace 13/19
Removing duplicates: 1 points, 3048 rays and 0 lines.
Removing redundancy: 1 points, 1307 rays and 0 lines.
After intersection:  1 points, 377 rays and 0 lines.
Intersecting halfspace 14/19
Removing duplicates: 1 points, 1057 rays and 0 lines.
Removing redundancy: 1 points, 946 rays and 0 lines.
After intersection:  1 points, 535 rays and 0 lines.
Intersecting halfspace 15/19
Removing duplicates: 1 points, 4087 rays and 0 lines.
Removing redundancy: 1 points, 3134 rays and 0 lines.
After intersection:  1 points, 2532 rays and 0 lines.
Intersecting halfspace 16/19

The doubledescription function implemented in Polyhedra isn't as efficient as CDD yet. It doesn't even have the same complexity because currently, we do not track adjacency information (see https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/26). Therefore, when intersecting a V-rep with an halfspace, CDD splits the points in the ones on the left and right and computes the new points at the intersection by taking the intersection of each edge made of a point of the right and a point of the left which are adjacent. Because we have no adjacency information, we consider all the edge with one point on the right and one on the left even if they are not adjacent. That gives quadratic complexity in the number of points instead of linear complexity. Because of this we also generate many redundant points and we need to remove them by calling removeredundancy.