Merge pull request #2907 from JuliaReach/schillic/1312

#1312 - Add more exponentiation backends

Author: schillic

Project.toml

@@ -24,6 +24,7 @@ CDDLib = "0.6, 0.7"
 Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 Documenter = "0.23, 0.24, 0.25, 0.26, 0.27"
 Expokit = "0.2"
+ExponentialUtilities = "1"
 ExprTools = "0.1"
 GLPK = "0.11, 0.12, 0.13, 0.14, 0.15"
 GLPKMathProgInterface = "0.4, 0.5"
@@ -52,6 +53,7 @@ CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4"
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
 IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
@@ -68,4 +70,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [targets]
-test = ["CDDLib", "Distributions", "Documenter", "Expokit", "GR", "IntervalConstraintProgramming", "IntervalMatrices", "Makie", "Optim", "Polyhedra", "RecipesBase", "StaticArrays", "Symbolics", "TaylorModels", "Test"]
+test = ["CDDLib", "Distributions", "Documenter", "Expokit", "ExponentialUtilities", "GR", "IntervalConstraintProgramming", "IntervalMatrices", "Makie", "Optim", "Polyhedra", "RecipesBase", "StaticArrays", "Symbolics", "TaylorModels", "Test"]

docs/src/lib/utils.md

@@ -54,6 +54,7 @@ _leq_trig
 _random_zero_sum_vector
 rectify
 require(::Symbol)
+require(::Vector{Symbol})
 reseed
 same_block_structure
 substitute

src/Initialization/init_ExponentialUtilities.jl

@@ -0,0 +1 @@
+eval(load_exponentialutilities())

src/LazyOperations/ExponentialMap.jl

@@ -10,14 +10,6 @@ export SparseMatrixExp,
        get_column,
        get_columns
 
-function load_expokit()
-return quote
-
-using .Expokit: expmv
-
-end end  # quote / load_expokit
-
-
 # --- SparseMatrixExp & ExponentialMap ---
 
 """
@@ -39,7 +31,7 @@ julia> using SparseArrays
 
 julia> A = sprandn(100, 100, 0.1);
 
-julia> using Expokit
+julia> using ExponentialUtilities
 
 julia> E = SparseMatrixExp(A);
 
@@ -64,10 +56,12 @@ julia> get_columns(E, [10]); # same as get_column(E, 10) but a 100x1 matrix is r
 ### Notes
 
 This type is provided for use with very large and very sparse matrices.
-The evaluation of the exponential matrix action over vectors relies on the
-[Expokit](https://github.com/acroy/Expokit.jl) package. Hence, you will have to
-install and load this optional dependency to have access to the functionality
-of `SparseMatrixExp`.
+The evaluation of the exponential matrix action over vectors relies on external
+packages such as
+[ExponentialUtilities](https://github.com/SciML/ExponentialUtilities.jl) or
+[Expokit](https://github.com/acroy/Expokit.jl).
+Hence, you will have to install and load such an optional dependency to have
+access to the functionality of `SparseMatrixExp`.
 """
 struct SparseMatrixExp{N, MN<:AbstractSparseMatrix{N}} <: AbstractMatrix{N}
     M::MN
@@ -92,38 +86,39 @@ function size(spmexp::SparseMatrixExp, ax::Int)
     return size(spmexp.M, ax)
 end
 
-function get_column(spmexp::SparseMatrixExp{N}, j::Int) where {N}
-    require(:Expokit; fun_name="get_column")
-
+function get_column(spmexp::SparseMatrixExp{N}, j::Int;
+                    backend=get_exponential_backend()) where {N}
     n = size(spmexp, 1)
     aux = zeros(N, n)
     aux[j] = one(N)
-    return expmv(one(N), spmexp.M, aux)
+    return _expmv(backend, one(N), spmexp.M, aux)
 end
 
-function get_columns(spmexp::SparseMatrixExp{N}, J::AbstractArray) where {N}
-    require(:Expokit; fun_name="get_columns")
-
+function get_columns(spmexp::SparseMatrixExp{N}, J::AbstractArray;
+                     backend=get_exponential_backend()) where {N}
     n = size(spmexp, 1)
     aux = zeros(N, n)
     res = zeros(N, n, length(J))
     @inbounds for (k, j) in enumerate(J)
         aux[j] = one(N)
-        res[:, k] = expmv(one(N), spmexp.M, aux)
+        res[:, k] = _expmv(backend, one(N), spmexp.M, aux)
         aux[j] = zero(N)
     end
     return res
 end
 
 """
-    get_row(spmexp::SparseMatrixExp{N}, i::Int) where {N}
+    get_row(spmexp::SparseMatrixExp{N}, i::Int;
+            [backend]=get_exponential_backend()) where {N}
 
 Return a single row of a sparse matrix exponential.
 
 ### Input
 
-- `spmexp` -- sparse matrix exponential
-- `i`      -- row index
+- `spmexp`  -- sparse matrix exponential
+- `i`       -- row index
+- `backend` -- (optional; default: `get_exponential_backend()`) exponentiation
+               backend
 
 ### Output
 
@@ -135,25 +130,23 @@ This function uses Julia's `transpose` function to create the result.
 The result is of type `Transpose`; in Julia versions older than v0.7, the result
 was of type `RowVector`.
 """
