No method matching difference(::VPolygon, ::VPolygon)

[ferrolho]

How can I subtract a set from another?

I am trying to use https://github.com/JuliaReach/LazySets.jl/blob/48974d76e8d8691bdfdd6d2a24dbbe37015149a8/src/ConcreteOperations/difference.jl#L38 as in the snippet

board = VPolygon([0.0 1.0 1.0 0.0
                  0.0 0.0 1.0 1.0])

shape₁ = VPolygon([0.2 0.5 0.7
                   0.2 0.9 0.1])

result = board \ shape₁

but get the following error:

ERROR: MethodError: no method matching difference(::VPolygon{Float64,Array{Float64,1}}, ::VPolygon{Float64,Array{Float64,1}})
Stacktrace:
 [1] \(::VPolygon{Float64,Array{Float64,1}}, ::VPolygon{Float64,Array{Float64,1}}) at /home/henrique/.julia/packages/LazySets/zvXxX/src/ConcreteOperations/difference.jl:38
 [2] top-level scope at REPL[7]:1

[mforets]

Let me remind that are two non-equivalent notions of set difference: Minkowski (or Pontryagin) difference on one hand, implemented in this library as minkowski_difference, and set difference implemented in this library as difference(X, Y) for which X \ Y is meant to be an alias and which is defined as the set of points which belong to X but not to Y. If you are interested in the latter seems like it is not implemented yet unfortunately. I don't remember that we have an open issue for it.

[mforets]

If your application is OK with an approximate result that's another story, it is not hard to overapproximate (cover) the shape with the set union of boxes, then compute the difference between the board and each element in the shape. In this case, we can estimate the error in the approximation as well and refine the approximation to the desired accuracy.

[ferrolho]

Thank you for your replies. I understand the difference, but I am looking for the latter. Picture poking a circular hole through a square piece of paper; that's the kind of operation I am looking for: result = square \ circle.

I'm having a look at Overapproximations right now. Is this what you expect to work?

[mforets]

I'm having a look at Overapproximations right now. Is this what you expect to work?

Yes exactly. I tried covering with boxes, the result is shown below and the code is available in this notebook (the implementation is not very clean due to some missing methods..). I dont understand why there are missing blue rectangles on the edges.

[ferrolho]

Thank you for taking the time to prepare a notebook. I will explore this workaround and see where it can get me.

I think I will leave this issue open for now; it can serve as a reminder that concrete difference of sets is not yet implemented.

[mforets]

How about using complements? If I didn't make any silly mistake (... it's late here :moon: ), the following formula applies:

Which is easy to implement:

complement(H::HalfSpace) = HalfSpace(-H.a, -H.b)

function _difference(X, Y)
    out = [intersection(X, complement(H)) for H in constraints_list(Y)]
    return UnionSetArray(out)
end

In fact I think the approach works with rather general assumptions on X and Y, such as Y being polyhedral.

The result is shown below, and I've updated the previous notebook Difference_VPolygon with the new function.

To hide the black lines delimiting each constraint we can use lw=0.0,

[schillic]

How about using complements?

Yes, this idea works as long as Y is a convex polyhedron. For more general convex Ys you can try to over- or underapproximate Y with such a polyhedron first (depending on whether you prefer an under- or overapproximation of the result).

EDIT: I added code and a plot to compute an underapproximation of the difference. We do not have good support to compute an overapproximation of the difference (which requires computing an underapproximation of Y), but we can think about it if you rather want that.

function _difference_ua(X, Y; ε=1e-3)
    Y_oa = overapproximate(Y, ε)
    return _difference(X, Y_oa)
end

shape₂ = Ball2([0.5, 0.5], 0.2)
result = _difference_ua(board, shape₂)

[ferrolho]

Thank you @mforets and @schillic for both of your interesting replies.

As suggested above, I have tried using the following snippet:

complement(H::HalfSpace) = HalfSpace(-H.a, -H.b)

