Opened 4 months ago

Last modified 6 weeks ago

#30319 new defect

wrong intersection of polytopes

Reported by: vdelecroix Owned by:
Priority: critical Milestone: sage-9.3
Component: geometry Keywords:
Cc: jipilab Merged in:
Authors: Reviewers:
Report Upstream: Reported upstream. Developers acknowledge bug. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description

As reported in this sage-devel thread, the intersection of polytopes is sometimes very wrong

a = Polyhedron([[0, -1, 1], [1, -1, 1], [1, 1, -1]])
b = Polyhedron([[0.0, -0.5, 1.5], [1.0, -0.5, 1.5], [1.0, 1.5, -0.5]])
c = a.intersection(b)

Here c should have been empty but is equal to b.

Change History (14)

comment:1 Changed 4 months ago by mkoeppe

  • Priority changed from major to critical

comment:2 Changed 4 months ago by gh-kliem

Looks like a cdd bug.

sage: a = Polyhedron([[0, -1, 1], [1, -1, 1], [1, 1, -1]])
sage: b = Polyhedron([[0.0, -0.5, 1.5], [1.0, -0.5, 1.5], [1.0, 1.5, -0.5]])
sage: a = a.change_ring(RDF)
sage: new_eqns = a.equations() + b.equations()
sage: new_ieqs = a.inequalities() + b.inequalities()
sage: Polyhedron(eqns=new_eqns, ieqs=new_ieqs[1:2]+new_ieqs[4:5], verbose=True)
---- CDD input -----
H-representation
linearity 2 1 2
begin
 5 4 real
 0.0 0.0 1.0 1.0
 -1.0 0.0 1.0 1.0
 -1.0 2.0 -1.0 0.0
 -1.0 4.0 -2.0 0.0
 1 0 0 0
end

---- CDD output -----
Canonicalize the matrix.
Implicit linearity rows are:

Redundant rows are:1 3 5 

Nonredundant representation:
The new row positions are as follows (orig:new).
Each redundant row has the new number 0.
Each deleted duplicated row has a number nagative of the row that
represents its equivalence class.
 1:0 2:1 3:0 4:2 5:0
H-representation
linearity 1  1
begin
 2 4 real
 -1  0  1  1
 -1  4 -2  0
end

size = 5 x 4
Number Type = real

---- CDD input -----
Canonicalize the matrix.
Implicit linearity rows are:

Redundant rows are:1 3 5 

Nonredundant representation:
The new row positions are as follows (orig:new).
Each redundant row has the new number 0.
Each deleted duplicated row has a number nagative of the row that
represents its equivalence class.
 1:0 2:1 3:0 4:2 5:0
H-representation
linearity 1  1
begin
 2 4 real
 -1  0  1  1
 -1  4 -2  0
end

---- CDD output -----
The second representation:
V-representation
linearity 1  3
begin
 3 4 real
  0  1  0  0
  1  7.500000000E-01  1  0
  0 -5.000000000E-01 -1  1
end

Vertex incidence
ecd_file: Incidence of generators and inequalities
begin
  3    3
 1 -2 : 2 
 2 -2 : 3 
 3 -3 : 
end

Vertex adjacency
ead_file: Adjacency of generators
begin
  3    3
 1 1 : 2 
 2 1 : 1 
 3 0 : 
end

The first (input) representation
H-representation
linearity 1  1
begin
 2 4 real
 -1  0  1  1
 -1  4 -2  0
end

Facet incidence
icd_file: Incidence of inequalities and generators
begin
  3    3
 1 -3 : 
 2 -2 : 1 
 3 -2 : 2 
end

Facet adjacency
iad_file: Adjacency of inequalities
begin
  3    3
 1 -2 : 1 
 2 -2 : 2 
 3 -2 : 3 
end

size = 2 x 4
Number Type = real

---- CDD input -----
The second representation:
V-representation
linearity 1  3
begin
 3 4 real
  0  1  0  0
  1  7.500000000E-01  1  0
  0 -5.000000000E-01 -1  1
end
---- CDD output -----
The second representation:
H-representation
linearity 1  3
begin
 3 4 real
  1  0  0  0
 -1  4 -2  0
 -1  0  1  1
end

size = 3 x 4
Number Type = real

A 2-dimensional polyhedron in RDF^3 defined as the convex hull of 1 vertex, 1 ray, 1 line

