From 9f58ef97cabe785eca1ea5190f120404b2d64133 Mon Sep 17 00:00:00 2001
From: schillic <git@christianschilling.net>
Date: Sat, 22 Apr 2023 14:58:21 +0200
Subject: [PATCH] fix Pade approximation and test _exp

---
 Project.toml                         |  2 +-
 src/Discretization/exponentiation.jl | 28 +++++++++++++++++-----------
 src/Initialization/init.jl           |  3 +++
 src/Initialization/init_Expokit.jl   |  1 +
 test/discretization.jl               | 24 ++++++++++++++++++++++++
 test/runtests.jl                     |  2 ++
 6 files changed, 48 insertions(+), 12 deletions(-)
 create mode 100644 src/Initialization/init_Expokit.jl
 create mode 100644 test/discretization.jl

diff --git a/Project.toml b/Project.toml
index c80d4abae0..3cfed2ea61 100644
--- a/Project.toml
+++ b/Project.toml
@@ -76,4 +76,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["CDDLib", "DifferentialEquations", "FastExpm", "Flowstar", "JuMP", "Plots", "Polyhedra", "SparseArrays", "Symbolics", "Test"]
+test = ["CDDLib", "DifferentialEquations", "Expokit", "FastExpm", "Flowstar", "JuMP", "Plots", "Polyhedra", "SparseArrays", "Symbolics", "Test"]
diff --git a/src/Discretization/exponentiation.jl b/src/Discretization/exponentiation.jl
index e90ed3e7c4..71bb571d2a 100644
--- a/src/Discretization/exponentiation.jl
+++ b/src/Discretization/exponentiation.jl
@@ -113,7 +113,19 @@ _alias(alg::Val{:pade}) = PadeExp
 @inline _exp(A::AbstractMatrix, ::LazyExpAlg) = SparseMatrixExp(A)
 
 # pade approximants (requires Expokit.jl)
-@inline _exp(A::SparseMatrixCSC, ::PadeExpAlg) = padm(A)
+@inline function _exp(::AbstractMatrix, ::PadeExpAlg)
+    throw(ArgumentError("algorithm requires a sparse matrix"))
+end
+
+@inline function _exp(::SparseMatrixCSC, ::PadeExpAlg)
+    @requires Expokit
+end
+
+function load_expokit_pade()
+    return quote
+        @inline _exp(A::SparseMatrixCSC, ::PadeExpAlg) = Expokit.padm(A)
+    end
+end  # quote / load_expokit_pade
 
 function _exp(A::AbstractMatrix, alg::IntervalExpAlg)
     return exp_overapproximation(_interval_matrix(A), one(eltype(A)), alg.order)
@@ -165,21 +177,15 @@ A matrix or a lazy wrapper of the matrix exponential, depending on `alg`.
 
 ### Notes
 
-If the algorithm `LazyExp` is used, evaluations of the action of the matrix
-exponential are done with the `expmv` implementation from `Expokit`
-(but see `LazySets#1312` for the planned generalization to other backends).
+If the algorithm `LazyExp` is used, actions of the matrix exponential are
+evaluated with an external library such as `ExponentialUtilities.jl` or
+`Expokit.jl`.
 """
 function _exp(A::AbstractMatrix, δ, alg::AbstractExpAlg=BaseExp)
     n = checksquare(A)
     return _exp(A * δ, alg)
 end
 
-function _exp(A::AbstractMatrix, δ, alg::PadeExpAlg)
-    n = checksquare(A)
-    @requires Expokit
-    return _exp(A * δ, alg)
-end
-
 @inline function _P_2n(A::AbstractMatrix{N}, δ, n) where {N}
     return [Matrix(A*δ) Matrix(δ*I, n, n);
             zeros(n, 2*n)]::Matrix{N}
@@ -431,7 +437,7 @@ function _Φ₂_inv(A::IdentityMultiple, δ, alg, Φ=nothing)
     return IdentityMultiple(α, size(A, 1))
 end
 
-function _Eplus(A::SparseMatrixCSC{N, D}, X0::AbstractHyperrectangle{N}, δt; m=min(30, size(A, 1)), tol=1e-7) where {N, D}    
+function _Eplus(A::SparseMatrixCSC{N, D}, X0::AbstractHyperrectangle{N}, δt; m=min(30, size(A, 1)), tol=1e-7) where {N, D}
     n = dim(X0)
     A2 = A * A; # fast if A sparse
     V = symmetric_interval_hull(A2 * X0)
diff --git a/src/Initialization/init.jl b/src/Initialization/init.jl
index d744f001c4..5c2bd5d5f9 100644
--- a/src/Initialization/init.jl
+++ b/src/Initialization/init.jl
@@ -121,6 +121,9 @@ function __init__()
     # sparse dynamic representation of multivariate polynomials
     @require DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" include("init_DynamicPolynomials.jl")
 
+    # exponentiation methods using Krylov subspace and Padé approximations
+    @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("init_Expokit.jl")
+
     # exponentiation methods using Krylov subspace approximations
     @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" include("init_ExponentialUtilities.jl")
 
diff --git a/src/Initialization/init_Expokit.jl b/src/Initialization/init_Expokit.jl
new file mode 100644
index 0000000000..5129b9d075
--- /dev/null
+++ b/src/Initialization/init_Expokit.jl
@@ -0,0 +1 @@
+eval(load_expokit_pade())
diff --git a/test/discretization.jl b/test/discretization.jl
new file mode 100644
index 0000000000..012d22070f
--- /dev/null
+++ b/test/discretization.jl
@@ -0,0 +1,24 @@
+using SparseArrays
+using ReachabilityAnalysis: _exp
+
+@testset "Exponentiation" begin
+    A = [10.0 0; 0 20]
+    δ = 0.1
+    expAδ = [exp(1) 0; 0 exp(2)]
+
+    # default algorithm
+    @test _exp(A, δ) == expAδ
+
+    # base algorithm
+    @test _exp(A, δ, RA.BaseExpAlg()) == expAδ
+
+    # Krylov algorithm
+    @test _exp(sparse(A), δ, RA.LazyExpAlg()) == expAδ
+
+    # interval matrix
+    @test all((_exp(A, δ, RA.IntervalExpAlg()) .- IntervalMatrix(expAδ)) .<= 1e-4)
+
+    # Padé approximation
+    @test_throws ArgumentError _exp(A, δ, RA.PadeExpAlg())
+    @test _exp(sparse(A), δ, RA.PadeExpAlg()) ≈ expAδ
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 917b4a087d..8caf817f35 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,6 +13,7 @@ using StaticArrays
 using Polyhedra, CDDLib # for VREP algorithm
 import TaylorModels # for TMJets algorithm
 import Flowstar # for FLOWSTAR algorithm
+import Expokit # for Padé approximation
 
 # fix namespace conflicts with Polyhedra
 using LazySets: dim, HalfSpace, Interval, Line2D, translate, project
@@ -35,6 +36,7 @@ include("models/generated/Brusselator.jl")
 include("models/hybrid/thermostat.jl")
 
 include("solve.jl")
+include("discretization.jl")
 include("reachsets.jl")
 include("flowpipes.jl")
 include("traces.jl")