JuliaPolyhedra / Polyhedra.jl

Polyhedral Computation Interface
Other
172 stars 27 forks source link

Assertion error in linear combination of rays #317

Closed Zinoex closed 1 year ago

Zinoex commented 1 year ago

I have been hitting an AssertionError when finding the H-rep from a V-rep. Specifically, when taking a linear combination of two rays on Line 361 of doubledescription.jl. I have narrowed it down to the following MVE. The concrete values for value1 and value2 are value1 = 4.0034642267983145e-13, value2 = 0.9999999999843894.

p = vrep([
           [-0.5387494623496795, -0.25960487027515805, -0.18239728626874152, 0.5038770788686975],
           [-0.5387493092049437, -0.25907155192885517, -0.1823972714041875, 0.5084592514097651],
           [-0.534924422025761, -0.029637292043959318, -0.1721646090548099, 0.5032870086786094],
           [-0.5349242688810253, -0.029103973697656493, -0.17216459419025582, 0.5078691812196772],
           [-0.5356218325215752, -0.03227549385025563, -0.1825609308207058, 0.5034125548692323],
           [-0.5356216793768396, -0.03174217550395278, -0.18256091595615176, 0.5079947274103],
       ])

doubledescription(p, nothing)

I've been trying to understand what is going on in the method to see if I could fix it myself including comparing the method to cddllib, and to my understanding, cddlib does abs(value2) * r1 + abs(value1) * r2, which different from what the comment in this method says. Is one method favorable to another and how can I fix the issue above?

blegat commented 1 year ago

Thanks for taking the time to create a MWE. There are two rays: r1, r2. r1 is part of the final representation and r2 is cut-off by halfspace 3. It is trying to add a ray that belongs to halfspace 3 by combining r1 and r2. The issue is that r1 belongs to the halfspace 3 (which you can see since value1 is 4e-13. So this https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/1cc19df35ea62efc3fefa162f709f04d09ea4c94/src/doubledescription.jl#L397 should have stopped add_intersection! because there is no ray to add here! Why did isin return false while r1 belongs to the hyperplane of the 3rd halfspace ? This is because when we checked it here: https://github.com/JuliaPolyhedra/Polyhedra.jl/blob/1cc19df35ea62efc3fefa162f709f04d09ea4c94/src/doubledescription.jl#L321, it returned -1.531648132849239e-7 which is not considered to be zero. But then, the scalar product because 4.0034642267983145e-13 after the ray r1 was line_projected with the line cut-off by the hyperplane of the 1st halfspace. This line_project should not alter whether the ray belongs to the 3rd hyperplane since the line should belong to it but the scalar product with the line is 1.486455444843493e-8 so... Changing the tolerance might help here but what we should really do is something more numerically stable like SVD-based for these cutoff lines

Zinoex commented 1 year ago

Thank you for taking the time to review this and update the show message. Since I was only interested in the normal vectors of the half-space constraints, I was able to center and rescale the V-representation to avoid the numerical issues that you have concluded. Changing the tolerances may have worked too. Again thank you for taking the time.