Outsource Tetrahedron to its own module

docs/src/lib/sets/SparsePolynomialZonotope.md

@@ -6,6 +6,11 @@ CurrentModule = LazySets.SparsePolynomialZonotopeModule
 
 ```@docs
 SparsePolynomialZonotope
+```
+
+## Operations
+
+```@docs
 center(::SparsePolynomialZonotope)
 expmat(::SparsePolynomialZonotope)
 genmat_dep(::SparsePolynomialZonotope)

docs/src/lib/sets/Tetrahedron.md

@@ -1,11 +1,16 @@
 ```@meta
-CurrentModule = LazySets
+CurrentModule = LazySets.TetrahedronModule
 ```
 
 # [Tetrahedron](@id def_Tetrahedron)
 
 ```@docs
 Tetrahedron
-σ(::AbstractVector, ::Tetrahedron)
+```
+
+## Operations
+
+```@docs
 ∈(::AbstractVector, ::Tetrahedron)
+σ(::AbstractVector, ::Tetrahedron)
 ```

src/LazySets.jl

@@ -164,7 +164,9 @@ const PolynomialZonotope = SimpleSparsePolynomialZonotope
 
 include("Sets/VPolygon.jl")
 include("Sets/VPolytope.jl")
-include("Sets/Tetrahedron.jl")
+
+include("Sets/Tetrahedron/TetrahedronModule.jl")
+@reexport using ..TetrahedronModule: Tetrahedron
 
 include("Sets/ZeroSet/ZeroSetModule.jl")
 @reexport using ..ZeroSetModule: ZeroSet

src/Sets/Tetrahedron/Tetrahedron.jl

@@ -0,0 +1,46 @@
+"""
+    Tetrahedron{N, VN<:AbstractVector{N}, VT<:AbstractVector{VN}} <: AbstractPolytope{N}
+
+Type that represents a tetrahedron in vertex representation.
+
+### Fields
+
+- `vertices` -- list of vertices
+
+### Examples
+
+A tetrahedron can be constructed by passing the list of vertices.
+The following builds the tetrahedron with edge length 2 from the [wikipedia page Tetrahedron](https://en.wikipedia.org/wiki/Tetrahedron):
+
+```jldoctest
+julia> vertices = [[1, 0, -1/sqrt(2)], [-1, 0, -1/sqrt(2)], [0, 1, 1/sqrt(2)], [0, -1, 1/sqrt(2)]];
+
+julia> T = Tetrahedron(vertices);
+
+julia> dim(T)
+3
+
+julia> zeros(3) ∈ T
+true
+
+julia> σ(ones(3), T)
+3-element Vector{Float64}:
+ 0.0
+ 1.0
+ 0.7071067811865475
+```
+"""
+struct Tetrahedron{N,VN<:AbstractVector{N},VT<:AbstractVector{VN}} <: AbstractPolytope{N}
+    vertices::VT
+
+    function Tetrahedron(vertices::VT) where {N,VN<:AbstractVector{N},VT<:AbstractVector{VN}}
+        @assert length(vertices) == 4 "a tetrahedron requires four vertices"
+        return new{N,VN,VT}(vertices)
+    end
+end
+
+# constructor from rectangular matrix
+function Tetrahedron(vertices_matrix::MT) where {N,MT<:AbstractMatrix{N}}
+    vertices = collect(eachcol(vertices_matrix))
+    return Tetrahedron(vertices)
+end

src/Sets/Tetrahedron/TetrahedronModule.jl

@@ -0,0 +1,24 @@
+module TetrahedronModule
+
+using Reexport
+
+using ..LazySets: AbstractPolytope, VPolytope
+using Random: AbstractRNG, GLOBAL_RNG
+using ReachabilityBase.Comparison: isapproxzero
+using LinearAlgebra: dot, cross
+
+@reexport import ..API: constraints_list, dim, isoperationtype, rand, ∈, σ
+@reexport using ..API
+
+export Tetrahedron
+
+include("Tetrahedron.jl")
+
+include("constraints_list.jl")
+include("dim.jl")
+include("isoperationtype.jl")
+include("rand.jl")
+include("in.jl")
+include("support_vector.jl")
+
+end  # module

src/Sets/Tetrahedron/constraints_list.jl

@@ -0,0 +1,3 @@
+function constraints_list(T::Tetrahedron)
+    return constraints_list(convert(VPolytope, T))
+end

src/Sets/Tetrahedron/dim.jl

@@ -0,0 +1,3 @@
+function dim(::Tetrahedron)
+    return 3
+end

src/Sets/Tetrahedron/in.jl

@@ -0,0 +1,40 @@
+"""
+    ∈(x::AbstractVector, T::Tetrahedron)
+
+Check whether a given point is contained in a tetrahedron.
+
+### Input
+
+- `x` -- point/vector
+- `P` -- tetrahedron in vertex representation
+
+### Output
+
+`true` iff ``x ∈ P``.
+
+### Algorithm
+
+For each plane of the tetrahedron, we check if the point `x` is on the same side as the remaining vertex.
+We need to check this for each plane [1].
+
+[1] https://stackoverflow.com/questions/25179693/how-to-check-whether-the-point-is-in-the-tetrahedron-or-not
+"""
+function ∈(x::AbstractVector, T::Tetrahedron)
+    v = T.vertices
+    return same_side(v[1], v[2], v[3], v[4], x) &&
+           same_side(v[4], v[1], v[2], v[3], x) &&
+           same_side(v[3], v[4], v[1], v[2], x) &&
+           same_side(v[2], v[3], v[4], v[1], x)
+end
+
+# Return `true` iff point `x` lies in the same half-space as `v4` with respect to the hyperplane `H` determined by `v1`, `v2` and `v3`.
+function same_side(v1, v2, v3, v4, x)
+    normal = cross(v2 - v1, v3 - v1)
+    dotx = dot(normal, x - v1)
+    if isapproxzero(dotx)
+        return true
+    end
+    dotv4 = dot(normal, v4 - v1)
+    # If the point `x` lies in `H` then `dotx` is zero.
+    return signbit(dotv4) == signbit(dotx)
+end

src/Sets/Tetrahedron/isoperationtype.jl

@@ -0,0 +1 @@
+isoperationtype(::Type{<:Tetrahedron}) = false

src/Sets/Tetrahedron/rand.jl

@@ -0,0 +1,6 @@
+function rand(::Type{Tetrahedron}; N::Type{<:Real}=Float64, rng::AbstractRNG=GLOBAL_RNG,
+              seed::Union{Int,Nothing}=nothing)
+    P = rand(VPolytope; N=N, dim=3, rng=rng, seed=seed, num_vertices=4)
+    vertices = P.vertices
+    return Tetrahedron(vertices)
+end

src/Sets/Tetrahedron/support_vector.jl

@@ -0,0 +1,21 @@
+"""
+    σ(d::AbstractVector, P::VPolytope)
+
+Return a support vector of a tetrahedron in a given direction.
+
+### Input
+
+- `d` -- direction
+- `P` -- tetrahedron
+
+### Output
+
+A support vector in the given direction.
+
+### Algorithm
+
+Currently falls back to the `VPolytope` implementation.
+"""
+function σ(d::AbstractVector, T::Tetrahedron)
+    return σ(d, convert(VPolytope, T))
+end

test/Sets/Tetrahedron.jl

@@ -10,6 +10,10 @@ for N in [Float64, Float32, Rational{Int}]
     Tmat = Tetrahedron(vmat)
     @test Tmat == T
 
+    # constraints_list
+    c = constraints_list(T)
+    @test isequivalent(HPolytope(c), VPolytope(vertices))
+
     # support vector
     @test σ(ones(3), T) == N[0, 1, 1 / sqrt(2)]
 
@@ -23,3 +27,10 @@ for N in [Float64, Float32, Rational{Int}]
     T = Tetrahedron([N[0, 0, 0], N[0, 1, 0], N[1, 0, 0], N[1, 0, 1]])
     @test zeros(N, 3) ∈ T
 end
+
+for N in [Float64, Float32]
+    # rand
+    @test rand(Tetrahedron; N=N) isa Tetrahedron{N}
+end
+
+@test !isoperationtype(Tetrahedron)