Skip to content

Commit

Permalink
Made the benchmarks nicer once more. Soon to be merged!
Browse files Browse the repository at this point in the history
  • Loading branch information
@nathanaelbosch
nathanaelbosch committed Nov 22, 2023
1 parent 2317147 commit de00edb
Show file tree
Hide file tree
Showing 55 changed files with 13,721 additions and 12,400 deletions.
31 changes: 30 additions & 1 deletion benchmarks/hodgkinhuxley.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
- [Results are similar for fixed time steps.](@ref hh_fixed_steps)


```julia
```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, SciMLBase, OrdinaryDiffEq, Plots, SimpleUnPack
using ProbNumDiffEq
Expand Down Expand Up @@ -111,6 +111,33 @@ wp = WorkPrecisionSet(
)

plot(wp, title="Adaptive steps - no smoothing", color=colors)


_ref_setups = [
"Tsit5" => Dict(:alg => Tsit5())
"Vern7" => Dict(:alg => Vern7())
"RadauIIA5" => Dict(:alg => RadauIIA5())
]
ref_labels = first.(_ref_setups)
ref_setups = last.(_ref_setups)
ref_wp_final = WorkPrecisionSet(
prob, abstols ./ 1000, reltols ./ 1000, ref_setups;
names = ref_labels,
appxsol = test_sol,
dense = false,
save_everystep = false,
maxiters = Int(1e7),
)
ref_wp_dense = WorkPrecisionSet(
prob, abstols ./ 1000, reltols ./ 1000, ref_setups;
names = ref_labels,
appxsol = test_sol,
dense = true,
save_everystep = true,
maxiters = Int(1e7),
)

plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

### With smoothing
Expand Down Expand Up @@ -146,13 +173,15 @@ wp = WorkPrecisionSet(
)

plot(wp, title="Adaptive steps - with smoothing", color=colors)
plot!(ref_wp_dense, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

```@raw html
<details><summary>Interoplation errors (L2):</summary>
```
```julia
plot(wp, x=:L2, title="Adaptive steps - with smoothing", color=colors)
plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand Down
52 changes: 42 additions & 10 deletions benchmarks/lotkavolterra.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
Benchmark adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/NonStiffODE/LotkaVolterra_wpd/).

```julia
```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
using ProbNumDiffEq
Expand All @@ -32,6 +32,14 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)

function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
```

```julia
Expand Down Expand Up @@ -75,6 +83,32 @@ wp = WorkPrecisionSet(
)

plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))

_ref_setups = [
"Tsit5" => Dict(:alg => Tsit5())
"Vern7" => Dict(:alg => Vern7())
"RadauIIA5" => Dict(:alg => RadauIIA5())
]
ref_labels = first.(_ref_setups)
ref_setups = last.(_ref_setups)
ref_wp_final = WorkPrecisionSet(
prob, abstols, reltols, ref_setups;
names = ref_labels,
appxsol = test_sol,
dense = false,
save_everystep = false,
maxiters = Int(1e7),
)
ref_wp_dense = WorkPrecisionSet(
prob, abstols, reltols, ref_setups;
names = ref_labels,
appxsol = test_sol,
dense = true,
save_everystep = true,
maxiters = Int(1e7),
)

plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

```@raw html
Expand Down Expand Up @@ -105,6 +139,7 @@ wp = WorkPrecisionSet(
)

plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
plot!(ref_wp_dense, x=:l2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand All @@ -115,6 +150,7 @@ plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
```
```julia
plot(wp, x=:L2, palette=Plots.palette([:blue, :red], length(_setups)))
plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand Down Expand Up @@ -147,6 +183,7 @@ wp = WorkPrecisionSet(
)

plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))
plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