-function get_row(spmexp::SparseMatrixExp{N}, i::Int) where {N}
-    require(:Expokit; fun_name="get_row")
-
+function get_row(spmexp::SparseMatrixExp{N}, i::Int;
+                 backend=get_exponential_backend()) where {N}
     n = size(spmexp, 1)
     aux = zeros(N, n)
     aux[i] = one(N)
-    return transpose(expmv(one(N), transpose(spmexp.M), aux))
+    return transpose(_expmv(backend, one(N), transpose(spmexp.M), aux))
 end
 
-function get_rows(spmexp::SparseMatrixExp{N}, I::AbstractArray{Int}) where {N}
-    require(:Expokit; fun_name="get_rows")
-
+function get_rows(spmexp::SparseMatrixExp{N}, I::AbstractArray{Int};
+                  backend=get_exponential_backend()) where {N}
     n = size(spmexp, 1)
     aux = zeros(N, n)
     res = zeros(N, length(I), n)
     Mtranspose = transpose(spmexp.M)
     @inbounds for (k, i) in enumerate(I)
         aux[i] = one(N)
-        res[k, :] = expmv(one(N), Mtranspose, aux)
+        res[k, :] = _expmv(backend, one(N), Mtranspose, aux)
         aux[i] = zero(N)
     end
     return res
@@ -287,14 +280,16 @@ function dim(em::ExponentialMap)
 end
 
 """
-    σ(d::AbstractVector, em::ExponentialMap)
+    σ(d::AbstractVector, em::ExponentialMap; [backend]=get_exponential_backend())
 
 Return the support vector of the exponential map.
 
 ### Input
 
-- `d`  -- direction
-- `em` -- exponential map
+- `d`       -- direction
+- `em`      -- exponential map
+- `backend` -- (optional; default: `get_exponential_backend()`) exponentiation
+               backend
 
 ### Output
 
@@ -305,28 +300,25 @@ If the direction has norm zero, the result depends on the wrapped set.
 
 If ``E = \\exp(M)⋅S``, where ``M`` is a matrix and ``S`` is a set, it
 follows that ``σ(d, E) = \\exp(M)⋅σ(\\exp(M)^T d, S)`` for any direction ``d``.
-
-We allow sparse direction vectors, but will convert them to dense vectors to be
-able to use `expmv`.
 """
-function σ(d::AbstractVector, em::ExponentialMap)
-    require(:Expokit; fun_name="σ")
-
+function σ(d::AbstractVector, em::ExponentialMap;
+           backend=get_exponential_backend())
     N = promote_type(eltype(d), eltype(em))
-    d_dense = d isa Vector ? d : Vector(d)
-    v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v   <- exp(M^T) * d
-    return expmv(one(N), em.spmexp.M, σ(v, em.X)) # res <- exp(M) * σ(v, S)
+    v = _expmv(backend, one(N), transpose(em.spmexp.M), d)  # v   <- exp(M^T) * d
+    return _expmv(backend, one(N), em.spmexp.M, σ(v, em.X)) # res <- exp(M) * σ(v, S)
 end
 
 """
-    ρ(d::AbstractVector, em::ExponentialMap)
+    ρ(d::AbstractVector, em::ExponentialMap; [backend]=get_exponential_backend())
 
 Return the support function of the exponential map.
 
 ### Input
 
-- `d`  -- direction
-- `em` -- exponential map
+- `d`       -- direction
+- `em`      -- exponential map
+- `backend` -- (optional; default: `get_exponential_backend()`) exponentiation
+               backend
 
 ### Output
 
@@ -336,16 +328,11 @@ The support function in the given direction.
 
 If ``E = \\exp(M)⋅S``, where ``M`` is a matrix and ``S`` is a set, it
 follows that ``ρ(d, E) = ρ(\\exp(M)^T d, S)`` for any direction ``d``.
-
-We allow sparse direction vectors, but will convert them to dense vectors to be
-able to use `expmv`.
 """
-function ρ(d::AbstractVector, em::ExponentialMap)
-    require(:Expokit; fun_name="ρ")
-
+function ρ(d::AbstractVector, em::ExponentialMap;
+           backend=get_exponential_backend())
     N = promote_type(eltype(d), eltype(em))
-    d_dense = d isa Vector ? d : Vector(d)
-    v = expmv(one(N), transpose(em.spmexp.M), d_dense) # v <- exp(M^T) * d
+    v = _expmv(backend, one(N), transpose(em.spmexp.M), d) # v <- exp(M^T) * d
     return ρ(v, em.X)
 end
 
