#375 - Intersection of a halfspace with polyhedron: line search method

[mforets]

Closes #375.

[mforets]

Some examples:

using LazySets, Optim, Compat
p = Ball2(zeros(2), 1.0)
h = HalfSpace([1.0, 1.0], 0.0) # x + y <= 0

plot(p, aspectratio=1, alpha=0.2, 1e-3)
plot!(p ∩ h, 1e-3, aspectratio=1, alpha=0.2, color="red")

[mforets]

This one seems wrong, but i don't know why:

p = VPolygon([rand(2) for i in 1:8])
plot(p, aspectratio=1, alpha=0.2, 1e-3)
h = HalfSpace([1.0, 0.0], 0.5) # x <= 0.5
plot!(p ∩ h, 1e-3, aspectratio=1, alpha=0.2, color="red")

[schillic]

How come that σ is defined for you? I can only plot the box approximation.

[mforets]

sorry, i forgot to push that part. just a second.

[mforets]

do you agree about my definition of σ?

[mforets]

Hmm it is not a problem with VPolygon. Here is a similar error (i only get the overapproximation) using a Ball2:

p = Ball2(ones(2), 1.0)
plot(p, aspectratio=1, alpha=0.2, 1e-3)
h = HalfSpace([1.0, 0.0], 1.0) # x <= 1
plot!(p ∩ h, 1e-3, aspectratio=1, alpha=0.2, color="red")

[mforets]

Here i get

julia> σ([1., 1], p∩h)
2-element Array{Float64,1}:
 1.5
 1.5

but i think it should be [1.0, 2.0].

Also

julia> σ([-1., 1], p∩h)
2-element Array{Float64,1}:
 -0.707107
  0.707107

gives a point which does not belong to p.

[mforets]

ok the problem is that σ does not work like ρ(ℓvec, cap) * ℓvec in general, since this condition is not sufficient for the point to belong to the set p. ...

for instance to test that the support function is correct, consider the same example but using template directions (which evaluates the supp func) instead of iterative refinement (which evaluates the supp vector):

using LazySets.Approximations

p = Ball2(ones(2), 1.0)
h = HalfSpace([1.0, 0.0], 1.0) # x <= 1
pint = overapproximate(p∩h, BoxDiagDirections(2))
pint = convert(HPolygon, pint)

plot(p, aspectratio=1, alpha=0.2, 1e-3)
plot!(pint, aspectratio=1, alpha=0.2, color="red")

[schillic]

ok the problem is that σ does not work like ρ(ℓvec, cap) * ℓvec in general, since this condition is not sufficient for the point to belong to the set p. ...

Yes, this does not work in general.

julia> b = BallInf([0., 0.], 1.); d = [1., 0.5];

julia> ρ(d, b)
1.5

julia> ρ(d, b) * d
2-element Array{Float64,1}:
 1.5 
 0.75

julia> ρ(d, b) * d ∈ b
false

julia> σ(d, b)
2-element Array{Float64,1}:
 1.0
 1.0

[mforets]

Yes, thanks for checking. I don't have a formula for the support vector, so i propose to move on with the computation of supp func for now.

It remains to add a add a unit test and the function to the docs. Regarding the Optim package, shall we do something similar to what is done in function load_polyhedra() ? Where do i get the number in

function __init__()
    @require Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029" load_polyhedra()
    @require Optim = "????" load_optim()
end

The theorem on f(λ) can be found in Rockafellar's Variational Anaysis. I'd like to check there, i suspect one can extract the point in the set from the minimizer..

[schillic]

A preliminary idea for the support vector σ(d, p ∩ h): Let x := σ(d, p). Either x ∈ h holds, in which case you can return x, or we know that σ(d, p ∩ h) lies on the hyperplane described by d and ρ(d, p). So you can shift x along that hyperplane until you hit the hyperplane of h. Moving along a hyperplane is not unique, though...

[schillic]

Where do i get the number

I look it up in .julia/environments/v0.7.

[schillic]

Do we use "ℓ" for directions now? (I am not against that choice.)

I would also add an option to skip the emptiness check. In practice we want to do that in advance before computing ρ. The name of that option is tricky: A possible name is_known_to_be_nonempty is a bit overkill. Maybe skip_emptiness_check?

[mforets]

Do we use "ℓ" for directions now? (I am not against that choice.)

i was using it because i was following the notation from the paper. but i shall switch back to d.

I would also add an option to skip the emptiness check.

good idea. maybe emptiness_check.

[schillic]

Here is my late review 🐌

[schillic]

The check_intersection argument is missing here.

[schillic]

new line

[schillic]

Should the second expression be the symmetric case?

[mforets]

Thanks for the comments. I put them in #681