```@raw html
Expand Down Expand Up @@ -177,6 +214,7 @@ wp = WorkPrecisionSet(
)

plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
plot!(ref_wp_dense, x=:l2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand All @@ -187,6 +225,7 @@ plot(wp, x=:l2, palette=Plots.palette([:blue, :red], length(_setups)))
```
```julia
plot(wp, x=:L2, palette=Plots.palette([:blue, :red], length(_setups)))
plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand Down Expand Up @@ -225,6 +264,7 @@ wp = WorkPrecisionSet(
)

plot(wp, color=[1 1 1 1 2 2 2 2])
plot!(ref_wp_final, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

```@raw html
Expand Down Expand Up @@ -261,6 +301,7 @@ wp = WorkPrecisionSet(
)

plot(wp, x=:L2, color=[1 1 1 1 2 2 2 2])
plot!(ref_wp_dense, x=:L2, color=:gray, alpha=0.7, linestyle=:dash)
```
```@raw html
</details>
Expand All @@ -269,15 +310,6 @@ plot(wp, x=:L2, color=[1 1 1 1 2 2 2 2])
## [`EK0`](@ref) vs. [`EK1`](@ref): Calibration
```julia
plot(wp, x=:final, y=:chi2_final, color=[1 1 1 1 2 2 2 2], yguide="Chi-squared (final)")

# Should be distributed according to a Chi-squared distribution:
function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
plot_chisq_interval!(2)
```

Expand Down
49 changes: 33 additions & 16 deletions benchmarks/orego.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,14 @@


!!! note "Summary"
TODO
- [**The `EK1` is able to solve mass-matrix DAEs.** To achieve low error, use order 4 or higher.](@ref orego_results)
- The order-to-error-tolerance heuristic holds: lower tolerance level ``\rightarrow`` higher order.


Adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DAE/OregoDAE/).

```julia
```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots
using ModelingToolkit
Expand All @@ -21,9 +22,15 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
```

### OREGO problem definition
function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
```

```julia
@variables t y1(t)=1.0 y2(t)=2.0 y3(t)=3.0
Expand All @@ -47,22 +54,22 @@ ref_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14)
plot(ref_sol, title="OREGO Solution", legend=false, xticks=:auto)
```

## EK1 accross orders
## [`EK1` across orders](@id orego_results)

```julia
DENSE = false;
SAVE_EVERYSTEP = false;

_setups = [
"EK1($order)" => Dict(:alg => EK1(order=order, smooth=DENSE))
for order in 2:4
for order in 2:6
]

labels = first.(_setups)
setups = last.(_setups)

abstols = 1.0 ./ 10.0 .^ (4:8)
reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = 1.0 ./ 10.0 .^ (6:10)
reltols = 1.0 ./ 10.0 .^ (3:7)

wp = WorkPrecisionSet(
mmprob, abstols, reltols, setups;
Expand All @@ -75,21 +82,31 @@ wp = WorkPrecisionSet(
)

plot(wp, palette=Plots.palette([:blue, :red], length(_setups)))


_ref_setups = [
"Rosenbrock23" => Dict(:alg => Rosenbrock23())
"Rodas4P" => Dict(:alg => Rodas4P())
"RadauIIA" => Dict(:alg => RadauIIA5())
]
ref_labels = first.(_ref_setups)
ref_setups = last.(_ref_setups)
ref_wp = WorkPrecisionSet(
mmprob, abstols, reltols, ref_setups;
names = ref_labels,
appxsol = ref_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
maxiters = Int(1e7),
)
plot!(ref_wp, x=:final, color=:gray, alpha=0.7, linestyle=:dash)
```

### Calibration

```julia
plot(wp; x=:final, y=:chi2_final, yguide="Chi-squared (final)",
palette=Plots.palette([:blue, :red], length(_setups)))

function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
plot_chisq_interval!(3)
```

Expand Down
59 changes: 32 additions & 27 deletions benchmarks/pleiades.jmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@