@@ -354,14 +341,16 @@ function concretize(em::ExponentialMap)
 end
 
 """
-    ∈(x::AbstractVector, em::ExponentialMap)
+    ∈(x::AbstractVector, em::ExponentialMap; [backend]=get_exponential_backend())
 
 Check whether a given point is contained in an exponential map of a set.
 
 ### Input
 
-- `x`  -- point/vector
-- `em` -- exponential map of a set
+- `x`       -- point/vector
+- `em`      -- exponential map of a set
+- `backend` -- (optional; default: `get_exponential_backend()`) exponentiation
+               backend
 
 ### Output
 
@@ -387,24 +376,26 @@ julia> [1.0, 1.0] ∈ em
 true
 ```
 """
-function ∈(x::AbstractVector, em::ExponentialMap)
-    require(:Expokit; fun_name="∈")
-
+function ∈(x::AbstractVector, em::ExponentialMap;
+           backend=get_exponential_backend())
     @assert length(x) == dim(em) "a vector of length $(length(x)) is " *
         "incompatible with a set of dimension $(dim(em))"
     N = promote_type(eltype(x), eltype(em))
-    y = expmv(-one(N), em.spmexp.M, x)
+    y = _expmv(backend, -one(N), em.spmexp.M, x)
     return y ∈ em.X
 end
 
 """
-    vertices_list(em::ExponentialMap{N}) where {N}
+    vertices_list(em::ExponentialMap{N};
+                  [backend]=get_exponential_backend()) where {N}
 
 Return the list of vertices of a (polytopic) exponential map.
 
 ### Input
 
-- `em` -- exponential map
+- `em`      -- exponential map
+- `backend` -- (optional; default: `get_exponential_backend()`) exponentiation
+               backend
 
 ### Output
 
@@ -415,16 +406,15 @@ A list of vertices.
 We assume that the underlying set `X` is polytopic.
 Then the result is just the exponential map applied to the vertices of `X`.
 """
-function vertices_list(em::ExponentialMap{N}) where {N}
-    require(:Expokit; fun_name="vertices_list")
-
+function vertices_list(em::ExponentialMap{N};
+                       backend=get_exponential_backend()) where {N}
     # collect low-dimensional vertices lists
     vlist_X = vertices_list(em.X)
 
     # create resulting vertices list
     vlist = Vector{Vector{N}}(undef, length(vlist_X))
     @inbounds for (i, v) in enumerate(vlist_X)
-        vlist[i] = expmv(one(N), em.spmexp.M, v)
+        vlist[i] = _expmv(backend, one(N), em.spmexp.M, v)
     end
 
     return vlist
@@ -558,14 +548,17 @@ function dim(eprojmap::ExponentialProjectionMap)
 end
 
 """
-    σ(d::AbstractVector, eprojmap::ExponentialProjectionMap)
+    σ(d::AbstractVector, eprojmap::ExponentialProjectionMap;
+      [backend]=get_exponential_backend())
 
 Return the support vector of a projection of an exponential map.
 
 ### Input
 
 - `d`        -- direction
 - `eprojmap` -- projection of an exponential map
+- `backend`  -- (optional; default: `get_exponential_backend()`) exponentiation
+                backend
 
 ### Output
 
@@ -577,22 +570,17 @@ If the direction has norm zero, the result depends on the wrapped set.
 If ``S = (L⋅M⋅R)⋅X``, where ``L`` and ``R`` are matrices, ``M`` is a matrix
 exponential, and ``X`` is a set, it follows that
 ``σ(d, S) = L⋅M⋅R⋅σ(R^T⋅M^T⋅L^T⋅d, X)`` for any direction ``d``.
-
-We allow sparse direction vectors, but will convert them to dense vectors to be
-able to use `expmv`.
 """
-function σ(d::AbstractVector, eprojmap::ExponentialProjectionMap)
-    require(:Expokit; fun_name="σ")
-
-    d_dense = d isa Vector ? d : Vector(d)
-    daux = transpose(eprojmap.projspmexp.L) * d_dense
+function σ(d::AbstractVector, eprojmap::ExponentialProjectionMap;
+           backend=get_exponential_backend())
+    daux = transpose(eprojmap.projspmexp.L) * d
     N = promote_type(eltype(d), eltype(eprojmap))
-    aux1 = expmv(one(N), transpose(eprojmap.projspmexp.spmexp.M), daux)
+    aux1 = _expmv(backend, one(N), transpose(eprojmap.projspmexp.spmexp.M), daux)
     daux = _At_mul_B(eprojmap.projspmexp.R, aux1)
     svec = σ(daux, eprojmap.X)
 
     aux2 = eprojmap.projspmexp.R * svec
-    daux = expmv(one(N), eprojmap.projspmexp.spmexp.M, aux2)
+    daux = _expmv(backend, one(N), eprojmap.projspmexp.spmexp.M, aux2)
     return eprojmap.projspmexp.L * daux
 end
 

src/LazySets.jl

@@ -51,6 +51,7 @@ include("Utils/require.jl")
 include("Utils/helper_functions.jl")
 include("Utils/macros.jl")
 include("Utils/iterators.jl")
+include("Utils/matrix_exponential.jl")
 
 # ==================
 # Abstract set types