Merge pull request #770 from JuliaReach/schillic/init

Fix optional packages in Exponentiation module

src/Discretization/Exponentiation.jl

@@ -136,11 +136,15 @@ _exp(A::AbstractMatrix, ::LazyExpAlg) = SparseMatrixExp(A)
 
 # pade approximants (requires Expokit.jl)
 _exp(::AbstractMatrix, ::PadeExpAlg) = throw(ArgumentError("algorithm requires a sparse matrix"))
-_exp(::SparseMatrixCSC, ::PadeExpAlg) = require(@__MODULE__, Expokit)
+
+function _exp(A::SparseMatrixCSC, alg::PadeExpAlg)
+    require(@__MODULE__, :Expokit)
+    return _exp_pade(A, alg)
+end
 
 function load_expokit_pade()
     return quote
-        @inline Exponentiation._exp(A::SparseMatrixCSC, ::PadeExpAlg) = Expokit.padm(A)
+        _exp_pade(A::SparseMatrixCSC, ::PadeExpAlg) = Expokit.padm(A)
     end
 end  # quote / load_expokit_pade
 
@@ -457,7 +461,7 @@ function _Eplus(A::SparseMatrixCSC{N,D}, X0::AbstractHyperrectangle{N}, δt; m=m
     v = V.radius
     Aabs = copy(abs.(A))
 
-    require(@__MODULE__, ExponentialUtilities)
+    require(@__MODULE__, :ExponentialUtilities)
     Pv = _phiv(Aabs, v, 1, δt; m, tol)
     return E⁺ = Hyperrectangle(zeros(n), Pv)
 end
@@ -471,4 +475,15 @@ function elementwise_abs(A::SparseMatrixCSC)
 end
 elementwise_abs(A::IdentityMultiple) = IdentityMultiple(abs(A.M.λ), size(A, 1))
 
+# =====================
+# Optional dependencies
+# =====================
+
+using Requires
+
+function __init__()
+    @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" eval(load_expokit_pade())
+    @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" include("init_ExponentialUtilities.jl")
+end
+
 end # module

src/Discretization/init_ExponentialUtilities.jl

@@ -0,0 +1,29 @@
+using .ExponentialUtilities: phiv!, arnoldi!, KrylovSubspace
+
+function _phiv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)
+    # initialization of the krylov subspace
+    TA, Tb = eltype(A), eltype(b)
+    T = promote_type(TA, Tb)
+    Ks = KrylovSubspace{T,real(T)}(length(b), m)
+    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)
+
+    out = Matrix{Float64}(undef, size(A, 1), 3)
+    phiv!(out, NSTEPS * dt, Ks, 2)
+
+    return view(out, :, 3) .* dt^2
+end
+
+# # Compute out <- exp(A * NSTEPS * dt) * b  (currently not used)
+# using .ExponentialUtilities: expv!
+# function _expv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)
+#     # initialization of the krylov subspace
+#     TA, Tb = eltype(A), eltype(b)
+#     T = promote_type(TA, Tb)
+#     Ks = KrylovSubspace{T,real(T)}(length(b), m)
+#     arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)
+#
+#     out = similar(b)
+#     expv!(out, NSTEPS * dt, Ks)
+#
+#     return out
+# end

src/Initialization/init.jl

@@ -115,15 +115,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")
 
-    # exponentiation methods for large sparse matrices
-    @require FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2" include("init_FastExpm.jl")
-
     # external reachability backend: Flow* C++ library
     @require Flowstar = "a8054ddd-9dca-4d20-8ffe-ae96ec1541f1" include("init_Flowstar.jl")
 

src/Initialization/init_ExponentialUtilities.jl

@@ -1,31 +1,4 @@
-using .ExponentialUtilities: expv!, phiv!, arnoldi!, KrylovSubspace
-
-# Compute out <- exp(A * NSTEPS * dt) * b
-function _expv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)
-    # initialization of the krylov subspace
-    TA, Tb = eltype(A), eltype(b)
-    T = promote_type(TA, Tb)
-    Ks = KrylovSubspace{T,real(T)}(length(b), m)
-    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)
-
-    out = similar(b)
-    expv!(out, NSTEPS * dt, Ks)
-
-    return out
-end
-
-function _phiv(A, b, NSTEPS, dt; hermitian=false, m=min(30, size(A, 1)), tol=1e-7)
-    # initialization of the krylov subspace
-    TA, Tb = eltype(A), eltype(b)
-    T = promote_type(TA, Tb)
-    Ks = KrylovSubspace{T,real(T)}(length(b), m)
-    arnoldi!(Ks, A, b; m=m, ishermitian=hermitian, tol=tol)
-
-    out = Matrix{Float64}(undef, size(A, 1), 3)
-    phiv!(out, NSTEPS * dt, Ks, 2)
-
-    return view(out, :, 3) .* dt^2
-end
+using .ExponentialUtilities: expv!, arnoldi!, KrylovSubspace
 
 eval(load_krylov_LGG09_homog())
 eval(load_krylov_LGG09_inhomog())

test/Project.toml

@@ -4,7 +4,6 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
-FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
 Flowstar = "a8054ddd-9dca-4d20-8ffe-ae96ec1541f1"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -26,7 +25,6 @@ DifferentialEquations = "6.17, 7"
 DynamicPolynomials = "0.3 - 0.5"
 Expokit = "0.2"
 ExponentialUtilities = "1"
-FastExpm = "1"
 Flowstar = "0.2.4"
 IntervalArithmetic = "0.16 - 0.20, =0.20.9"  # new versions require updates and are incompatible with dependencies
 JuMP = "0.21 - 0.23, 1"