For some reason cdd thinks that row 1 is redundant. If you drop an inequality, everything is fine.

Last edited 2 months ago by gh-kliem (previous) (diff)

comment:3 Changed 4 months ago by tmonteil

  • Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.

Could someone please report it upstream (i do not have a github account) ?

comment:4 Changed 4 months ago by gh-kliem

  • Report Upstream changed from Not yet reported upstream; Will do shortly. to Reported upstream. No feedback yet.

comment:5 Changed 2 months ago by dimpase

Note that with data in comment:2 one has

sage: Polyhedron(eqns=new_eqns)                                                                                                      
The empty polyhedron in RDF^3

so the emptiness should have been detected earlier.

comment:6 Changed 2 months ago by gh-kliem

Are you proposing that we test whether the equations are infeasible ourselves as a workaround?

Sounds like a workaround that should work for us.

Is cdd specifically excluding this case?

comment:7 Changed 2 months ago by gh-kliem

Something like

new_eqns_matrix = matrix(new_eqns)
feasible_directions = new_eqns_matrix.right_kernel_matrix()
if feasible_directions.column(0).is_zero():
    initialize_empty_polyhedron()

should work.

comment:8 Changed 2 months ago by dimpase

On the other hand, with RDF data and the internal representation used by cdd, it is not so hard to come up with (ugly) examples on which it can fail to produce a meaningful answer.

I'd say that computing c like in the ticket description should throw an error, as a is exact and b is not. Or at the very least print a warning like

sage: a = Polyhedron([[0, -1, 1], [1, -1, 1], [1, 1, -1]],base_ring=RDF) 
....: b = Polyhedron([[0.0, -0.5, 1.5], [1.0, -0.5, 1.5], [1.0, 1.5, -0.5]])                                                         
sage: a.intersection(b)                                                                                                              
/mnt/opt/Sage/sage-dev/local/lib/python3.8/site-packages/sage/geometry/polyhedron/backend_cdd.py:151: UserWarning: This polyhedron data is numerically complicated; cdd could not convert between the inexact V and H representation without loss of data. The resulting object might show inconsistencies.
  warn("This polyhedron data is numerically complicated; cdd could not convert between the inexact V and H representation without loss of data. The resulting object might show inconsistencies.")
A 2-dimensional polyhedron in RDF^3 defined as the convex hull of 3 vertices

comment:9 Changed 2 months ago by dimpase

it's also weird that no warning is printed in this case:

sage: a = Polyhedron([[0, -1, 1], [1, -1, 1], [1, 1, -1]])
....: b = Polyhedron([[0.0, -0.5, 1.5], [1.0, -0.5, 1.5], [1.0, 1.5, -0.5]])      
....: b.intersection(a)                                  
A 2-dimensional polyhedron in RDF^3 defined as the convex hull of 3 vertices

as if seeing the exact a to intersect b with convinces cdd that everything should be fine.

comment:10 follow-up: Changed 2 months ago by gh-kliem

The behavior of coercing base rings is the correct and expected behaviour, I would say.

The failure of cdd has nothing to do with inexactness:

sage: a = Polyhedron([[0, -1, 1], [1, -1, 1], [1, 1, -1]], backend='cdd')                                                                                                           
sage: b = Polyhedron([[0, -1/2, 3/2], [1, -1/2, 3/2], [1, 3/2, -1/2]], backend='cdd')                                                                                               
sage: a.intersection(b)                                                                                                                                                             
A 2-dimensional polyhedron in QQ^3 defined as the convex hull of 3 vertices

comment:11 Changed 2 months ago by gh-kliem

cdd for some reason determines that one of the equations is redundant. After this the calculation is correct. For this reason, no warning is printed. Because converting from Vrepresention to Hrepresentation works.

The output of cdd is numerical consistent, but incorrect.

comment:12 in reply to: ↑ 10 Changed 2 months ago by dimpase

Replying to gh-kliem:

The behavior of coercing base rings is the correct and expected behaviour, I would say.

The failure of cdd has nothing to do with inexactness:

OK, I see, you're right - I thought it's an inexactness issue, such as computing rank of an RDF matrix, sorry.

comment:13 Changed 2 months ago by gh-kliem

  • Report Upstream changed from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

comment:14 Changed 6 weeks ago by mkoeppe

  • Milestone changed from sage-9.2 to sage-9.3
Note: See TracTickets for help on using tickets.