!!! note "Summary"
Pleiades is a medium-dimensional, non-stiff, second-order ODE. We see that:
- The `EK0` is _much_ faster than the `EK1` as it scales linearly with the ODE dimension.
- If the problem is a second-order ODE, _implement it as a second-order ODE_!
- [**The `EK0` is _much_ faster than the `EK1` as it scales linearly with the ODE dimension.**](@ref pleiades_results)
- [**If the problem is a second-order ODE, _implement it as a second-order ODE_!**](@ref pleiades_results)


```julia
```julia, results="hidden"
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots, ODEInterfaceDiffEq
using ModelingToolkit
Expand All @@ -20,9 +20,15 @@ Plots.theme(
margin=5Plots.mm,
xticks=10.0 .^ (-16:1:16),
)
```

### Pleiades problem definition
function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
```

```julia
# first-order ODE
Expand Down Expand Up @@ -83,7 +89,7 @@ plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legen
scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7)
```

## EK0 vs EK1 & first-order vs. second-order
## [`EK0` vs `EK1` & first-order vs. second-order](@id pleiades_results)
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
Expand All @@ -97,9 +103,6 @@ _setups = [
"EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1)
"EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2)
"EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2)
"Classic: Tsit5" => Dict(:alg => Tsit5(), :prob_choice => 1)
"Classic: RadauIIA5" => Dict(:alg => RadauIIA5(), :prob_choice => 1)
"Classic: DPRKN6" => Dict(:alg => DPRKN6(), :prob_choice => 2)
]

labels = first.(_setups)
Expand All @@ -118,30 +121,32 @@ wp = WorkPrecisionSet(
maxiters = Int(1e7),
)

plot(wp,
color=[1 1 1 1 2 2 2 2 3 3 3],
linestyle=[:solid :solid :dash :dash :solid :solid :dash :dash :solid :solid :dash],
alpha=[1 1 1 1 1 1 1 1 0.5 0.5 0.5],
color = [1 1 1 1 2 2 2 2]
linestyle = [:solid :solid :dash :dash :solid :solid :dash :dash]
plot(wp; color, linestyle)

_ref_setups = [
"Classic: Tsit5" => Dict(:alg => Tsit5(), :prob_choice => 1)
"Classic: RadauIIA5" => Dict(:alg => RadauIIA5(), :prob_choice => 1)
"Classic: DPRKN6" => Dict(:alg => DPRKN6(), :prob_choice => 2)
]
ref_labels = first.(_ref_setups)
ref_setups = last.(_ref_setups)
ref_wp = WorkPrecisionSet(
probs, abstols ./ 1000, reltols ./ 1000, ref_setups;
names = ref_labels,
appxsol = ref_sols,
dense = false,
save_everystep = false,
maxiters = Int(1e7),
)
plot!(ref_wp, x=:final, color=:gray, alpha=0.7, linestyle=[:solid :solid :dash])
```

### Calibration

```julia
# We can only evaluate calibration for the PN solvers
_wp = WorkPrecisionSet(wp.wps[1:end-3], wp.N-3, wp.abstols, wp.reltols, wp.prob, wp.setups[1:end-3],
wp.names[1:end-3], wp.error_estimate, wp.numruns)
plot(_wp; x=:final, y=:chi2_final,
color=[1 1 1 1 2 2 2 2],
linestyle=[:solid :solid :dash :dash :solid :solid :dash :dash],
)
function plot_chisq_interval!(df, q=0.01)
dist = Chisq(df)
low, high, mid = quantile(dist, [q, 1-q])..., mean(dist)
hline!([low, high], linestyle=:dash, color=:black, label="",
fill_between=[high nothing], fillcolor=:green, fillalpha=0.15)
hline!([mid], linestyle=:solid, color=:black, label="")
end
plot(wp; x=:final, y=:chi2_final, color, linestyle)
plot_chisq_interval!(length(u0)*2)
```

Expand Down
Loading

0 comments on commit de00edb

Please sign in to comment.