diff --git a/benchmarks/lotkavolterra.jmd b/benchmarks/lotkavolterra.jmd
index 0d63a718f..d5719ed12 100644
--- a/benchmarks/lotkavolterra.jmd
+++ b/benchmarks/lotkavolterra.jmd
@@ -175,44 +175,50 @@ wp = WorkPrecisionSet(
plot(wp, color=[2 2 2 3 3 3], xticks = 10.0 .^ (-16:1:5))
```
-## TaylorModeInit vs ClassicSolverInit
+## Comparison of the different initialization schemes
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
-_setups = [
- "EK1(2) TaylorInit" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(3) TaylorInit" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(5) TaylorInit" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(2) Tsit5Init" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(3) Tsit5Init" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(5) Tsit5Init" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(2) Tsit5Init+ddu" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
- "EK1(3) Tsit5Init+ddu" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
- "EK1(5) Tsit5Init+ddu" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
-]
-
-labels = first.(_setups)
-setups = last.(_setups)
-
abstols = 1.0 ./ 10.0 .^ (4:14)
reltols = 1.0 ./ 10.0 .^ (1:11)
-wp = WorkPrecisionSet(
- prob, abstols, reltols, setups;
- names = labels,
- #print_names = true,
- appxsol = test_sol,
- dense = DENSE,
- save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
- maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
+orders = (2, 3, 5, 8)
+ps = []
+for o in orders
+ _setups = [
+ "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+ "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+ "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+ "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit()))
+ ]
+
+ labels = first.(_setups)
+ setups = last.(_setups)
+
+ wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+ )
+
+ p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+ push!(ps, p)
+end
+plot(
+ ps...,
+ layout=(length(orders), 1),
+ size = (1000, length(orders)*300),
+ xlabel=["" "" "" "Error"],
)
-
-plot(wp, color=[2 2 2 4 4 4 5 5 5], xticks = 10.0 .^ (-16:1:5))
```
diff --git a/benchmarks/vanderpol.jmd b/benchmarks/vanderpol.jmd
index 72c836e56..f6992ede9 100644
--- a/benchmarks/vanderpol.jmd
+++ b/benchmarks/vanderpol.jmd
@@ -64,6 +64,90 @@ wp = WorkPrecisionSet(
plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
```
+## Comparison of the different initialization schemes
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+orders = (3, 5, 8)
+ps = []
+for o in orders
+ _setups = [
+ "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+ "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+ "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+ # "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit())) # unstable
+ ]
+
+ labels = first.(_setups)
+ setups = last.(_setups)
+
+ wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+ )
+
+ p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+ push!(ps, p)
+end
+plot(
+ ps...,
+ layout=(length(orders), 1),
+ size = (1000, length(orders)*300),
+ xlabel=["" "" "" "Error"],
+)
+```
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (2:5)
+
+_setups = [
+ "EK1(8) TaylorInit(8)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(8)))
+ "EK1(8) TaylorInit(7)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(7)))
+ "EK1(8) TaylorInit(6)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(6)))
+ "EK1(8) TaylorInit(5)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(5)))
+ "EK1(8) TaylorInit(4)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(4)))
+ "EK1(8) TaylorInit(3)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(3)))
+ "EK1(8) TaylorInit(2)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(2)))
+ # "EK1(8) TaylorInit(1)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(1))) # fails, see above
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+)
+
+plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
+```
+
+
## Solving the first- vs second-order ODE
```julia
diff --git a/docs/Project.toml b/docs/Project.toml
index 82cd2b515..4bf36e8fe 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,3 +12,4 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
ProbNumDiffEq = "bf3e78b0-7d74-48a5-b855-9609533b56a5"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/docs/src/benchmarks/figures/lotkavolterra_2_1.svg b/docs/src/benchmarks/figures/lotkavolterra_2_1.svg
index 1cd567d89..1c65f4707 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_2_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_2_1.svg
@@ -1,148 +1,148 @@
diff --git a/docs/src/benchmarks/figures/lotkavolterra_3_1.svg b/docs/src/benchmarks/figures/lotkavolterra_3_1.svg
index 6732a66e5..af5fc0f55 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_3_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_3_1.svg
@@ -1,318 +1,318 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/lotkavolterra_4_1.svg b/docs/src/benchmarks/figures/lotkavolterra_4_1.svg
index 8c1532e7f..a79cbf88d 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_4_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_4_1.svg
@@ -1,320 +1,320 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/lotkavolterra_5_1.svg b/docs/src/benchmarks/figures/lotkavolterra_5_1.svg
index 04c2927a6..6f1183f7c 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_5_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_5_1.svg
@@ -1,380 +1,380 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/lotkavolterra_6_1.svg b/docs/src/benchmarks/figures/lotkavolterra_6_1.svg
index 3ef8c6070..ba78c9332 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_6_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_6_1.svg
@@ -1,366 +1,364 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/lotkavolterra_7_1.svg b/docs/src/benchmarks/figures/lotkavolterra_7_1.svg
index 17783dfca..8bd7fe31e 100644
--- a/docs/src/benchmarks/figures/lotkavolterra_7_1.svg
+++ b/docs/src/benchmarks/figures/lotkavolterra_7_1.svg
@@ -1,394 +1,1144 @@
-
diff --git a/docs/src/benchmarks/figures/vanderpol_2_1.svg b/docs/src/benchmarks/figures/vanderpol_2_1.svg
index 14bc5c40c..7b0e72aba 100644
--- a/docs/src/benchmarks/figures/vanderpol_2_1.svg
+++ b/docs/src/benchmarks/figures/vanderpol_2_1.svg
@@ -1,140 +1,140 @@
diff --git a/docs/src/benchmarks/figures/vanderpol_3_1.svg b/docs/src/benchmarks/figures/vanderpol_3_1.svg
index 822317d0e..79360d4f3 100644
--- a/docs/src/benchmarks/figures/vanderpol_3_1.svg
+++ b/docs/src/benchmarks/figures/vanderpol_3_1.svg
@@ -1,227 +1,227 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/vanderpol_4_1.svg b/docs/src/benchmarks/figures/vanderpol_4_1.svg
index bc33aed36..8802fbb69 100644
--- a/docs/src/benchmarks/figures/vanderpol_4_1.svg
+++ b/docs/src/benchmarks/figures/vanderpol_4_1.svg
@@ -1,140 +1,501 @@
-
diff --git a/docs/src/benchmarks/figures/vanderpol_5_1.svg b/docs/src/benchmarks/figures/vanderpol_5_1.svg
index ee17c679b..032fd7beb 100644
--- a/docs/src/benchmarks/figures/vanderpol_5_1.svg
+++ b/docs/src/benchmarks/figures/vanderpol_5_1.svg
@@ -1,256 +1,165 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/benchmarks/figures/vanderpol_6_1.svg b/docs/src/benchmarks/figures/vanderpol_6_1.svg
new file mode 100644
index 000000000..4b36cb9bf
--- /dev/null
+++ b/docs/src/benchmarks/figures/vanderpol_6_1.svg
@@ -0,0 +1,140 @@
+
+
diff --git a/docs/src/benchmarks/figures/vanderpol_7_1.svg b/docs/src/benchmarks/figures/vanderpol_7_1.svg
new file mode 100644
index 000000000..3bcf19a8c
--- /dev/null
+++ b/docs/src/benchmarks/figures/vanderpol_7_1.svg
@@ -0,0 +1,256 @@
+
+
diff --git a/docs/src/benchmarks/lotkavolterra.md b/docs/src/benchmarks/lotkavolterra.md
index ae584b8f7..932edba74 100644
--- a/docs/src/benchmarks/lotkavolterra.md
+++ b/docs/src/benchmarks/lotkavolterra.md
@@ -198,44 +198,50 @@ plot(wp, color=[2 2 2 3 3 3], xticks = 10.0 .^ (-16:1:5))
-## TaylorModeInit vs ClassicSolverInit
+## Comparison of the different initialization schemes
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
-_setups = [
- "EK1(2) TaylorInit" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(3) TaylorInit" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(5) TaylorInit" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=TaylorModeInit()))
- "EK1(2) Tsit5Init" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(3) Tsit5Init" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(5) Tsit5Init" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit()))
- "EK1(2) Tsit5Init+ddu" => Dict(:alg => EK1(order=2, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
- "EK1(3) Tsit5Init+ddu" => Dict(:alg => EK1(order=3, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
- "EK1(5) Tsit5Init+ddu" => Dict(:alg => EK1(order=5, smooth=DENSE, initialization=ClassicSolverInit(init_on_ddu=true)))
-]
-
-labels = first.(_setups)
-setups = last.(_setups)
-
abstols = 1.0 ./ 10.0 .^ (4:14)
reltols = 1.0 ./ 10.0 .^ (1:11)
-wp = WorkPrecisionSet(
- prob, abstols, reltols, setups;
- names = labels,
- #print_names = true,
- appxsol = test_sol,
- dense = DENSE,
- save_everystep = SAVE_EVERYSTEP,
- numruns = 10,
- maxiters = Int(1e7),
- timeseries_errors = false,
- verbose = false,
+orders = (2, 3, 5, 8)
+ps = []
+for o in orders
+ _setups = [
+ "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+ "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+ "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+ "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit()))
+ ]
+
+ labels = first.(_setups)
+ setups = last.(_setups)
+
+ wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+ )
+
+ p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+ push!(ps, p)
+end
+plot(
+ ps...,
+ layout=(length(orders), 1),
+ size = (1000, length(orders)*300),
+ xlabel=["" "" "" "Error"],
)
-
-plot(wp, color=[2 2 2 4 4 4 5 5 5], xticks = 10.0 .^ (-16:1:5))
```
@@ -275,7 +281,6 @@ Platform Info:
Environment:
JULIA_NUM_THREADS = auto
JULIA_STACKTRACE_MINIMAL = true
- JULIA_IMAGE_THREADS = 1
```
@@ -302,7 +307,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Project.toml`
[65888b18] ParameterizedFunctions v5.16.0
[91a5bcdd] Plots v1.39.0
[bf3e78b0] ProbNumDiffEq v0.13.0 `~/.julia/dev/ProbNumDiffEq`
- [0bca4576] SciMLBase v2.6.0
+ [0bca4576] SciMLBase v2.7.3
[505e40e9] SciPyDiffEq v0.2.1
[90137ffa] StaticArrays v1.6.5
[c3572dad] Sundials v4.20.1
@@ -329,7 +334,6 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[ec485272] ArnoldiMethod v0.2.0
[c9d4266f] ArrayAllocators v0.3.0
[4fba245c] ArrayInterface v7.5.1
- [30b0a656] ArrayInterfaceCore v0.1.29
[6e4b80f9] BenchmarkTools v1.3.2
[e2ed5e7c] Bijections v0.1.6
[d1d4a3ce] BitFlags v0.1.7
@@ -366,14 +370,14 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[864edb3b] DataStructures v0.18.15
[e2d170a0] DataValueInterfaces v1.0.0
[8bb1440f] DelimitedFiles v1.9.1
- [2b5f629d] DiffEqBase v6.136.0
+ [2b5f629d] DiffEqBase v6.138.0
[459566f4] DiffEqCallbacks v2.33.1
[f3b72e0c] DiffEqDevTools v2.39.0
[77a26b50] DiffEqNoiseProcess v5.19.0
[163ba53b] DiffResults v1.1.0
[b552c78f] DiffRules v1.15.1
[b4f34e82] Distances v0.10.10
- [31c24e10] Distributions v0.25.102
+ [31c24e10] Distributions v0.25.103
[ffbed154] DocStringExtensions v0.9.3
⌅ [5b8099bc] DomainSets v0.6.7
[fa6b7ba4] DualNumbers v0.6.8
@@ -442,7 +446,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[7ed4a6bd] LinearSolve v2.15.0
[2ab3a3ac] LogExpFunctions v0.3.26
[e6f89c97] LoggingExtras v1.0.3
- [bdcacae8] LoopVectorization v0.12.165
+ [bdcacae8] LoopVectorization v0.12.166
[10e44e05] MATLAB v0.8.4
[e2752cbe] MATLABDiffEq v1.2.0
[33e6dc65] MKL v0.6.1
@@ -461,7 +465,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[2774e3e8] NLsolve v4.5.1
[77ba4419] NaNMath v1.0.2
⌅ [356022a1] NamedDims v0.2.50
- [8913a72c] NonlinearSolve v2.6.0
+ [8913a72c] NonlinearSolve v2.6.1
[54ca160b] ODEInterface v0.5.0
[09606e27] ODEInterfaceDiffEq v3.13.3
[6fd5a793] Octavian v0.3.27
@@ -514,11 +518,11 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[fdea26ae] SIMD v3.4.5
[94e857df] SIMDTypes v0.1.0
[476501e8] SLEEFPirates v0.6.42
- [0bca4576] SciMLBase v2.6.0
+ [0bca4576] SciMLBase v2.7.3
[e9a6253c] SciMLNLSolve v0.1.9
[c0aeaf25] SciMLOperators v0.3.6
[505e40e9] SciPyDiffEq v0.2.1
- [6c6a2e73] Scratch v1.2.0
+ [6c6a2e73] Scratch v1.2.1
[91c51154] SentinelArrays v1.4.0
[efcf1570] Setfield v1.1.1
[1277b4bf] ShiftedArrays v2.0.0
@@ -560,7 +564,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[a759f4b9] TimerOutputs v0.5.23
[c751599d] ToeplitzMatrices v0.8.2
[0796e94c] Tokenize v0.5.25
- [3bb67fe8] TranscodingStreams v0.10.1
+ [3bb67fe8] TranscodingStreams v0.10.2
[a2a6695c] TreeViews v0.3.0
[d5829a12] TriangularSolve v0.1.20
[410a4b4d] Tricks v0.1.8
@@ -579,7 +583,6 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[1b915085] WinReg v1.0.0
[ddb6d928] YAML v0.4.9
[c2297ded] ZMQ v1.2.2
- [700de1a5] ZygoteRules v0.2.4
[0518478a] deSolveDiffEq v0.1.1
[6e34b625] Bzip2_jll v1.0.8+0
[83423d85] Cairo_jll v1.16.1+1
diff --git a/docs/src/benchmarks/vanderpol.md b/docs/src/benchmarks/vanderpol.md
index 1c0a5df39..2fa57f023 100644
--- a/docs/src/benchmarks/vanderpol.md
+++ b/docs/src/benchmarks/vanderpol.md
@@ -75,6 +75,96 @@ plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^
+## Comparison of the different initialization schemes
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (6:11)
+reltols = 1.0 ./ 10.0 .^ (3:8)
+
+orders = (3, 5, 8)
+ps = []
+for o in orders
+ _setups = [
+ "EK1($o) TaylorInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=TaylorModeInit(o)))
+ "EK1($o) ForwardDiffInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ForwardDiffInit(o)))
+ "EK1($o) SimpleInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=SimpleInit()))
+ # "EK1($o) ClassicSolverInit" => Dict(:alg => EK1(order=o, smooth=DENSE, initialization=ClassicSolverInit())) # unstable
+ ]
+
+ labels = first.(_setups)
+ setups = last.(_setups)
+
+ wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+ )
+
+ p = plot(wp, color=[2 4 5 6], xticks = 10.0 .^ (-16:1:5))
+ push!(ps, p)
+end
+plot(
+ ps...,
+ layout=(length(orders), 1),
+ size = (1000, length(orders)*300),
+ xlabel=["" "" "" "Error"],
+)
+```
+
+
+
+```julia
+DENSE = false;
+SAVE_EVERYSTEP = false;
+
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (2:5)
+
+_setups = [
+ "EK1(8) TaylorInit(8)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(8)))
+ "EK1(8) TaylorInit(7)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(7)))
+ "EK1(8) TaylorInit(6)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(6)))
+ "EK1(8) TaylorInit(5)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(5)))
+ "EK1(8) TaylorInit(4)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(4)))
+ "EK1(8) TaylorInit(3)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(3)))
+ "EK1(8) TaylorInit(2)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(2)))
+ # "EK1(8) TaylorInit(1)" => Dict(:alg => EK1(order=8, smooth=DENSE, initialization=TaylorModeInit(1))) # fails, see above
+]
+
+labels = first.(_setups)
+setups = last.(_setups)
+
+wp = WorkPrecisionSet(
+ prob, abstols, reltols, setups;
+ names = labels,
+ #print_names = true,
+ appxsol = test_sol,
+ dense = DENSE,
+ save_everystep = SAVE_EVERYSTEP,
+ numruns = 10,
+ maxiters = Int(1e7),
+ timeseries_errors = false,
+ verbose = false,
+)
+
+plot(wp, palette=Plots.palette([:blue, :red], length(_setups)), xticks = 10.0 .^ (-16:1:5))
+```
+
+
+
+
+
+
## Solving the first- vs second-order ODE
```julia
@@ -91,7 +181,7 @@ test_sol2 = solve(prob2, RadauIIA5(), abstol=1/10^14, reltol=1/10^14, dense=fals
plot(test_sol2, title="Van der Pol Solution (2nd order)", legend=false, ylims=(-2.5, 2.5))
```
-
+
```julia
DENSE = false;
@@ -130,7 +220,7 @@ wp = WorkPrecisionSet(
plot(wp, color=[1 1 1 1 2 2 2 2], xticks = 10.0 .^ (-16:1:5))
```
-
+
@@ -164,7 +254,6 @@ Platform Info:
Environment:
JULIA_NUM_THREADS = auto
JULIA_STACKTRACE_MINIMAL = true
- JULIA_IMAGE_THREADS = 1
```
@@ -191,7 +280,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Project.toml`
[65888b18] ParameterizedFunctions v5.16.0
[91a5bcdd] Plots v1.39.0
[bf3e78b0] ProbNumDiffEq v0.13.0 `~/.julia/dev/ProbNumDiffEq`
- [0bca4576] SciMLBase v2.6.0
+ [0bca4576] SciMLBase v2.7.3
[505e40e9] SciPyDiffEq v0.2.1
[90137ffa] StaticArrays v1.6.5
[c3572dad] Sundials v4.20.1
@@ -218,7 +307,6 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[ec485272] ArnoldiMethod v0.2.0
[c9d4266f] ArrayAllocators v0.3.0
[4fba245c] ArrayInterface v7.5.1
- [30b0a656] ArrayInterfaceCore v0.1.29
[6e4b80f9] BenchmarkTools v1.3.2
[e2ed5e7c] Bijections v0.1.6
[d1d4a3ce] BitFlags v0.1.7
@@ -255,14 +343,14 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[864edb3b] DataStructures v0.18.15
[e2d170a0] DataValueInterfaces v1.0.0
[8bb1440f] DelimitedFiles v1.9.1
- [2b5f629d] DiffEqBase v6.136.0
+⌃ [2b5f629d] DiffEqBase v6.138.0
[459566f4] DiffEqCallbacks v2.33.1
[f3b72e0c] DiffEqDevTools v2.39.0
[77a26b50] DiffEqNoiseProcess v5.19.0
[163ba53b] DiffResults v1.1.0
[b552c78f] DiffRules v1.15.1
[b4f34e82] Distances v0.10.10
- [31c24e10] Distributions v0.25.102
+ [31c24e10] Distributions v0.25.103
[ffbed154] DocStringExtensions v0.9.3
⌅ [5b8099bc] DomainSets v0.6.7
[fa6b7ba4] DualNumbers v0.6.8
@@ -328,10 +416,10 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[50d2b5c4] Lazy v0.15.1
[1d6d02ad] LeftChildRightSiblingTrees v0.2.0
[d3d80556] LineSearches v7.2.0
- [7ed4a6bd] LinearSolve v2.15.0
+⌃ [7ed4a6bd] LinearSolve v2.15.0
[2ab3a3ac] LogExpFunctions v0.3.26
[e6f89c97] LoggingExtras v1.0.3
- [bdcacae8] LoopVectorization v0.12.165
+ [bdcacae8] LoopVectorization v0.12.166
[10e44e05] MATLAB v0.8.4
[e2752cbe] MATLABDiffEq v1.2.0
[33e6dc65] MKL v0.6.1
@@ -350,7 +438,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[2774e3e8] NLsolve v4.5.1
[77ba4419] NaNMath v1.0.2
⌅ [356022a1] NamedDims v0.2.50
- [8913a72c] NonlinearSolve v2.6.0
+⌃ [8913a72c] NonlinearSolve v2.6.1
[54ca160b] ODEInterface v0.5.0
[09606e27] ODEInterfaceDiffEq v3.13.3
[6fd5a793] Octavian v0.3.27
@@ -371,7 +459,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[995b91a9] PlotUtils v1.3.5
[91a5bcdd] Plots v1.39.0
[e409e4f3] PoissonRandom v0.4.4
- [f517fe37] Polyester v0.7.8
+⌃ [f517fe37] Polyester v0.7.8
[1d0040c9] PolyesterWeave v0.2.1
⌅ [f27b6e38] Polynomials v3.2.13
[2dfb63ee] PooledArrays v1.4.3
@@ -403,23 +491,23 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[fdea26ae] SIMD v3.4.5
[94e857df] SIMDTypes v0.1.0
[476501e8] SLEEFPirates v0.6.42
- [0bca4576] SciMLBase v2.6.0
+ [0bca4576] SciMLBase v2.7.3
[e9a6253c] SciMLNLSolve v0.1.9
- [c0aeaf25] SciMLOperators v0.3.6
+⌃ [c0aeaf25] SciMLOperators v0.3.6
[505e40e9] SciPyDiffEq v0.2.1
- [6c6a2e73] Scratch v1.2.0
- [91c51154] SentinelArrays v1.4.0
+ [6c6a2e73] Scratch v1.2.1
+⌃ [91c51154] SentinelArrays v1.4.0
[efcf1570] Setfield v1.1.1
[1277b4bf] ShiftedArrays v2.0.0
[992d4aef] Showoff v1.0.3
[777ac1f9] SimpleBufferStream v1.1.0
- [727e6d20] SimpleNonlinearSolve v0.1.23
+⌃ [727e6d20] SimpleNonlinearSolve v0.1.23
[699a6c99] SimpleTraits v0.9.4
[ce78b400] SimpleUnPack v1.1.0
[66db9d55] SnoopPrecompile v1.0.3
[b85f4697] SoftGlobalScope v1.1.0
[a2af1166] SortingAlgorithms v1.2.0
- [47a9eef4] SparseDiffTools v2.9.2
+⌃ [47a9eef4] SparseDiffTools v2.9.2
[e56a9233] Sparspak v0.3.9
[276daf66] SpecialFunctions v2.3.1
[928aab9d] SpecialMatrices v3.0.0
@@ -449,7 +537,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[a759f4b9] TimerOutputs v0.5.23
[c751599d] ToeplitzMatrices v0.8.2
[0796e94c] Tokenize v0.5.25
- [3bb67fe8] TranscodingStreams v0.10.1
+ [3bb67fe8] TranscodingStreams v0.10.2
[a2a6695c] TreeViews v0.3.0
[d5829a12] TriangularSolve v0.1.20
[410a4b4d] Tricks v0.1.8
@@ -468,7 +556,6 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[1b915085] WinReg v1.0.0
[ddb6d928] YAML v0.4.9
[c2297ded] ZMQ v1.2.2
- [700de1a5] ZygoteRules v0.2.4
[0518478a] deSolveDiffEq v0.1.1
[6e34b625] Bzip2_jll v1.0.8+0
[83423d85] Cairo_jll v1.16.1+1
@@ -507,7 +594,7 @@ Status `~/.julia/dev/ProbNumDiffEq/benchmarks/Manifest.toml`
[efe28fd5] OpenSpecFun_jll v0.5.5+0
[91d4177d] Opus_jll v1.3.2+0
[30392449] Pixman_jll v0.42.2+0
- [c0090381] Qt6Base_jll v6.5.2+2
+⌃ [c0090381] Qt6Base_jll v6.5.2+2
[f50d1b31] Rmath_jll v0.4.0+0
⌅ [fb77eaff] Sundials_jll v5.2.1+0
[a44049a8] Vulkan_Loader_jll v1.3.243+0
diff --git a/docs/src/initialization.md b/docs/src/initialization.md
index 2c45b74ed..9c7392146 100644
--- a/docs/src/initialization.md
+++ b/docs/src/initialization.md
@@ -63,6 +63,8 @@ See also [[1]](@ref initrefs).
```@docs
TaylorModeInit
+ForwardDiffInit
+SimpleInit
ClassicSolverInit
```
diff --git a/src/ProbNumDiffEq.jl b/src/ProbNumDiffEq.jl
index f8fc33edb..7cb7b99a0 100644
--- a/src/ProbNumDiffEq.jl
+++ b/src/ProbNumDiffEq.jl
@@ -58,7 +58,7 @@ include("diffusions.jl")
export FixedDiffusion, DynamicDiffusion, FixedMVDiffusion, DynamicMVDiffusion
include("initialization/common.jl")
-export TaylorModeInit, ClassicSolverInit
+export TaylorModeInit, ClassicSolverInit, SimpleInit, ForwardDiffInit
include("algorithms.jl")
export EK0, EK1
@@ -69,7 +69,8 @@ include("caches.jl")
include("checks.jl")
-include("initialization/taylormode.jl")
+include("initialization/simpleinit.jl")
+include("initialization/autodiffinit.jl")
include("initialization/classicsolverinit.jl")
include("solution.jl")
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 740df81dd..08d42e19b 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -8,7 +8,7 @@ abstract type AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm end
smooth=true,
prior=IWP(order),
diffusionmodel=DynamicDiffusion(),
- initialization=TaylorModeInit())
+ initialization=TaylorModeInit(prior.num_derivatives))
**Gaussian ODE filter with zeroth-order vector field linearization.**
@@ -49,7 +49,7 @@ EK0(;
prior=IWP(order),
diffusionmodel=DynamicDiffusion(),
smooth=true,
- initialization=TaylorModeInit(),
+ initialization=TaylorModeInit(prior.num_derivatives),
) = EK0(prior, diffusionmodel, smooth, initialization)
_unwrap_val(::Val{B}) where {B} = B
@@ -60,7 +60,7 @@ _unwrap_val(B) = B
smooth=true,
prior=IWP(order),
diffusionmodel=DynamicDiffusion(),
- initialization=TaylorModeInit(),
+ initialization=TaylorModeInit(prior.num_derivatives),
kwargs...)
**Gaussian ODE filter with first-order vector field linearization.**
@@ -103,7 +103,7 @@ EK1(;
prior::PT=IWP(order),
diffusionmodel::DT=DynamicDiffusion(),
smooth=true,
- initialization::IT=TaylorModeInit(),
+ initialization::IT=TaylorModeInit(prior.num_derivatives),
chunk_size=Val{0}(),
autodiff=Val{true}(),
diff_type=Val{:forward},
diff --git a/src/initialization/autodiffinit.jl b/src/initialization/autodiffinit.jl
new file mode 100644
index 000000000..3a7920fff
--- /dev/null
+++ b/src/initialization/autodiffinit.jl
@@ -0,0 +1,119 @@
+function initial_update!(integ, cache, init::AutodiffInitializationScheme)
+ @unpack u, f, p, t = integ
+ @unpack d, q, q, x, Proj = cache
+ D = d * (q + 1)
+
+ @unpack x_tmp, K1, C_Dxd, C_DxD, C_dxd, measurement = cache
+ if size(K1, 2) != d
+ K1 = K1[:, 1:d]
+ end
+
+ if f isa ODEFunction &&
+ f.f isa SciMLBase.FunctionWrappersWrappers.FunctionWrappersWrapper
+ f = ODEFunction(SciMLBase.unwrapped_f(f), mass_matrix=f.mass_matrix)
+ end
+
+ f_derivatives = get_derivatives(init, u, f, p, t)
+ integ.stats.nf += init.order
+ @assert length(f_derivatives) == init.order + 1
+
+ # This is hacky and should definitely be removed. But it also works so 🤷
+ MM = if f.mass_matrix isa UniformScaling
+ f.mass_matrix
+ else
+ _MM = copy(f.mass_matrix)
+ if any(iszero.(diag(_MM)))
+ _MM = typeof(promote(_MM[1], 1e-20)[1]).(_MM)
+ _MM .+= 1e-20I(d)
+ end
+ _MM
+ end
+
+ for (o, df) in zip(0:q, f_derivatives)
+ if f isa DynamicalODEFunction
+ @assert df isa ArrayPartition
+ df = df[2, :]
+ end
+
+ df = view(df, :)
+
+ H = if o == 0
+ Proj(o)
+ else
+ MM * Proj(o)
+ end
+ init_condition_on!(x, H, df, cache)
+ end
+end
+
+"""
+ Compute initial derivatives of an IIP ODEProblem with TaylorIntegration.jl
+"""
+function get_derivatives(
+ init::TaylorModeInit, u, f::SciMLBase.AbstractODEFunction{true}, p, t)
+ q = init.order
+ tT = Taylor1(typeof(t), q)
+ tT[0] = t
+ uT = similar(u, Taylor1{eltype(u)})
+ @inbounds @simd ivdep for i in eachindex(u)
+ uT[i] = Taylor1(u[i], q)
+ end
+ duT = zero(uT)
+ uauxT = similar(uT)
+ TaylorIntegration.jetcoeffs!(f, tT, uT, duT, uauxT, p)
+ # return hcat([evaluate.(differentiate.(uT, i)) for i in 0:q]...)'
+ return [evaluate.(differentiate.(uT, i)) for i in 0:q]
+end
+
+function get_derivatives(
+ init::ForwardDiffInit, u, f::SciMLBase.AbstractODEFunction{true}, p, t)
+ q = init.order
+ _f(u) = (du = copy(u); f(du, u, p, t); du)
+
+ out = [u]
+ push!(out, _f(u))
+
+ f_n = _f
+ for _ in 2:q
+ f_n = forwarddiff_oop_vectorfield_derivative_iteration(f_n, _f)
+ push!(out, f_n(u))
+ end
+
+ return out
+end
+
+function forwarddiff_oop_vectorfield_derivative_iteration(f_n, f_0)
+ function df(u)
+ J = ForwardDiff.jacobian(f_n, u)
+ return J * f_0(u)
+ end
+ return df
+end
+
+# Curently unused but potentially a future upgrade to the version above:
+function iip_forwarddiff_get_derivatives!(
+ out, init::ForwardDiffInit, u, f::SciMLBase.AbstractODEFunction{true}, p, t)
+ q = init.order
+ d = length(u)
+ _f(du, u) = f(du, u, p, t)
+
+ out[1:d] .= u
+ @views _f(out[d+1:2d], u)
+
+ f_n = _f
+ for o in 2:q
+ f_n = forwarddiff_iip_vectorfield_derivative_iteration(f_n, _f)
+ @views f_n(out[o*d+1:(o+1)*d], u)
+ end
+
+ return out
+end
+
+function forwarddiff_iip_vectorfield_derivative_iteration(f_n, f_0)
+ function df(du, u)
+ J = ForwardDiff.jacobian(f_n, du, u)
+ f_0(du, u)
+ _matmul!(du, J, du)
+ end
+ return df
+end
diff --git a/src/initialization/common.jl b/src/initialization/common.jl
index d5c9fa0ba..ed8d5accf 100644
--- a/src/initialization/common.jl
+++ b/src/initialization/common.jl
@@ -1,9 +1,20 @@
abstract type InitializationScheme end
+abstract type AutodiffInitializationScheme <: InitializationScheme end
"""
- TaylorModeInit()
+ SimpleInit()
-Exact initialization via Taylor-mode automatic differentiation.
+Simple initialization, only with the given initial value and derivative.
+
+The remaining derivatives are set to zero with unit covariance (unless specified otherwise
+by setting a custom [`FixedDiffusion`](@ref)).
+"""
+struct SimpleInit <: InitializationScheme end
+
+"""
+ TaylorModeInit(order)
+
+Exact initialization via Taylor-mode automatic differentiation up to order `order`.
**This is the recommended initialization method!**
@@ -13,12 +24,45 @@ via Taylor-mode automatic differentiation.
In some special cases it can happen that TaylorIntegration.jl is incompatible with the
given problem (typically because the problem definition does not allow for elements of type
- `Taylor`). If this happens, try [`ClassicSolverInit`](@ref).
+ `Taylor`). If this happens, try one of [`SimpleInit`](@ref), [`ForwardDiffInit`](@ref)
+(for low enough orders), [`ClassicSolverInit`](@ref).
# References
* [kraemer20stableimplementation](@cite) Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
"""
-struct TaylorModeInit <: InitializationScheme end
+struct TaylorModeInit <: AutodiffInitializationScheme
+ order::Int64
+ TaylorModeInit(order::Int64) = begin
+ if order < 1
+ throw(ArgumentError("order must be >= 1"))
+ end
+ new(order)
+ end
+end
+TaylorModeInit() = begin
+ throw(ArgumentError("order must be specified"))
+end
+
+"""
+ ForwardDiffInit(order)
+
+Exact initialization via ForwardDiff.jl up to order `order`.
+
+**Warning:** This does not scale well to high orders!
+For orders > 3, [`TaylorModeInit`](@ref) most likely performs better.
+"""
+struct ForwardDiffInit <: AutodiffInitializationScheme
+ order::Int64
+ ForwardDiffInit(order::Int64) = begin
+ if order < 1
+ throw(ArgumentError("order must be >= 1"))
+ end
+ new(order)
+ end
+end
+ForwardDiffInit() = begin
+ throw(ArgumentError("order must be specified"))
+end
"""
ClassicSolverInit(; alg=OrdinaryDiffEq.Tsit5(), init_on_ddu=false)
@@ -44,10 +88,14 @@ optionally the second derivative can also be set via automatic differentiation b
* [kraemer20stableimplementation](@cite) Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
* [schober16probivp](@cite) Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
"""
-Base.@kwdef struct ClassicSolverInit{ALG} <: InitializationScheme
- alg::ALG = Tsit5()
- init_on_ddu::Bool = false
+struct ClassicSolverInit{ALG} <: InitializationScheme
+ alg::ALG
+ init_on_ddu::Bool
end
+ClassicSolverInit(alg::DiffEqBase.AbstractODEAlgorithm=AutoVern7(Rodas4())) =
+ ClassicSolverInit(; alg, init_on_ddu=false)
+ClassicSolverInit(; alg=AutoVern7(Rodas4()), init_on_ddu=false) =
+ ClassicSolverInit(alg, init_on_ddu)
"""
initial_update!(integ, cache[, init::InitializationScheme])
diff --git a/src/initialization/simpleinit.jl b/src/initialization/simpleinit.jl
new file mode 100644
index 000000000..e410b5adb
--- /dev/null
+++ b/src/initialization/simpleinit.jl
@@ -0,0 +1,35 @@
+function initial_update!(integ, cache, init::SimpleInit)
+ @unpack u, f, p, t = integ
+ @unpack x, d, Proj = cache
+ du = integ.uprev
+
+ if f isa ODEFunction &&
+ f.f isa SciMLBase.FunctionWrappersWrappers.FunctionWrappersWrapper
+ f = ODEFunction(SciMLBase.unwrapped_f(f), mass_matrix=f.mass_matrix)
+ end
+
+ # This is hacky and should definitely be removed. But it also works so 🤷
+ MM = if f.mass_matrix isa UniformScaling
+ f.mass_matrix
+ else
+ _MM = copy(f.mass_matrix)
+ if any(iszero.(diag(_MM)))
+ _MM = typeof(promote(_MM[1], 1e-20)[1]).(_MM)
+ _MM .+= 1e-20I(d)
+ end
+ _MM
+ end
+
+ f(du, u, p, t)
+ integ.stats.nf += 1
+
+ if f isa DynamicalODEFunction
+ @assert u isa ArrayPartition
+ u = u[2, :]
+ @assert du isa ArrayPartition
+ du = du[2, :]
+ end
+
+ init_condition_on!(x, Proj(0), view(u, :), cache)
+ init_condition_on!(x, MM * Proj(1), view(du, :), cache)
+end
diff --git a/src/initialization/taylormode.jl b/src/initialization/taylormode.jl
deleted file mode 100644
index 57e873453..000000000
--- a/src/initialization/taylormode.jl
+++ /dev/null
@@ -1,63 +0,0 @@
-function initial_update!(integ, cache, init::TaylorModeInit)
- @unpack u, f, p, t = integ
- @unpack d, q, q, x, Proj = cache
- D = d * (q + 1)
-
- @unpack x_tmp, K1, C_Dxd, C_DxD, C_dxd, measurement = cache
- if size(K1, 2) != d
- K1 = K1[:, 1:d]
- end
-
- if f isa ODEFunction &&
- f.f isa SciMLBase.FunctionWrappersWrappers.FunctionWrappersWrapper
- f = ODEFunction(SciMLBase.unwrapped_f(f), mass_matrix=f.mass_matrix)
- end
-
- f_derivatives = taylormode_get_derivatives(u, f, p, t, q)
- integ.stats.nf += q
- @assert length(0:q) == length(f_derivatives)
-
- # This is hacky and should definitely be removed. But it also works so 🤷
- MM = if f.mass_matrix isa UniformScaling
- f.mass_matrix
- else
- _MM = copy(f.mass_matrix)
- if any(iszero.(diag(_MM)))
- _MM = typeof(promote(_MM[1], 1e-20)[1]).(_MM)
- _MM .+= 1e-20I(d)
- end
- _MM
- end
-
- for (o, df) in zip(0:q, f_derivatives)
- if f isa DynamicalODEFunction
- @assert df isa ArrayPartition
- df = df[2, :]
- end
-
- df = view(df, :)
-
- H = if o == 0
- Proj(o)
- else
- MM * Proj(o)
- end
- init_condition_on!(x, H, df, cache)
- end
-end
-
-"""
- Compute initial derivatives of an IIP ODEProblem with TaylorIntegration.jl
-"""
-function taylormode_get_derivatives(u, f::SciMLBase.AbstractODEFunction{true}, p, t, q)
- tT = Taylor1(typeof(t), q)
- tT[0] = t
- uT = similar(u, Taylor1{eltype(u)})
- @inbounds @simd ivdep for i in eachindex(u)
- uT[i] = Taylor1(u[i], q)
- end
- duT = zero(uT)
- uauxT = similar(uT)
- TaylorIntegration.jetcoeffs!(f, tT, uT, duT, uauxT, p)
- return [evaluate.(differentiate.(uT, i)) for i in 0:q]
-end
diff --git a/test/complexity.jl b/test/complexity.jl
index a20a46d41..31482dea5 100644
--- a/test/complexity.jl
+++ b/test/complexity.jl
@@ -44,7 +44,7 @@ using Test, SafeTestsets
lr_ek1 = linregress(log.(dims_ek1), log.(times_ek1))
@test_skip slope(lr_ek1)[1] ≈ 2 atol = 0.2
# This is what we would actually expect, not sure what's going wrong:
- @test_broken slope(lr_ek1)[1] ≈ 3 atol = 0.1
+ @test_skip slope(lr_ek1)[1] ≈ 3 atol = 0.1
end
@testset "Order 3 + Taylor-init + no smoothing" begin
@@ -73,7 +73,7 @@ using Test, SafeTestsets
lr_ek1 = linregress(log.(dims_ek1), log.(times_ek1))
@test_skip slope(lr_ek1)[1] ≈ 2 atol = 0.5
# This is what we would actually expect, not sure what's going wrong:
- @test_broken slope(lr_ek1)[1] ≈ 3 atol = 0.1
+ @test_skip slope(lr_ek1)[1] ≈ 3 atol = 0.1
end
@testset "Order 3 with smoothing and everyting" begin
@@ -99,6 +99,6 @@ using Test, SafeTestsets
lr_ek1 = linregress(log.(dims_ek1), log.(times_ek1))
@test_skip slope(lr_ek1)[1] ≈ 2 atol = 0.2
# This is what we would actually expect, not sure what's going wrong:
- @test_broken slope(lr_ek1)[1] ≈ 3 atol = 0.1
+ @test_skip slope(lr_ek1)[1] ≈ 3 atol = 0.1
end
end
diff --git a/test/correctness.jl b/test/correctness.jl
index 5a6d593cc..22e52c508 100644
--- a/test/correctness.jl
+++ b/test/correctness.jl
@@ -17,6 +17,7 @@ CONSTANT_ALGS = (
EK0(order=3, smooth=false, diffusionmodel=FixedMVDiffusion()) => 1e-7,
EK0(order=3, smooth=false, diffusionmodel=DynamicMVDiffusion()) => 1e-8,
EK0(order=3, smooth=false, initialization=ClassicSolverInit()) => 1e-7,
+ EK0(order=3, smooth=false, initialization=SimpleInit()) => 1e-5,
EK0(
order=3,
smooth=false,
@@ -28,6 +29,7 @@ CONSTANT_ALGS = (
EK1(order=5, smooth=false) => 1e-11,
EK1(order=3, smooth=false, diffusionmodel=FixedDiffusion()) => 1e-8,
EK1(order=3, smooth=false, initialization=ClassicSolverInit()) => 1e-7,
+ EK1(order=3, smooth=false, initialization=SimpleInit()) => 1e-5,
# smoothing
EK0(order=3, smooth=true) => 1e-8,
EK0(order=3, smooth=true, diffusionmodel=FixedDiffusion()) => 2e-8,
@@ -43,12 +45,14 @@ ADAPTIVE_ALGS = (
EK0(order=8) => 2e-5,
EK0(order=3, diffusionmodel=DynamicMVDiffusion()) => 5e-5,
EK0(order=3, initialization=ClassicSolverInit()) => 5e-5,
+ EK0(order=3, initialization=SimpleInit()) => 1e-4,
EK0(order=3, diffusionmodel=DynamicMVDiffusion(), initialization=ClassicSolverInit()) => 4e-5,
EK1(order=2) => 2e-5,
EK1(order=3) => 1e-5,
EK1(order=5) => 1e-6,
EK1(order=8) => 5e-6,
EK1(order=3, initialization=ClassicSolverInit()) => 1e-5,
+ EK1(order=3, initialization=SimpleInit()) => 1e-4,
)
PROBS = (
diff --git a/test/mass_matrix.jl b/test/mass_matrix.jl
index ca1a6cfa8..ea74a4bda 100644
--- a/test/mass_matrix.jl
+++ b/test/mass_matrix.jl
@@ -66,7 +66,16 @@ end
f = ODEFunction(rober, mass_matrix=M)
prob = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e-2), (0.04, 3e7, 1e4))
- ref = solve(prob, EK1(order=3))
- sol = solve(prob, RadauIIA5())
+ ref = solve(prob, RadauIIA5())
+ sol = solve(prob, EK1(order=3))
+ @test sol[end] ≈ ref[end] rtol = 1e-8
+
+ sol = solve(prob, EK1(order=3, initialization=ForwardDiffInit(3)))
+ @test sol[end] ≈ ref[end] rtol = 1e-8
+
+ sol = solve(prob, EK1(order=3, initialization=ClassicSolverInit(RadauIIA5())))
+ @test sol[end] ≈ ref[end] rtol = 1e-8
+
+ sol = solve(prob, EK1(order=3, initialization=SimpleInit()))
@test sol[end] ≈ ref[end] rtol = 1e-8
end
diff --git a/test/secondorderode.jl b/test/secondorderode.jl
index 5addf4be9..83cc7b325 100644
--- a/test/secondorderode.jl
+++ b/test/secondorderode.jl
@@ -47,19 +47,22 @@ end
end
end
-@testset "ClassicSolverInit for SecondOrderODEProblems" begin
+_name(structinstance) = typeof(structinstance).name.wrapper
+@testset "Diffusion: $(_name(DIFFUSION))" for DIFFUSION in (
+ FixedDiffusion(), DynamicDiffusion(), FixedMVDiffusion(), DynamicMVDiffusion())
@test_nowarn solve(
prob_iip,
- EK1(initialization=ClassicSolverInit(), smooth=false),
+ EK0(diffusionmodel=FixedDiffusion(), smooth=false),
save_everystep=false,
dense=false,
)
end
-@testset "Fixed Diffusion" begin
+@testset "Initialization: $(_name(INIT))" for INIT in (
+ TaylorModeInit(3), ForwardDiffInit(3), SimpleInit(), ClassicSolverInit())
@test_nowarn solve(
prob_iip,
- EK0(diffusionmodel=FixedDiffusion(), smooth=false),
+ EK1(initialization=INIT, smooth=false),
save_everystep=false,
dense=false,
)
diff --git a/test/state_init.jl b/test/state_init.jl
index 9fe8be214..4330e5804 100644
--- a/test/state_init.jl
+++ b/test/state_init.jl
@@ -33,23 +33,26 @@ import ODEProblemLibrary: prob_ode_fitzhughnagumo, prob_ode_pleiades
f!(du, u, p, t) = (du .= f(u, p, t))
prob = ODEProblem{true,true}(f!, u0, tspan)
- @testset "`taylormode_get_derivatives`" begin
- dfs = ProbNumDiffEq.taylormode_get_derivatives(
+ @testset "Exact `get_dervatives`: $INIT" for INIT in (
+ TaylorModeInit(q), ForwardDiffInit(q))
+ dfs = ProbNumDiffEq.get_derivatives(
+ INIT,
prob.u0,
prob.f,
prob.p,
prob.tspan[1],
- q,
)
@test length(dfs) == q + 1
@test true_init_states ≈ vcat(dfs...)
end
- @testset "Taylormode: `initial_update!`" begin
- integ = init(prob, EK0(order=q))
- ProbNumDiffEq.initial_update!(integ, integ.cache, TaylorModeInit())
- x = integ.cache.x
- @test reshape(x.μ, :, 2)'[:] ≈ true_init_states
+ _name(structinstance) = typeof(structinstance).name.wrapper
+ @testset "`init` and `initial_update!` with $(_name(INIT))" for INIT in (
+ TaylorModeInit(q), ForwardDiffInit(q), SimpleInit())
+ _q = INIT isa SimpleInit ? 1 : INIT.order
+ integ = init(prob, EK0(order=q, initialization=INIT))
+ dus = integ.cache.x.μ
+ @test reshape(dus, :, 2)'[1:d*(_q+1)] ≈ true_init_states[1:d*(_q+1)]
end
@testset "Low-order exact init via ClassiSolverInit: `initial_update!`" begin
@@ -121,3 +124,18 @@ end
end
end
end
+
+@testset "Errors" begin
+ @test_throws ArgumentError ForwardDiffInit()
+ @test_nowarn ForwardDiffInit(1)
+ @test_throws ArgumentError ForwardDiffInit(0)
+ @test_throws ArgumentError ForwardDiffInit(-1)
+ @test_throws ArgumentError TaylorModeInit()
+ @test_nowarn TaylorModeInit(1)
+ @test_throws ArgumentError TaylorModeInit(0)
+ @test_throws ArgumentError TaylorModeInit(-1)
+ @test_nowarn SimpleInit()
+ @test_nowarn ClassicSolverInit()
+ @test_nowarn ClassicSolverInit(Tsit5())
+ @test_throws MethodError ClassicSolverInit(3)
+end
diff --git a/test/stats.jl b/test/stats.jl
index e9ec43800..2362f5e1f 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -1,4 +1,5 @@
using ProbNumDiffEq
+using OrdinaryDiffEq
using Test
using LinearAlgebra
@@ -6,7 +7,7 @@ const q = 3
@testset "stats.nf testing $alg" for alg in (
EK0(prior=IWP(q), smooth=false),
- EK0(prior=IWP(q), smooth=false, initialization=ClassicSolverInit()),
+ EK0(prior=IWP(q), smooth=false, initialization=ClassicSolverInit(alg=Tsit5())),
EK1(prior=IWP(q), smooth=false),
EK1(prior=IWP(q), smooth=false, autodiff=false),
EK1(prior=IOUP(q, -1), smooth=false),