#326 - Template directions

src/Approximations/Approximations.jl

@@ -16,7 +16,9 @@ export approximate,
        norm,
        overapproximate,
        radius,
-       symmetric_interval_hull
+       symmetric_interval_hull,
+       BoxDirections,
+       OctDirections
 
 const TOL(N::Type{Float64}) = eps(N)
 const TOL(N::Type{Float32}) = eps(N)
@@ -29,6 +31,7 @@ const DIR_SOUTH(N) = [zero(N), -one(N)]
 
 include("iterative_refinement.jl")
 include("box_approximations.jl")
+include("template_directions.jl")
 include("overapproximate.jl")
 include("decompositions.jl")
 

src/Approximations/overapproximate.jl

@@ -121,6 +121,31 @@ function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
     return Zonotope(center, generators)
 end
 
+"""
+    overapproximate(X::LazySet, dir::AbstractDirections)::HPolytope
+
+Overapproximating a set with template directions.
+
+### Input
+
+- `X`   -- set
+- `dir` -- direction representation
+
+### Output
+
+A `HPolytope` overapproximating the set `X` with the directions from `dir`.
+"""
+function overapproximate(X::LazySet{N},
+                         dir::AbstractDirections{N})::HPolytope{N} where N
+    halfspaces = Vector{LazySets.LinearConstraint{N}}()
+    sizehint!(halfspaces, length(dir))
+    H = HPolytope(halfspaces)
+    for d in dir
+        addconstraint!(H, LazySets.LinearConstraint(d, ρ(d, X)))
+    end
+    return H
+end
+
 @require IntervalArithmetic begin
 
 """

src/Approximations/template_directions.jl

@@ -0,0 +1,143 @@
+import LazySets.dim
+
+"""
+    AbstractDirections{N}
+
+Abstract type for template direction representations.
+
+### Notes
+
+All subtypes should implement the standard iterator methods from `Base` and the
+function `dim(d<:AbstractDirections)::Int`.
+"""
+abstract type AbstractDirections{N} end
+
+# box directions
+
+"""
+    UnitVector{T} <: AbstractVector{T}
+
+A lazy unit vector with arbitrary one-element.
+
+### Fields
+
+- `i` -- index of non-zero entry
+- `n` -- vector length
+- `v` -- non-zero entry
+"""
+struct UnitVector{T} <: AbstractVector{T}
+    i::Int
+    n::Int
+    v::T
+end
+
+function Base.getindex(e::UnitVector{T}, i::Int) where T
+    @boundscheck @assert 1 <= i <= e.n
+    return i == e.i ? e.v : zero(T)
+end
+
+Base.size(e::UnitVector) = (e.n,)
+
+"""
+    BoxDirections{N} <: AbstractDirections{N}
+
+Box direction representation.
+
+### Fields
+
+- `n` -- dimension
+"""
+struct BoxDirections{N} <: AbstractDirections{N}
+    n::Int
+end
+
+# constructor for type Float64
+BoxDirections(n::Int) = BoxDirections{Float64}(n)
+
+Base.eltype(::Type{BoxDirections{N}}) where N = AbstractVector{N}
+Base.start(bd::BoxDirections) = 1
+Base.next(bd::BoxDirections{N}, state) where N = (
+    UnitVector{N}(abs(state), bd.n, convert(N, sign(state))), # value
+    state == bd.n ? -bd.n : state + 1) # next state
+Base.done(bd::BoxDirections, state) = state == 0
+Base.length(bd::BoxDirections) = 2*bd.n
+
+"""
+    dim(bd::BoxDirections)::Int
+
+Returns the dimension of the generated directions.
+
+### Input
+
+- `bd` -- box direction representation
+
+### Output
+
+The dimension of the generated directions.
+"""
+function dim(bd::BoxDirections)::Int
+    return bd.n
+end
+
+# octagon directions
+
+"""
+    OctDirections{N} <: AbstractDirections{N}
+
+Octagon direction representation.
+
+### Fields
+
+- `n` -- dimension
+
+### Notes
+
+Octagon directions can be seen as the union of diagonal directions (all entries
+are ±1) and box directions (one entry is ±1, all other entries are 0).
+The iterator first enumerates all diagonal directions; then it enumerates all
+box directions.
+"""
+struct OctDirections{N} <: AbstractDirections{N}
+    n::Int
+end
+
+# constructor for type Float64
+OctDirections(n::Int) = OctDirections{Float64}(n)
+
+Base.eltype(::Type{OctDirections{N}}) where N = AbstractVector{N}
+Base.start(od::OctDirections{N}) where N = ones(N, od.n)
+function Base.next(od::OctDirections{N}, state::AbstractVector) where N
+    i = 1
+    while i <= od.n && state[i] < 0
+        state[i] = -state[i]
+        i = i+1
+    end
+    if i > od.n
+        return (copy(state), 1)
+    else
+        state[i] = -state[i]
+        return (copy(state), state)
+    end
+end
+function Base.next(od::OctDirections{N}, state::Int) where N
+    return next(BoxDirections{N}(od.n), state)
+end
+Base.done(od::OctDirections, state) = state == 0
+Base.length(od::OctDirections) = 2^od.n + 2*od.n
+
+"""
+    dim(od::OctDirections)::Int
+
+Returns the dimension of the generated directions.
+
+### Input
+
+- `od` -- octagon direction representation
+
+### Output
+
+The dimension of the generated directions.
+"""
+function dim(od::OctDirections)::Int
+    return od.n
+end

test/runtests.jl

@@ -37,6 +37,7 @@ include("to_N.jl")
 # Algorithms for approximation of convex sets using support vectors
 # =================================================================
 @time @testset "LazySets.Approximations.overapproximation" begin include("unit_overapproximate.jl") end
+@time @testset "LazySets.Approximations.template_directions" begin include("unit_template_directions.jl") end
 @time @testset "LazySets.Approximations.box_approximation" begin include("unit_box_approximation.jl") end
 @time @testset "LazySets.Approximations.ballinf_approximation" begin include("unit_ballinf_approximation.jl") end
 @time @testset "LazySets.Approximations.radiusdiameter" begin include("unit_radiusdiameter.jl") end

test/unit_template_directions.jl

@@ -0,0 +1,36 @@
+import LazySets.Approximations,
+       Approximations.UnitVector,
+       Approximations.BoxDirections,
+       Approximations.OctDirections,
+       Approximations.overapproximate
+
+for N in [Float64, Float32, Rational{Int}]
+    # unit vector
+    for n in 1:3
+        for i in 1:n
+            e = UnitVector(i, n, one(N))
+            for j in 1:n
+                @test e[j] == (i == j ? one(N) : zero(N))
+            end
+        end
+    end
+
+    # template direction approximation
+    for n in 2:3
+        B = BallInf(zeros(N, n), N(2.))
+        A = 2.*eye(N, n) + ones(N, n, n)
+        X = A * B
+
+        # box directions
+        dir = BoxDirections{N}(n)
+        @test dim(dir) == n
+        box = overapproximate(X, dir)
+        @test length(dir) == length(box.constraints) == 2*n
+
+        # octagon directions
+        dir = OctDirections{N}(n)
+        @test dim(dir) == n
+        oct = overapproximate(X, dir)
+        @test length(dir) == length(oct.constraints) == 2^n + 2*n
+    end
+end