From bb9b051c5990ce214c84c960a76b05e2e619b886 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 11 Oct 2022 23:19:41 +0200 Subject: [PATCH] revise iterative_refinement code --- docs/src/lib/approximations.md | 11 ++- docs/src/man/iterative_refinement.md | 18 ++--- src/Approximations/Approximations.jl | 2 + src/Approximations/iterative_refinement.jl | 88 ++++++++++++---------- src/Approximations/overapproximate.jl | 4 +- 5 files changed, 64 insertions(+), 59 deletions(-) diff --git a/docs/src/lib/approximations.md b/docs/src/lib/approximations.md index 8f280839f9..404d62836c 100644 --- a/docs/src/lib/approximations.md +++ b/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 diff --git a/docs/src/man/iterative_refinement.md b/docs/src/man/iterative_refinement.md index 3eac4b315a..ccf1f29736 100644 --- a/docs/src/man/iterative_refinement.md +++ b/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. diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 3203f8bf9d..687eb395fe 100644 --- a/src/Approximations/Approximations.jl +++ b/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 diff --git a/src/Approximations/iterative_refinement.jl b/src/Approximations/iterative_refinement.jl index b5e7d6e72c..1ca5a39155 100644 --- a/src/Approximations/iterative_refinement.jl +++ b/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,9 +28,9 @@ 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 @@ -38,42 +38,42 @@ Convert a local approximation to a linear constraint. ### 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,8 +168,7 @@ 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 @@ -176,41 +176,47 @@ function tohrep(Ω::PolygonalOverapproximation) 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 diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 739c7941ad..799517b661 100644 --- a/src/Approximations/overapproximate.jl +++ b/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