function _difference(X, Y)
    out = [intersection(X, complement(H)) for H in constraints_list(Y)]
    return UnionSetArray(out)
end

It worked very well for a single convex shape. Below I have plotted a "before and after" subtracting a triangle from a square.

There are two issues that followed:

1. Trying to get the area results in the following error:

   ERROR: LoadError: MethodError: no method matching area(::UnionSetArray{Float64,HPolytope{Float64,Array{Float64,1}}})
   Closest candidates are:
   area(::EmptySet{N}) where N at /home/henrique/.julia/packages/LazySets/G2PtZ/src/Sets/EmptySet.jl:377
   area(::LazySet{N}) where N at /home/henrique/.julia/packages/LazySets/G2PtZ/src/Interfaces/LazySet.jl:947

   Is there a way I can convert the resulting UnionSetArray into a LazySet for the 2nd area method to be dispatched?

2. Indeed, I would like to subtract convex shapes, but not only one single shape. For example, in the figure below, I'd like to subtract all the triangles to the square. However, if I try to call _difference (from the code snippet above) for a UnionSetArray and a VPolygon, I get the following error:

   MethodError: no method matching UnionSetArray(::Array{UnionSetArray{Float64,S} where S<:LazySet{Float64},1})
   Closest candidates are:
   UnionSetArray(::Array{S,1}) where {N<:Real, S<:LazySet{N}} at /home/henrique/.julia/packages/LazySets/G2PtZ/src/LazyOperations/UnionSetArray.jl:18
   UnionSetArray(::EmptySet{N}, ::EmptySet{N}) where N<:Real at /home/henrique/.julia/packages/LazySets/G2PtZ/src/Utils/macros.jl:45
   UnionSetArray(::Universe{#283#N}, ::EmptySet{#283#N}) where #283#N<:Real at /home/henrique/.julia/packages/LazySets/G2PtZ/src/Utils/macros.jl:258
   ...
   Stacktrace:
   [1] _difference(::UnionSetArray{Float64,HPolytope{Float64,Array{Float64,1}}}, ::VPolygon{Float64,Array{Float64,1}}) at /home/henrique/tests/test_lazy_sets.jl:31
   [2] top-level scope at /home/henrique/tests/test_lazy_sets.jl:35
   [3] include(::String) at ./client.jl:439
   in expression starting at /home/henrique/tests/test_lazy_sets.jl:35

Edit: I will add the coordinates of the triangles below for reference.

board = VPolygon([0.0 1.0 1.0 0.0
                  0.0 0.0 1.0 1.0])

shape₁ = VPolygon([0.2 0.5 0.7
                   0.2 0.9 0.1])

shape₂ = VPolygon([0.8 0.1 0.2
                   0.8 0.4 0.8])

shape₃ = VPolygon([0.4 0.3 0.1
                   0.1 0.1 0.8])

shape₄ = VPolygon([0.7 0.9 0.9
                   0.6 0.4 0.6])

[schillic]

1. Computing the area of a union is tricky if the components overlay. I guess you would essentially have to do the work twice for that. If you are ultimately only interested in the area of the difference, I think it is better to instead compute the area of the intersection because that will be convex:

   area(board \ triangle) = area(board) - area(board ∩ triangle)

   (The philosophy to realize this idea in LazySets would be to define a lazy Difference type and then define area(::Difference) as given above.)

2. With multiple objects, if you go sequentially, after the first difference you are left with a non-convex set, and that complicates things considerably. I think if you are only interested in the area, you can again use a similar idea, but need to take care of overlays. Here is the code for two triangles.

   area(board \ (t1 ∪ t2)) = area(board) - area(board ∩ t1) - area(board ∩ t2) + area(board ∩ t1 ∩ t2)

   (The generalization to and arbitrary number of sets alternates - and +:

   area(board) - ∑_i area(board ∩ ti) + ∑_i≠j area(board ∩ ti ∩ tj) - ∑_i≠j≠k area(board ∩ ti ∩ tj ∩ tk) + ...

   An efficient implementation would need to track intersections.)