==/issubset on HPolytopes

[tomerarnon]

Several behaviors observed:

• == returning false negatives (fixed by defines == generically to work on LazySet pairs #604)
• issubset returning false positives

== (fixed by #604)

x = HPolytope(rand(2,2), rand(2))
y = deepcopy(x)

x == y  # returns false

x.constraints[1] == y.constraints[1] # also returns false, which suggests HalfSpaces also do not have a comparison operator

Defining:

 ==(h1::HalfSpace, h2::HalfSpace) = h1.a == h2.a && h1.b == h2.b
 ==(p1::HPolytope, p2::HPolytope) = p1.constraints == p2.constraints

solves this. I am not familiar enough with the package to know whether this definition is problematic, or whether a more general one is possible, etc.

issubset

Also, given two HPolytopes (these two were generated using rand() values and happened to exhibit this behavior. Not all pairs of HPolytopes do this, and I don't yet know why)

julia> x1 = HPolytope{Float64}(HalfSpace[HalfSpace([0.972936, 0.696308], 0.991854), HalfSpace([0.739235, 0.788761], 0.209161)]);
julia> x2 = HPolytope{Float64}(HalfSpace[HalfSpace([5.81193, 8.24438], 0.388707), HalfSpace([1.0198, 3.2731], 0.303999)]);

julia> issubset(x1, x2)
true
julia> issubset(x2, x1)
true

x1 and x2 cannot both be subsets of one another, since they are not equal.

I cannot find the general reason for this. For some values of x1 and x2, issubset(x1, x2) works but

julia> issubset(x2, x1)
ERROR: AssertionError: ispoint == zero(T) || ispoint == one(T)
Stacktrace:
 [1] isrowpoint(::Ptr{CDDLib.Cdd_MatrixData{Float64}}, ::Int64, ::Type{Float64}) at /Users/tomer/.julia/v0.6/CDDLib/src/matrix.jl:234
 [2] CDDLib.CDDGeneratorMatrix{2,Float64,Float64}(::Ptr{CDDLib.Cdd_MatrixData{Float64}}) at /Users/tomer/.julia/v0.6/CDDLib/src/matrix.jl:116
 [3] copygenerators(::CDDLib.CDDPolyhedra{2,Float64,Float64}) at /Users/tomer/.julia/v0.6/CDDLib/src/polyhedra.jl:109
 [4] getext(::CDDLib.CDDPolyhedron{2,Float64}) at /Users/tomer/.julia/v0.6/CDDLib/src/polyhedron.jl:52
 [5] removevredundancy!(::CDDLib.CDDPolyhedron{2,Float64}) at /Users/tomer/.julia/v0.6/CDDLib/src/polyhedron.jl:244
 [6] #vertices_list#81(::CDDLib.CDDLibrary, ::Polyhedra.#removevredundancy!, ::Function, ::LazySets.HPolytope{Float64}) at /Users/tomer/.julia/v0.6/LazySets/src/HPolytope.jl:370
 [7] issubset(::LazySets.HPolytope{Float64}, ::LazySets.HPolytope{Float64}, ::Bool) at /Users/tomer/.julia/v0.6/LazySets/src/is_subset.jl:182
 [8] issubset(::LazySets.HPolytope{Float64}, ::LazySets.HPolytope{Float64}) at /Users/tomer/.julia/v0.6/LazySets/src/is_subset.jl:180

[schillic]

About ==:
The comparison with == is not implemented yet (see #378). The default implementation falls back to ===, which is object identity. Let us know if you need this feature, then we can prioritize.

[tomerarnon]

No, == is not terribly critical for me. I was using it to test other features to see that they worked as I expected. I.e. performing some transformations and then comparing the result to the one I expected. Surprising results led to visual inspection led to discovering the == issue. Unfortunately I searched the open issues for "isequal" rather than "=="

However, given the recommendation in the docs that is mentioned in your linked issue, I wonder if the following definition could be useful as a catchall:

function ==(a::L, b::L) where L<:LazySet 
    for field in fieldnames(a)
        getfield(a, field) == getfield(b, field) || return false
    end
    return true
end

It is certainly less efficient than a per-type definition that leverages knowledge of the type, but would easily fill the void until the complete set of == is supported. This does assume that all fields internal to all LazySets are either other LazySets or types that already have == defined on them correctly.

[mforets]

However, given the recommendation in the docs that is mentioned in your linked issue, I wonder if the following definition could be useful as a catchall:

Looks good to me as well. This function could be in LazySet.jl and needs one to add import Base.== on top of that file. Let us know if you plan to send a PR with this new definition.

[tomerarnon]

Sure thing. Submitting the pull request for that now 😄

[schillic]

About issubset:
The implementation for two polytopic sets X ⊆ Y checks that each vertex of X is contained in Y. The problem with your example is that a polytope with just two half-spaces is unbounded. In this case the implementation only returns a single vertex, which happens to satisfy all half-space constraints here.

Hence I would say it is a misuse here. For efficiency I would not require that we check for boundedness on our side. We could provide a function that checks that explicitly (I created #605).

[tomerarnon]

Ah I see. I hadn't considered that the set was not bounded. Thanks for the explanation.

I agree that checking boundedness in issubset is overkill, but I am not sure that, in general, comparing two unbounded polytopes is necessarily a misuse. That said, I don't have any reasonable solution in mind for dealing with this.

[mforets]

I agree that it makes sense to define set inclusion for unbounded polyhedra. For that use case i would go for adding a new type Polyhedron, and keep the "polytopes" types to be bounded by assumption. We can define specialized functions under this hypothesis.

[schillic]

I created #606.