Merge pull request #3121 from JuliaReach/schillic/iterative_refinement

Revise iterative_refinement code

docs/src/lib/approximations.md

@@ -61,15 +61,14 @@ ballinf_approximation
 ## Iterative refinement
 
 ```@docs
+overapproximate_hausdorff
 LocalApproximation
 PolygonalOverapproximation
-new_approx(S::ConvexSet, p1::VN, d1::VN,
-           p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}
-addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN, p2::VN, d2::VN) where {N<:Real, VN<:AbstractVector{N}}
-refine(::LocalApproximation, ::ConvexSet)
+new_approx
+addapproximation!
+refine(::LocalApproximation, ::LazySet)
 tohrep(::PolygonalOverapproximation)
-_approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}
-constraint(::LocalApproximation)
+convert(::Type{HalfSpace}, ::LocalApproximation)
 ```
 
 ## Template directions

docs/src/man/iterative_refinement.md

@@ -68,7 +68,7 @@ Hausdorff distance.
 
 
 ```@example example_iterative_refinement
-import LazySets.Approximations:PolygonalOverapproximation, addapproximation!
+using LazySets.Approximations: PolygonalOverapproximation, addapproximation!
 
 Ω = PolygonalOverapproximation(b)
 p1, d1, p2, d2 = [1.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.0, 1.0]
@@ -98,7 +98,7 @@ To illustrate this, first let's add the remaining three approximations to `Ω`
 along the canonical directions, to build a box overapproximation of `b`.
 
 ```@example example_iterative_refinement
-import LazySets.Approximations: refine, tohrep
+using LazySets.Approximations: refine, tohrep
 
 plot(b, 1e-3, aspectratio=1, alpha=0.3)
 
@@ -161,7 +161,7 @@ point epsilon in the given numerical precision.
 ## Algorithm
 
 Having presented the individual steps, we give the pseudocode of the iterative
-refinement algorithm, see `_approximate(S, ε)`.
+refinement algorithm, see `overapproximate_hausdorff(S, ε)`.
 
 The algorithm consists of the following steps:
 
@@ -187,7 +187,7 @@ As a final example consider the iterative refinement of the ball `b` for
 different values of the approximation threshold `ε`.
 
 ```@example example_iterative_refinement
-import LazySets.Approximations:overapproximate, _approximate
+using LazySets.Approximations: overapproximate_hausdorff
 
 p0 = plot(b, 1e-6, aspectratio=1)
 p1 = plot!(p0, overapproximate(b, 1.), alpha=0.4, aspectratio=1)
@@ -205,15 +205,15 @@ Meanwhile, the number of constraints of the polygonal overapproximation
 increases, in this example by a power of 2 when the error is divided by a factor 10.
 
 ```@example example_iterative_refinement
