Use FiniteHorizonGramians.jl

Project.toml

@@ -10,6 +10,8 @@ DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+FiniteHorizonGramians = "b59a298d-d283-4a37-9369-85a9f9a111a5"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 GaussianDistributions = "43dcc890-d446-5863-8d1a-14597580bb8d"
@@ -45,6 +47,8 @@ DiffEqDevTools = "2"
 ExponentialUtilities = "1"
 FastBroadcast = "0.2"
 FastGaussQuadrature = "0.5, 1"
+FillArrays = "1.9"
+FiniteHorizonGramians = "0.1"
 ForwardDiff = "0.10"
 FunctionWrappersWrappers = "0.1.3"
 GaussianDistributions = "0.5"

benchmarks/hodgkinhuxley.jmd

@@ -34,7 +34,8 @@ Plots.theme(
 αh(V, VT) = 0.128 * exp(-(V - VT - 17) / 18)
 βh(V, VT) = 4 / (1 + exp(-(V - VT - 40) / 5))
 
-Inj(t) = (5 <= t <= 40) ? 500one(t) : zero(t)
+const current_tspan = (5, 40)
+Inj(t) = (current_tspan[1] <= t <= current_tspan[2]) ? 500one(t) : zero(t)
 
 function f(du, u, p, t)
     @unpack gNa, gK, ENa, EK, area, C, Eleak, VT, gleak = p
@@ -88,14 +89,13 @@ _setups = [
   "EK0(3)" => Dict(:alg=>EK0(order=3, smooth=DENSE))
   "EK1(2)" => Dict(:alg=>EK1(order=2, smooth=DENSE))
   "EK1(3)" => Dict(:alg=>EK1(order=3, smooth=DENSE))
-  "EK1(5)" => Dict(:alg=>EK1(order=5, smooth=DENSE))
+  "RosenbrockExpEK1(2)" => Dict(:alg=>RosenbrockExpEK(order=2, smooth=DENSE))
   "RosenbrockExpEK1(3)" => Dict(:alg=>RosenbrockExpEK(order=3, smooth=DENSE))
-  "RosenbrockExpEK1(5)" => Dict(:alg=>RosenbrockExpEK(order=5, smooth=DENSE))
 ]
 
 labels = first.(_setups)
 setups = last.(_setups)
-colors = [1 1 2 2 2 3 3]
+colors = [1 1 2 2 3 3]
 
 abstols = 1.0 ./ 10.0 .^ (6:10)
 reltols = 1.0 ./ 10.0 .^ (3:7)
@@ -108,6 +108,7 @@ wp = WorkPrecisionSet(
     save_everystep = SAVE_EVERYSTEP,
     maxiters = Int(1e7),
     numruns = 5,
+    tstops = current_tspan,
 )
 
 plot(wp, title="Adaptive steps - no smoothing", color=colors)
@@ -127,6 +128,7 @@ ref_wp_final = WorkPrecisionSet(
     dense = false,
     save_everystep = false,
     maxiters = Int(1e7),
+    tstops = current_tspan,
 )
 ref_wp_dense = WorkPrecisionSet(
     prob, abstols ./ 1000, reltols ./ 1000, ref_setups;
@@ -135,6 +137,7 @@ ref_wp_dense = WorkPrecisionSet(
     dense = true,
     save_everystep = true,
     maxiters = Int(1e7),
+    tstops = current_tspan,
 )
 
 plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
@@ -150,14 +153,13 @@ _setups = [
   "EK0(3)" => Dict(:alg=>EK0(order=3, smooth=DENSE))
   "EK1(2)" => Dict(:alg=>EK1(order=2, smooth=DENSE))
   "EK1(3)" => Dict(:alg=>EK1(order=3, smooth=DENSE))
-  "EK1(5)" => Dict(:alg=>EK1(order=5, smooth=DENSE))
+  "RosenbrockExpEK1(2)" => Dict(:alg=>RosenbrockExpEK(order=2, smooth=DENSE))
   "RosenbrockExpEK1(3)" => Dict(:alg=>RosenbrockExpEK(order=3, smooth=DENSE))
-  "RosenbrockExpEK1(5)" => Dict(:alg=>RosenbrockExpEK(order=5, smooth=DENSE))
 ]
 
 labels = first.(_setups)
 setups = last.(_setups)
-colors = [1 1 2 2 2 3 3]
+colors = [1 1 2 2 3 3]
 
 abstols = 1.0 ./ 10.0 .^ (6:10)
 reltols = 1.0 ./ 10.0 .^ (3:7)
@@ -170,6 +172,7 @@ wp = WorkPrecisionSet(
     save_everystep = SAVE_EVERYSTEP,
     maxiters = Int(1e7),
     numruns = 5,
+    tstops = current_tspan,
 )
 
 plot(wp, title="Adaptive steps - with smoothing", color=colors)