Skip to content

Commit

Permalink
Update Lotka-Volterra and add calibration plots
Browse files Browse the repository at this point in the history
  • Loading branch information
nathanaelbosch committed Nov 13, 2023
1 parent e412f4e commit 3266335
Show file tree
Hide file tree
Showing 19 changed files with 6,466 additions and 4,819 deletions.
1 change: 1 addition & 0 deletions benchmarks/Project.toml
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
[deps]
DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
Expand Down
118 changes: 43 additions & 75 deletions benchmarks/lotkavolterra.jmd
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ Benchmark adapted from
[SciMLBenchmarks.jl](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/NonStiffODE/LotkaVolterra_wpd/).

```julia
using LinearAlgebra, Statistics
using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Plots
using ProbNumDiffEq

Expand Down Expand Up @@ -48,7 +48,7 @@ plot(test_sol, title="Lotka-Volterra Solution", legend=false, xticks=:auto)

## EK0 accross orders

Final timepoint only:
### Final timepoint only
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
Expand All @@ -72,14 +72,12 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
verbose = false,
)

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

Discrete time-series errors (l2):
### Discrete time-series errors (l2)
```julia
DENSE = true;
SAVE_EVERYSTEP = true;
Expand All @@ -103,37 +101,19 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = true,
dense_errors = true,
verbose = false,
error_estimate = :l2,
)

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

Interpolation errors (L2):
### Interpolation errors (L2)
```julia
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = true,
dense_errors = true,
verbose = false,
error_estimate = :L2,
)

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

## EK1 accross orders

Last timepoint only:
### Final timepoint only
```julia
DENSE = false;
SAVE_EVERYSTEP = false;
Expand All @@ -157,14 +137,12 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
verbose = false,
)

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

Discrete time-series errors (l2):
### Discrete time-series errors (l2)
```julia
DENSE = true;
SAVE_EVERYSTEP = true;
Expand All @@ -188,35 +166,17 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = true,
dense_errors = true,
verbose = false,
error_estimate = :l2,
)

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

Interpolation errors (L2):
### Interpolation errors (L2)
```julia
wp = WorkPrecisionSet(
prob, abstols, reltols, setups;
names = labels,
appxsol = test_sol,
dense = DENSE,
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = true,
dense_errors = true,
verbose = false,
error_estimate = :L2,
)

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

## EK0 vs. EK1
## EK0 vs. EK1: Work-Precision

Final timepoint only:
```julia
Expand Down Expand Up @@ -248,9 +208,6 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
dense_errors = false,
verbose = false,
)

plot(wp, color=[1 1 1 1 2 2 2 2])
Expand Down Expand Up @@ -286,32 +243,51 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
dense_errors = true,
verbose = true,
error_estimate = :L2,
)

plot(wp, color=[1 1 1 1 2 2 2 2])
```

## EK0: Diffusion Comparison
## EK0 vs. EK1: Calibration
Final time-point:
```julia
plot(wp, x=:final, y=:chi2_final, color=[1 1 1 1 2 2 2 2])

# Should be distributed according to a Chi-squared distribution:
low, high, mid = quantile(Chisq(2), [0.01, 0.99]), mean(Chisq(2))
hline!([low, high], linestyle=:dash, color=:black, label="")
hline!([mid], linestyle=:solid, color=:black, label="")
```

Discrete time-series:
```julia
plot(wp, x=:l2, y=:chi2_steps, color=[1 1 1 1 2 2 2 2])
hline!([low, high], linestyle=:dash, color=:black, label="")
hline!([mid], linestyle=:solid, color=:black, label="")
```

Interpolation:
```julia
plot(wp, x=:L2, y=:chi2_interp, color=[1 1 1 1 2 2 2 2])
hline!([low, high], linestyle=:dash, color=:black, label="")
hline!([mid], linestyle=:solid, color=:black, label="")
```

## Comparison of the different diffusion models

### EK0

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

_setups = [
"EK0(2) Dynamic" => Dict(:alg => EK0(order=2, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK0(3) Dynamic" => Dict(:alg => EK0(order=3, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK0(0) Dynamic" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK0(2) Fixed" => Dict(:alg => EK0(order=2, smooth=DENSE, diffusionmodel=FixedDiffusion()))
"EK0(5) Dynamic" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK0(3) Fixed" => Dict(:alg => EK0(order=3, smooth=DENSE, diffusionmodel=FixedDiffusion()))
"EK0(5) Fixed" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=FixedDiffusion()))
"EK0(2) DynamicMV" => Dict(:alg => EK0(order=2, smooth=DENSE, diffusionmodel=DynamicMVDiffusion()))
"EK0(3) DynamicMV" => Dict(:alg => EK0(order=3, smooth=DENSE, diffusionmodel=DynamicMVDiffusion()))
"EK0(0) DynamicMV" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=DynamicMVDiffusion()))
"EK0(2) FixedMV" => Dict(:alg => EK0(order=2, smooth=DENSE, diffusionmodel=FixedMVDiffusion()))
"EK0(5) DynamicMV" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=DynamicMVDiffusion()))
"EK0(3) FixedMV" => Dict(:alg => EK0(order=3, smooth=DENSE, diffusionmodel=FixedMVDiffusion()))
"EK0(5) FixedMV" => Dict(:alg => EK0(order=5, smooth=DENSE, diffusionmodel=FixedMVDiffusion()))
]
Expand All @@ -330,24 +306,20 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
verbose = false,
)

plot(wp, color=[2 2 2 3 3 3 4 4 4 5 5 5])
```

## EK1: DynamicDiffusion vs FixedDiffusion
### EK1

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

_setups = [
"EK1(2) Dynamic" => Dict(:alg => EK1(order=2, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK1(3) Dynamic" => Dict(:alg => EK1(order=3, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK1(5) Dynamic" => Dict(:alg => EK1(order=5, smooth=DENSE, diffusionmodel=DynamicDiffusion()))
"EK1(2) Fixed" => Dict(:alg => EK1(order=2, smooth=DENSE, diffusionmodel=FixedDiffusion()))
"EK1(3) Fixed" => Dict(:alg => EK1(order=3, smooth=DENSE, diffusionmodel=FixedDiffusion()))
"EK1(5) Fixed" => Dict(:alg => EK1(order=5, smooth=DENSE, diffusionmodel=FixedDiffusion()))
]
Expand All @@ -366,8 +338,6 @@ wp = WorkPrecisionSet(
save_everystep = SAVE_EVERYSTEP,
numruns = 10,
maxiters = Int(1e7),
timeseries_errors = false,
verbose = false,
)

plot(wp, color=[2 2 2 3 3 3])
Expand Down Expand Up @@ -403,8 +373,6 @@ for o in orders
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), title = "Order $o")
Expand Down
10 changes: 5 additions & 5 deletions benchmarks/runall.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@ ENV["GKSwstype"] = "nul"

FILES = [
"lotkavolterra.jmd",
"vanderpol.jmd",
"hodgkinhuxley.jmd",
"rober.jmd",
"pleiades.jmd",
"multi-language-wrappers.jmd",
# "vanderpol.jmd",
# "hodgkinhuxley.jmd",
# "rober.jmd",
# "pleiades.jmd",
# "multi-language-wrappers.jmd",
]

filedir = @__DIR__
Expand Down
Loading

0 comments on commit 3266335

Please sign in to comment.