-h = ε ->  length(_approximate(b, ε).constraints)
+h = ε ->  length(overapproximate_hausdorff(b, ε).constraints)
 h(1.), h(0.1), h(0.01)
 ```
 
 !!! note
-    Actually, the plotting function for an arbitrary `ConvexSet` `plot(...)`,
-    called *recipe* in the context of
-    [Plots.jl](https://github.com/JuliaPlots/Plots.jl), is such that it receives
-    a numeric argument `ε` and the routine itself calls `overapproximate`.
+    Actually, the plotting function for an arbitrary convex `LazySet`,
+    `plot(...)` (called *recipe* in the context of
+    [Plots.jl](https://github.com/JuliaPlots/Plots.jl)), receives a numeric
+    argument `ε` and the routine itself calls `overapproximate`.
     However, some sets such as abstract polygons have their own plotting recipe
     and hence do not require the error threshold, since they are plotted exactly
     as the convex hull of their vertices.

src/Approximations/Approximations.jl

@@ -5,6 +5,8 @@ Module `Approximations.jl` -- polygonal approximation of sets.
 """
 module Approximations
 
+import Base: convert
+
 import IntervalArithmetic
 const IA = IntervalArithmetic
 

src/Approximations/iterative_refinement.jl

@@ -10,12 +10,12 @@ Type that represents a local approximation in 2D.
 - `p2`        -- second inner point
 - `d2`        -- second direction
 - `q`         -- intersection of the lines l1 ⟂ d1 at p1 and l2 ⟂ d2 at p2
-- `refinable` -- states if this approximation is refinable
+- `refinable` -- flag stating whether this approximation is refinable
 - `err`       -- error upper bound
 
 ### Notes
 
-The criteria for being refinable are determined in the method `new_approx`.
+The criteria for being refinable are determined in [`new_approx`](@ref).
 """
 struct LocalApproximation{N, VN<:AbstractVector{N}}
     p1::VN
@@ -28,52 +28,52 @@ struct LocalApproximation{N, VN<:AbstractVector{N}}
 end
 
 """
-    constraint(approx::LocalApproximation)
+    convert(::Type{HalfSpace}, approx::LocalApproximation)
 
-Convert a local approximation to a linear constraint.
+Convert a local approximation to a half-space.
 
 ### Input
 
 - `approx` -- local approximation
 
 ### Output
 
-A linear constraint.
+A half-space.
 """
-function constraint(approx::LocalApproximation)
-    return LinearConstraint(approx.d1, dot(approx.d1, approx.p1))
+function convert(::Type{HalfSpace}, approx::LocalApproximation)
+    return HalfSpace(approx.d1, dot(approx.d1, approx.p1))
 end
 
 """
-    PolygonalOverapproximation{N, SN<:ConvexSet{N}, VN<:AbstractVector{N}}
+    PolygonalOverapproximation{N, SN<:LazySet{N}, VN<:AbstractVector{N}}
 
-Type that represents the polygonal approximation of a convex set.
+Type that represents a polygonal overapproximation of a convex set.
 
 ### Fields
 
 - `S`            -- convex set
 - `approx_stack` -- stack of local approximations that still need to be examined
-- `constraints`  -- vector of linear constraints that are already finalized
+- `constraints`  -- vector of half-spaces that are already finalized
                     (i.e., they satisfy the given error bound)
 """
-struct PolygonalOverapproximation{N, SN<:ConvexSet{N}, VN<:AbstractVector{N}}
+struct PolygonalOverapproximation{N, SN<:LazySet{N}, VN<:AbstractVector{N}}
     S::SN
     approx_stack::Vector{LocalApproximation{N, VN}}
-    constraints::Vector{LinearConstraint{N, VN}}
+    constraints::Vector{HalfSpace{N, VN}}
 end
 
-function PolygonalOverapproximation(S::SN) where {N, SN<:ConvexSet{N}}
+function PolygonalOverapproximation(S::SN) where {N, SN<:LazySet{N}}
     empty_local_approx = Vector{LocalApproximation{N, Vector{N}}}()
-    empty_constraints = Vector{LinearConstraint{N,Vector{N}}}()
+    empty_constraints = Vector{HalfSpace{N,Vector{N}}}()
     return PolygonalOverapproximation(S, empty_local_approx, empty_constraints)
 end
 
 """
-    new_approx(S::ConvexSet, p1::VN, d1::VN,
+    new_approx(S::LazySet, p1::VN, d1::VN,
                p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}
 
 Create a `LocalApproximation` instance for the given excerpt of a polygonal
-approximation.
+overapproximation.
 
 ### Input
 
@@ -87,7 +87,7 @@ approximation.
 
 A local approximation of `S` in the given directions.
 """
-function new_approx(S::ConvexSet, p1::VN, d1::VN,
+function new_approx(S::LazySet, p1::VN, d1::VN,
                     p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}
     if norm(p1-p2, 2) <= _rtol(N)
         # this approximation cannot be refined and we set q = p1 by convention
@@ -129,10 +129,10 @@ function addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN,
 end
 
 """
-    refine(approx::LocalApproximation, S::ConvexSet)
+    refine(approx::LocalApproximation, S::LazySet)
 
-Refine a given local approximation of the polygonal approximation of a convex
-set by splitting along the normal direction of the approximation.
+Refine a given local approximation of the polygonal overapproximation of a
+convex set by splitting along the normal direction of the approximation.
 
 ### Input
 
@@ -143,7 +143,7 @@ set by splitting along the normal direction of the approximation.
 
 The tuple consisting of the refined right and left local approximations.
 """
-function refine(approx::LocalApproximation, S::ConvexSet)
+function refine(approx::LocalApproximation, S::LazySet)
     @assert approx.refinable
     ndir = normalize([approx.p2[2]-approx.p1[2], approx.p1[1]-approx.p2[1]])
     s = σ(ndir, S)
@@ -155,7 +155,8 @@ end
 """
     tohrep(Ω::PolygonalOverapproximation)
 
-Convert a polygonal overapproximation into a concrete polygon.
+Convert a polygonal overapproximation into a polygon in constraint
+representation.
 
 ### Input
 
@@ -167,50 +168,55 @@ A polygon in constraint representation.
 
 ### Algorithm
 
-Internally we keep the constraints sorted.
-Hence we do not need to use `addconstraint!` when creating the `HPolygon`.
+Internally, the constraints of `Ω` are already sorted.
 """
 function tohrep(Ω::PolygonalOverapproximation)
     # already finalized
     if isempty(Ω.approx_stack)
         return HPolygon(Ω.constraints, sort_constraints=false)
     end
     # some constraints not finalized yet
-    constraints = Ω.constraints
+    P = HPolygon(Ω.constraints, sort_constraints=false)
     for approx in Ω.approx_stack
-        push!(constraints, constraint(approx))
+        # the remaining constraints need to be sorted
+        addconstraint!(P, convert(HalfSpace, approx))
     end
-    return HPolygon(constraints, sort_constraints=false)
+    return P
 end
 
 """
-    _approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}
+    overapproximate_hausdorff(X::S, ε::Real) where {N<:AbstractFloat, S<:LazySet{N}}
 
-Return an ε-close approximation of the given 2D convex set (in terms of
-Hausdorff distance) as an inner and an outer approximation composed by sorted
-local `Approximation2D`.
+Return an ε-close overapproximation of the given 2D convex set (in terms of the
+Hausdorff distance) in the form of a polygon in constraint representation.
 
 ### Input
 
-- `S` -- 2D convex set
+- `X` -- 2D convex set
 - `ε` -- error bound
 
 ### Output
 
-An ε-close approximation of the given 2D convex set.
+A polygon in constraint representation.
 """
-function _approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}
+function overapproximate_hausdorff(X::S, ε::Real) where {N<:AbstractFloat, S<:LazySet{N}}
+    if !isconvextype(S)
+        error("ε-close overapproximation requires a convex set")
+    elseif dim(X) != 2
+        error("ε-close overapproximation requires a two-dimensional set")
+    end
+
     # initialize box directions
-    pe = σ(DIR_EAST(N), S)
-    pn = σ(DIR_NORTH(N), S)
-    pw = σ(DIR_WEST(N), S)
-    ps = σ(DIR_SOUTH(N), S)
+    pe = σ(DIR_EAST(N), X)
+    pn = σ(DIR_NORTH(N), X)
+    pw = σ(DIR_WEST(N), X)
+    ps = σ(DIR_SOUTH(N), X)
     east = dir_east(N, pe)
     north = dir_north(N, pe)
     west = dir_west(N, pe)
     south = dir_south(N, pe)
 
-    Ω = PolygonalOverapproximation(S)
+    Ω = PolygonalOverapproximation(X)
 
     # add constraints in reverse (i.e., clockwise) order to the stack
     addapproximation!(Ω, ps, south, pe, east)
@@ -224,7 +230,7 @@ function _approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}
 
         if !approx.refinable || approx.err <= ε
             # if the approximation is not refinable => continue
-            push!(Ω.constraints, constraint(approx))
+            push!(Ω.constraints, convert(HalfSpace, approx))
             continue
         end
 
@@ -248,5 +254,5 @@ function _approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}
         end
         push!(approx_stack, la1)
     end
-    return Ω
+    return tohrep(Ω)
 end

src/Approximations/overapproximate.jl

@@ -1,5 +1,3 @@
-import Base: convert
-
 using LazySets: block_to_dimension_indices,
                 substitute_blocks,
                 fast_interval_pow,
@@ -90,7 +88,7 @@ function overapproximate(S::LazySet{N},
         constraints[4] = LinearConstraint(DIR_SOUTH(N), ρ(DIR_SOUTH(N), S))
         return HPolygon(constraints, sort_constraints=false)
     else
-        P = tohrep(_approximate(S, ε))
+        P = overapproximate_hausdorff(S, ε)
         if prune
             remove_redundant_constraints!(P)
         end