Skip to content

Commit

Permalink
fix Pade approximation and test _exp
Browse files Browse the repository at this point in the history
  • Loading branch information
schillic committed Apr 22, 2023
1 parent 0896cbb commit 9f58ef9
Show file tree
Hide file tree
Showing 6 changed files with 48 additions and 12 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
28 changes: 17 additions & 11 deletions src/Discretization/exponentiation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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)
Expand Down
3 changes: 3 additions & 0 deletions src/Initialization/init.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down
1 change: 1 addition & 0 deletions src/Initialization/init_Expokit.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
eval(load_expokit_pade())
24 changes: 24 additions & 0 deletions test/discretization.jl
Original file line number Diff line number Diff line change
@@ -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
2 changes: 2 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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")
Expand Down

0 comments on commit 9f58ef9

Please sign in to comment.