From 64fef2275fde4faa5cd3bda6560a420d62752957 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Mon, 16 Oct 2023 08:07:59 +0000 Subject: [PATCH] build based on 06c9ea6 --- dev/.documenter-siteinfo.json | 2 +- dev/benchmarks/lotkavolterra/index.html | 2 +- .../multi-language-wrappers/index.html | 2 +- dev/benchmarks/pleiades/index.html | 2 +- dev/benchmarks/rober/index.html | 2 +- dev/benchmarks/vanderpol/index.html | 2 +- dev/diffusions/index.html | 4 +- dev/filtering/index.html | 8 +- dev/implementation/index.html | 2 +- dev/index.html | 2 +- dev/initialization/index.html | 4 +- dev/priors/index.html | 6 +- dev/references/index.html | 2 +- dev/solvers/index.html | 10 +- .../dae/{81b0c82b.svg => 23ce2a2c.svg} | 90 ++++++------ .../dae/{12e69d49.svg => fb9c5960.svg} | 134 +++++++++--------- dev/tutorials/dae/index.html | 6 +- .../{e22b23c3.svg => 136579f1.svg} | 72 +++++----- .../{5fcbcdb1.svg => 52e65c07.svg} | 76 +++++----- .../{27923a4c.svg => ac97650f.svg} | 64 ++++----- .../{373f8431.svg => e1d09d06.svg} | 64 ++++----- dev/tutorials/dynamical_odes/index.html | 10 +- .../{d844aea3.svg => 1c4f7b27.svg} | 74 +++++----- .../{d7a376a7.svg => 9b15c58f.svg} | 98 ++++++------- .../{5e6b2b1f.svg => df738690.svg} | 102 ++++++------- .../{5f074e4e.svg => e9c80a4b.svg} | 82 +++++------ .../exponential_integrators/index.html | 10 +- dev/tutorials/fenrir/17dc9379.svg | 126 ---------------- dev/tutorials/fenrir/3f2e3f62.svg | 126 ++++++++++++++++ dev/tutorials/fenrir/7974861f.svg | 131 +++++++++++++++++ dev/tutorials/fenrir/e5275bab.svg | 131 ----------------- dev/tutorials/fenrir/index.html | 14 +- .../{5e51c0bd.svg => 63d2c42c.svg} | 72 +++++----- .../{eeead470.svg => 97855bee.svg} | 80 +++++------ .../{08155c4a.svg => ddf784da.svg} | 72 +++++----- dev/tutorials/getting_started/index.html | 8 +- 36 files changed, 846 insertions(+), 846 deletions(-) rename dev/tutorials/dae/{81b0c82b.svg => 23ce2a2c.svg} (97%) rename dev/tutorials/dae/{12e69d49.svg => fb9c5960.svg} (99%) rename dev/tutorials/dynamical_odes/{e22b23c3.svg => 136579f1.svg} (98%) rename dev/tutorials/dynamical_odes/{5fcbcdb1.svg => 52e65c07.svg} (99%) rename dev/tutorials/dynamical_odes/{27923a4c.svg => ac97650f.svg} (88%) rename dev/tutorials/dynamical_odes/{373f8431.svg => e1d09d06.svg} (88%) rename dev/tutorials/exponential_integrators/{d844aea3.svg => 1c4f7b27.svg} (93%) rename dev/tutorials/exponential_integrators/{d7a376a7.svg => 9b15c58f.svg} (93%) rename dev/tutorials/exponential_integrators/{5e6b2b1f.svg => df738690.svg} (92%) rename dev/tutorials/exponential_integrators/{5f074e4e.svg => e9c80a4b.svg} (93%) delete mode 100644 dev/tutorials/fenrir/17dc9379.svg create mode 100644 dev/tutorials/fenrir/3f2e3f62.svg create mode 100644 dev/tutorials/fenrir/7974861f.svg delete mode 100644 dev/tutorials/fenrir/e5275bab.svg rename dev/tutorials/getting_started/{5e51c0bd.svg => 63d2c42c.svg} (96%) rename dev/tutorials/getting_started/{eeead470.svg => 97855bee.svg} (91%) rename dev/tutorials/getting_started/{08155c4a.svg => ddf784da.svg} (96%) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 550fc9309..a12f4c547 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-13T13:52:15","documenter_version":"1.1.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.3","generation_timestamp":"2023-10-16T08:07:35","documenter_version":"1.1.1"}} \ No newline at end of file diff --git a/dev/benchmarks/lotkavolterra/index.html b/dev/benchmarks/lotkavolterra/index.html index 0efc1a158..3748ea85e 100644 --- a/dev/benchmarks/lotkavolterra/index.html +++ b/dev/benchmarks/lotkavolterra/index.html @@ -597,4 +597,4 @@ [3f19e933] p7zip_jll v17.4.0+0 Info Packages marked with ⌃ and ⌅ have new versions available, but those wi th ⌅ are restricted by compatibility constraints from upgrading. To see why - use `status --outdated -m` + use `status --outdated -m` diff --git a/dev/benchmarks/multi-language-wrappers/index.html b/dev/benchmarks/multi-language-wrappers/index.html index 37b3402e7..cbda77104 100644 --- a/dev/benchmarks/multi-language-wrappers/index.html +++ b/dev/benchmarks/multi-language-wrappers/index.html @@ -688,4 +688,4 @@ [3f19e933] p7zip_jll v17.4.0+0 Info Packages marked with ⌃ and ⌅ have new versions available, but those wi th ⌅ are restricted by compatibility constraints from upgrading. To see why - use `status --outdated -m` + use `status --outdated -m` diff --git a/dev/benchmarks/pleiades/index.html b/dev/benchmarks/pleiades/index.html index f6ef311f7..8d1b4b663 100644 --- a/dev/benchmarks/pleiades/index.html +++ b/dev/benchmarks/pleiades/index.html @@ -525,4 +525,4 @@ [3f19e933] p7zip_jll v17.4.0+0 Info Packages marked with ⌃ and ⌅ have new versions available, but those wi th ⌅ are restricted by compatibility constraints from upgrading. To see why - use `status --outdated -m` + use `status --outdated -m` diff --git a/dev/benchmarks/rober/index.html b/dev/benchmarks/rober/index.html index dd37bcc64..ee2d85624 100644 --- a/dev/benchmarks/rober/index.html +++ b/dev/benchmarks/rober/index.html @@ -477,4 +477,4 @@ [3f19e933] p7zip_jll v17.4.0+0 Info Packages marked with ⌃ and ⌅ have new versions available, but those wi th ⌅ are restricted by compatibility constraints from upgrading. To see why - use `status --outdated -m` + use `status --outdated -m` diff --git a/dev/benchmarks/vanderpol/index.html b/dev/benchmarks/vanderpol/index.html index 15028e4fb..b055f763b 100644 --- a/dev/benchmarks/vanderpol/index.html +++ b/dev/benchmarks/vanderpol/index.html @@ -515,4 +515,4 @@ [3f19e933] p7zip_jll v17.4.0+0 Info Packages marked with ⌃ and ⌅ have new versions available, but those wi th ⌅ are restricted by compatibility constraints from upgrading. To see why - use `status --outdated -m` + use `status --outdated -m` diff --git a/dev/diffusions/index.html b/dev/diffusions/index.html index e25eea15b..d853002d2 100644 --- a/dev/diffusions/index.html +++ b/dev/diffusions/index.html @@ -5,9 +5,9 @@ \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1, \\ \text{d} Y^{(q)}(t) &= \textcolor{#389826}{A} Y(t) \ \text{d}t + \textcolor{#4063D8}{\Gamma} \ \text{d}W(t), \\ Y(0) &\sim \textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }. -\end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider choices relating to the "diffusion" $\textcolor{#4063D8}{\Gamma}$. If you're more interested in the drift matrix $\textcolor{#389826}{A}$ check out the Priors section, and for info on the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$ check out the Initialization section.

Diffusion and calibration

We call $\textcolor{#4063D8}{\Gamma}$ the "diffusion" parameter. Since it is typically not known we need to estimate it; this is called "calibration".

There are a few different choices for how to model and estimate $\textcolor{#4063D8}{\Gamma}$:

Or more compactly:

Isotropic:Diagonal (only for the EK0)
Time-varyingDynamicDiffusionDynamicMVDiffusion
Time-fixedFixedDiffusionFixedMVDiffusion

For more details on diffusions and calibration, check out this paper [6].

API

ProbNumDiffEq.DynamicDiffusion — Type
DynamicDiffusion()

Time-varying, isotropic diffusion, which is quasi-maximum-likelihood-estimated at each step.

This is the recommended diffusion when using adaptive step-size selection, and in particular also when solving stiff systems.

source
ProbNumDiffEq.FixedDiffusion — Type
FixedDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, isotropic diffusion, which is (optionally) quasi-maximum-likelihood-estimated.

This is the recommended diffusion when using fixed steps.

By default with calibrate=true, all covariances are re-scaled at the end of the solve with the MLE diffusion. Set calibrate=false to skip this step, e.g. when setting the initial_diffusion and then estimating the diffusion outside of the solver (e.g. with Fenrir.jl).

source
ProbNumDiffEq.DynamicMVDiffusion — Type
DynamicMVDiffusion()

Time-varying, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of DynamicDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [6] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source
ProbNumDiffEq.FixedMVDiffusion — Type
FixedMVDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of FixedDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [6] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source

References

[6]
+\end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider choices relating to the "diffusion" $\textcolor{#4063D8}{\Gamma}$. If you're more interested in the drift matrix $\textcolor{#389826}{A}$ check out the Priors section, and for info on the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$ check out the Initialization section.

Diffusion and calibration

We call $\textcolor{#4063D8}{\Gamma}$ the "diffusion" parameter. Since it is typically not known we need to estimate it; this is called "calibration".

There are a few different choices for how to model and estimate $\textcolor{#4063D8}{\Gamma}$:

  • FixedDiffusion assumes an isotropic, time-fixed $\textcolor{#4063D8}{\Gamma} = \sigma \cdot I_d$,
  • DynamicDiffusion assumes an isotropic, time-varying $\textcolor{#4063D8}{\Gamma}(t) = \sigma(t) \cdot I_d$ (recommended),
  • FixedMVDiffusion assumes a diagonal, time-fixed $\textcolor{#4063D8}{\Gamma} = \operatorname{diag}(\sigma_1, \dots, \sigma_d)$,
  • DynamicMVDiffusion assumes a diagonal, time-varying $\textcolor{#4063D8}{\Gamma}(t) = \operatorname{diag}(\sigma_1(t), \dots, \sigma_d(t))$.

Or more compactly:

Isotropic:Diagonal (only for the EK0)
Time-varyingDynamicDiffusionDynamicMVDiffusion
Time-fixedFixedDiffusionFixedMVDiffusion

For more details on diffusions and calibration, check out this paper [6].

API

ProbNumDiffEq.DynamicDiffusion — Type
DynamicDiffusion()

Time-varying, isotropic diffusion, which is quasi-maximum-likelihood-estimated at each step.

This is the recommended diffusion when using adaptive step-size selection, and in particular also when solving stiff systems.

source
ProbNumDiffEq.FixedDiffusion — Type
FixedDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, isotropic diffusion, which is (optionally) quasi-maximum-likelihood-estimated.

This is the recommended diffusion when using fixed steps.

By default with calibrate=true, all covariances are re-scaled at the end of the solve with the MLE diffusion. Set calibrate=false to skip this step, e.g. when setting the initial_diffusion and then estimating the diffusion outside of the solver (e.g. with Fenrir.jl).

source
ProbNumDiffEq.DynamicMVDiffusion — Type
DynamicMVDiffusion()

Time-varying, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of DynamicDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [6] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source
ProbNumDiffEq.FixedMVDiffusion — Type
FixedMVDiffusion(; initial_diffusion=1.0, calibrate=true)

Time-fixed, diagonal diffusion, which is quasi-maximum-likelihood-estimated at each step.

Only works with the EK0!

A multi-variate version of FixedDiffusion, where instead of an isotropic matrix, a diagonal matrix is estimated. This can be helpful to get more expressive posterior covariances when using the EK0, since the individual dimensions can be adjusted separately.

References

  • [6] Bosch et al, "Calibrated Adaptive Probabilistic ODE Solvers", AISTATS (2021)
source

References

[6]
N. Bosch, P. Hennig and F. Tronarp. Calibrated Adaptive Probabilistic ODE Solvers. In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, 3466–3474, PMLR (2021).
-
+ diff --git a/dev/filtering/index.html b/dev/filtering/index.html index d5632274b..f28d87a37 100644 --- a/dev/filtering/index.html +++ b/dev/filtering/index.html @@ -1,16 +1,16 @@ -Filtering and Smoothing · ProbNumDiffEq.jl
GitHub

Gaussian Filtering and Smoothing

Predict

ProbNumDiffEq.predict — Function
predict(x::Gaussian, A::AbstractMatrix, Q::AbstractMatrix)

Prediction step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$, compute and return $\mathcal{N}(A μ, A Σ A^T + Q)$.

See also the non-allocating square-root version predict!.

source
ProbNumDiffEq.predict! — Function
predict!(x_out, x_curr, Ah, Qh, cachemat)

In-place and square-root implementation of predict which saves the result into x_out.

Only works with PSDMatrices.PSDMatrix types as Ah, Qh, and in the covariances of x_curr and x_out (both of type Gaussian). To prevent allocations, a cache matrix cachemat of size $D \times 2D$ (where $D \times D$ is the size of Ah and Qh) needs to be passed.

See also: predict.

source

Update

ProbNumDiffEq.update — Function
update(x, measurement, H)

Update step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$ and a measurement $z = \mathcal{N}(\hat{z}, S)$, with $S = H Σ H^T$, compute

\[\begin{aligned} +Filtering and Smoothing · ProbNumDiffEq.jl

GitHub

Gaussian Filtering and Smoothing

Predict

ProbNumDiffEq.predict — Function
predict(x::Gaussian, A::AbstractMatrix, Q::AbstractMatrix)

Prediction step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$, compute and return $\mathcal{N}(A μ, A Σ A^T + Q)$.

See also the non-allocating square-root version predict!.

source
ProbNumDiffEq.predict! — Function
predict!(x_out, x_curr, Ah, Qh, cachemat)

In-place and square-root implementation of predict which saves the result into x_out.

Only works with PSDMatrices.PSDMatrix types as Ah, Qh, and in the covariances of x_curr and x_out (both of type Gaussian). To prevent allocations, a cache matrix cachemat of size $D \times 2D$ (where $D \times D$ is the size of Ah and Qh) needs to be passed.

See also: predict.

source

Update

ProbNumDiffEq.update — Function
update(x, measurement, H)

Update step in Kalman filtering for linear dynamics models.

Given a Gaussian $x = \mathcal{N}(μ, Σ)$ and a measurement $z = \mathcal{N}(\hat{z}, S)$, with $S = H Σ H^T$, compute

\[\begin{aligned} K &= Σ^P H^T S^{-1}, \\ μ^F &= μ + K (0 - \hat{z}), \\ Σ^F &= Σ - K S K^T, -\end{aligned}\]

and return an updated state \mathcal{N}(μ^F, Σ^F). Note that this assumes zero-measurements. When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating update!.

source
ProbNumDiffEq.update! — Function
update!(x_out, x_pred, measurement, H, K_cache, M_cache, S_cache)

In-place and square-root implementation of update which saves the result into x_out.

Implemented in Joseph Form to retain the PSDMatrix covariances:

\[\begin{aligned} +\end{aligned}\]

and return an updated state \mathcal{N}(μ^F, Σ^F). Note that this assumes zero-measurements. When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating update!.

source
ProbNumDiffEq.update! — Function
update!(x_out, x_pred, measurement, H, K_cache, M_cache, S_cache)

In-place and square-root implementation of update which saves the result into x_out.

Implemented in Joseph Form to retain the PSDMatrix covariances:

\[\begin{aligned} K &= Σ^P H^T S^{-1}, \\ μ^F &= μ + K (0 - \hat{z}), \\ \sqrt{Σ}^F &= (I - KH) \sqrt(Σ), -\end{aligned}\]

where $\sqrt{M}$ denotes the left square-root of a matrix M, i.e. $M = \sqrt{M} \sqrt{M}^T$.

To prevent allocations, write into caches K_cache and M_cache, both of size D × D, and S_cache of same type as measurement.Σ.

See also: update.

source

Smooth

ProbNumDiffEq.smooth — Function
smooth(x_curr, x_next_smoothed, A, Q)

Update step of the Kalman smoother, aka. Rauch-Tung-Striebel smoother, for linear dynamics models.

Given Gaussians $x_n = \mathcal{N}(μ_{n}, Σ_{n})$ and $x_{n+1} = \mathcal{N}(μ_{n+1}^S, Σ_{n+1}^S)$, compute

\[\begin{aligned} +\end{aligned}\]

where $\sqrt{M}$ denotes the left square-root of a matrix M, i.e. $M = \sqrt{M} \sqrt{M}^T$.

To prevent allocations, write into caches K_cache and M_cache, both of size D × D, and S_cache of same type as measurement.Σ.

See also: update.

source

Smooth

ProbNumDiffEq.smooth — Function
smooth(x_curr, x_next_smoothed, A, Q)

Update step of the Kalman smoother, aka. Rauch-Tung-Striebel smoother, for linear dynamics models.

Given Gaussians $x_n = \mathcal{N}(μ_{n}, Σ_{n})$ and $x_{n+1} = \mathcal{N}(μ_{n+1}^S, Σ_{n+1}^S)$, compute

\[\begin{aligned} μ_{n+1}^P &= A μ_n^F, \\ P_{n+1}^P &= A Σ_n^F A + Q, \\ G &= Σ_n^S A^T (Σ_{n+1}^P)^{-1}, \\ μ_n^S &= μ_n^F + G (μ_{n+1}^S - μ_{n+1}^P), \\ Σ_n^S &= (I - G A) Σ_n^F (I - G A)^T + G Q G^T + G Σ_{n+1}^S G^T, -\end{aligned}\]

and return a smoothed state \mathcal{N}(μ_n^S, Σ_n^S). When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating smooth!.

source
ProbNumDiffEq.smooth! — Function
smooth!(x_curr, x_next, Ah, Qh, cache, diffusion=1)

In-place and square-root implementation of smooth which overwrites x_curr.

Implemented in Joseph form to preserve square-root structure. It requires access to the solvers cache to prevent allocations.

See also: smooth.

source
+\end{aligned}\]

and return a smoothed state \mathcal{N}(μ_n^S, Σ_n^S). When called with ProbNumDiffEq.SquarerootMatrix type arguments it performs the update in Joseph / square-root form.

For better performance, we recommend to use the non-allocating smooth!.

source
ProbNumDiffEq.smooth! — Function
smooth!(x_curr, x_next, Ah, Qh, cache, diffusion=1)

In-place and square-root implementation of smooth which overwrites x_curr.

Implemented in Joseph form to preserve square-root structure. It requires access to the solvers cache to prevent allocations.

See also: smooth.

source
diff --git a/dev/implementation/index.html b/dev/implementation/index.html index 3df205dcf..df4f32de6 100644 --- a/dev/implementation/index.html +++ b/dev/implementation/index.html @@ -1,2 +1,2 @@ -Implementation via OrdinaryDiffEq.jl · ProbNumDiffEq.jl
GitHub

Solver Implementation via OrdinaryDiffEq.jl

ProbNumDiffEq.jl builds directly on OrdinaryDiffEq.jl to benefit from its iterator interface, flexible step-size control, and efficient Jacobian calculations. But, this requires extending non-public APIs. This page is meant to provide an overview on which parts exactly ProbNumDiffEq.jl builds on.

For more discussion on the pros and cons of building on OrdinaryDiffEq.jl, see this thread on discourse.

Building on OrdinaryDiffEq.jl

ProbNumDiffEq.jl shares most of OrdinaryDiffEq.jl's implementation. In particular:

  1. OrdinaryDiffEq.__init builds the cache and the integrator, and calls OrdinaryDiffEq.initialize!
  2. OrdinaryDiffEq.solve! implements the actual iterator structure, with
    • OrdinaryDiffEq.loopheader!
    • OrdinaryDiffEq.perform_step!
    • OrdinaryDiffEq.loopfooter!
    • OrdinaryDiffEq.postamble!

ProbNumDiffEq.jl builds around this structure and overloads some of the parts:

  • Algorithms: EK0/EK1 <: AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm
    • ./src/algorithms.jl provides the algorithms themselves
    • ./src/alg_utils.jl implements many traits (e.g. relating to autodiff, implicitness, step-size control)
  • Cache: EKCache <: AbstractODEFilterCache <: OrdinaryDiffEq.OrdinaryDiffEqCache
    • ./src/caches.jl implements the cache and its main constructor: OrdinaryDiffEq.alg_cache
  • Initialization and perform_step!: via OrdinaryDiffEq.initialize! and OrdinaryDiffEq.perform_step!. Implemented in ./src/perform_step.jl.
  • Custom postamble by overloading OrdinaryDiffEq.postamble! (which should always call OrdinaryDiffEq._postamble!). This is where we do the "smoothing" of the solution. Implemented in ./src/integrator_utils.jl.
  • Custom saving by overloading OrdinaryDiffEq.savevalues! (which should always call OrdinaryDiffEq._savevalues!). Implemented in ./src/integrator_utils.jl.

Building on DiffEqBase.jl

  • DiffEqBase.__init is currently overloaded to transform OOP problems into IIP problems (in ./src/solve.jl).
  • The solution object: ProbODESolution <: AbstractProbODESolution <: DiffEqBase.AbstractODESolution
    • ./src/solution.jl implements the main parts. Note that the main constructor DiffEqBase.build_solution is called by OrdinaryDiffEq.__init - so OrdinaryDiffEq.jl has control over its inputs.
    • There is also MeanProbODESolution <: DiffEqBase.AbstractODESolution: It allows handling the mean of a probabilistic ODE solution the same way one would handle any "standard" ODE solution - e.g. it is compatible with DiffEqDevTools.appxtrue.
    • AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation is the current interpolant, but it does not actually fully handle the interpolation right now. This part might be subject to change soon.
    • Plot recipe in ./src/solution_plotting.jl
    • Sampling in ./src/solution_sampling.jl
  • DiffEqBase.prepare_alg(::EK1{0}); closely follows a similar function implemented in OrdinaryDiffEq.jl ./src/alg_utils.jl
    • this also required DiffEqBase.remake(::EK1)

Other packages

  • DiffEqDevTools.appxtrue is overloaded to work with ProbODESolution (by just doing mean(sol)). This also enables DiffEqDevTools.WorkPrecision to work out of th box.
+Implementation via OrdinaryDiffEq.jl · ProbNumDiffEq.jl
GitHub

Solver Implementation via OrdinaryDiffEq.jl

ProbNumDiffEq.jl builds directly on OrdinaryDiffEq.jl to benefit from its iterator interface, flexible step-size control, and efficient Jacobian calculations. But, this requires extending non-public APIs. This page is meant to provide an overview on which parts exactly ProbNumDiffEq.jl builds on.

For more discussion on the pros and cons of building on OrdinaryDiffEq.jl, see this thread on discourse.

Building on OrdinaryDiffEq.jl

ProbNumDiffEq.jl shares most of OrdinaryDiffEq.jl's implementation. In particular:

  1. OrdinaryDiffEq.__init builds the cache and the integrator, and calls OrdinaryDiffEq.initialize!
  2. OrdinaryDiffEq.solve! implements the actual iterator structure, with
    • OrdinaryDiffEq.loopheader!
    • OrdinaryDiffEq.perform_step!
    • OrdinaryDiffEq.loopfooter!
    • OrdinaryDiffEq.postamble!

ProbNumDiffEq.jl builds around this structure and overloads some of the parts:

  • Algorithms: EK0/EK1 <: AbstractEK <: OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm
    • ./src/algorithms.jl provides the algorithms themselves
    • ./src/alg_utils.jl implements many traits (e.g. relating to autodiff, implicitness, step-size control)
  • Cache: EKCache <: AbstractODEFilterCache <: OrdinaryDiffEq.OrdinaryDiffEqCache
    • ./src/caches.jl implements the cache and its main constructor: OrdinaryDiffEq.alg_cache
  • Initialization and perform_step!: via OrdinaryDiffEq.initialize! and OrdinaryDiffEq.perform_step!. Implemented in ./src/perform_step.jl.
  • Custom postamble by overloading OrdinaryDiffEq.postamble! (which should always call OrdinaryDiffEq._postamble!). This is where we do the "smoothing" of the solution. Implemented in ./src/integrator_utils.jl.
  • Custom saving by overloading OrdinaryDiffEq.savevalues! (which should always call OrdinaryDiffEq._savevalues!). Implemented in ./src/integrator_utils.jl.

Building on DiffEqBase.jl

  • DiffEqBase.__init is currently overloaded to transform OOP problems into IIP problems (in ./src/solve.jl).
  • The solution object: ProbODESolution <: AbstractProbODESolution <: DiffEqBase.AbstractODESolution
    • ./src/solution.jl implements the main parts. Note that the main constructor DiffEqBase.build_solution is called by OrdinaryDiffEq.__init - so OrdinaryDiffEq.jl has control over its inputs.
    • There is also MeanProbODESolution <: DiffEqBase.AbstractODESolution: It allows handling the mean of a probabilistic ODE solution the same way one would handle any "standard" ODE solution - e.g. it is compatible with DiffEqDevTools.appxtrue.
    • AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation is the current interpolant, but it does not actually fully handle the interpolation right now. This part might be subject to change soon.
    • Plot recipe in ./src/solution_plotting.jl
    • Sampling in ./src/solution_sampling.jl
  • DiffEqBase.prepare_alg(::EK1{0}); closely follows a similar function implemented in OrdinaryDiffEq.jl ./src/alg_utils.jl
    • this also required DiffEqBase.remake(::EK1)

Other packages

  • DiffEqDevTools.appxtrue is overloaded to work with ProbODESolution (by just doing mean(sol)). This also enables DiffEqDevTools.WorkPrecision to work out of th box.
diff --git a/dev/index.html b/dev/index.html index de4fa4a9d..3629afcec 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Home · ProbNumDiffEq.jl
GitHub

Probabilistic Numerical Differential Equation Solvers

Banner

ProbNumDiffEq.jl provides probabilistic numerical solvers to the DifferentialEquations.jl ecosystem. The implemented ODE filters solve differential equations via Bayesian filtering and smoothing and compute not just a single point estimate of the true solution, but a posterior distribution that contains an estimate of its numerical approximation error.

For a short intro video, check out our poster presentation at JuliaCon2021.

Installation

Run Julia, enter ] to bring up Julia's package manager, and add the ProbNumDiffEq.jl package:

julia> ]
-(v1.9) pkg> add ProbNumDiffEq

Getting Started

For a quick introduction check out the "Solving ODEs with Probabilistic Numerics" tutorial.

Features

  • probdiffeq: Fast and feature-rich filtering-based probabilistic ODE solvers in JAX.
  • ProbNum: Probabilistic numerics in Python. It has not only probabilistic ODE solvers, but also probabilistic linear solvers, Bayesian quadrature, and many filtering and smoothing implementations.
  • Fenrir.jl: Parameter-inference in ODEs with probabilistic ODE solvers. This package builds on ProbNumDiffEq.jl to provide a negative marginal log-likelihood function, which can then be used with an optimizer or with MCMC for parameter inference.
+(v1.9) pkg> add ProbNumDiffEq

Getting Started

For a quick introduction check out the "Solving ODEs with Probabilistic Numerics" tutorial.

Features

diff --git a/dev/initialization/index.html b/dev/initialization/index.html index bb2549674..648872b2e 100644 --- a/dev/initialization/index.html +++ b/dev/initialization/index.html @@ -11,7 +11,7 @@ \end{aligned}\]

It is clear that this contains quite some information for $Y(0)$: The initial value $y_0$ and the vector field $f$ imply

\[\begin{aligned} Y^{(0)}(0) &= y_0, \\ Y^{(1)}(0) &= f(y_0, 0). -\end{aligned}\]

It turns out that we can also compute higher-order derivatives of $y$ with the chain rule, and then use these to better initialize $Y^{(i)}(0)$. This, done efficiently with Taylor-mode autodiff by using TaylorIntegration.jl, is what TaylorModeInit does. See also [1].

In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!

API

ProbNumDiffEq.TaylorModeInit — Type
TaylorModeInit()

Exact initialization via Taylor-mode automatic differentiation.

This is the recommended initialization method!

It uses TaylorIntegration.jl to efficiently compute the higher-order derivatives of the solution at the initial value, 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.

References

  • [4] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
source
ProbNumDiffEq.ClassicSolverInit — Type
ClassicSolverInit(; alg=OrdinaryDiffEq.Tsit5(), init_on_ddu=false)

Initialization via regression on a few steps of a classic ODE solver.

In a nutshell, instead of specifying $\mu_0$ exactly and setting $\Sigma_0=0$ (which is what TaylorModeInit does), use a classic ODE solver to compute a few steps of the solution, and then regress on the computed values (by running a smoother) to compute $\mu_0$ and $\Sigma_0$ as the mean and covariance of the smoothing posterior at time 0. See also [2].

The initial value and derivative are set directly from the given initial value problem; optionally the second derivative can also be set via automatic differentiation by setting init_on_ddu=true.

Arguments

  • alg: The solver to be used. Can be any solver from OrdinaryDiffEq.jl.
  • init_on_ddu: If true, the second derivative is also initialized exactly via automatic differentiation with ForwardDiff.jl.

References

  • [4] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
  • [5] Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
source

References

[4]
+\end{aligned}\]

It turns out that we can also compute higher-order derivatives of $y$ with the chain rule, and then use these to better initialize $Y^{(i)}(0)$. This, done efficiently with Taylor-mode autodiff by using TaylorIntegration.jl, is what TaylorModeInit does. See also [1].

In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!

API

ProbNumDiffEq.TaylorModeInit — Type
TaylorModeInit()

Exact initialization via Taylor-mode automatic differentiation.

This is the recommended initialization method!

It uses TaylorIntegration.jl to efficiently compute the higher-order derivatives of the solution at the initial value, 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.

References

  • [4] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
source
ProbNumDiffEq.ClassicSolverInit — Type
ClassicSolverInit(; alg=OrdinaryDiffEq.Tsit5(), init_on_ddu=false)

Initialization via regression on a few steps of a classic ODE solver.

In a nutshell, instead of specifying $\mu_0$ exactly and setting $\Sigma_0=0$ (which is what TaylorModeInit does), use a classic ODE solver to compute a few steps of the solution, and then regress on the computed values (by running a smoother) to compute $\mu_0$ and $\Sigma_0$ as the mean and covariance of the smoothing posterior at time 0. See also [2].

The initial value and derivative are set directly from the given initial value problem; optionally the second derivative can also be set via automatic differentiation by setting init_on_ddu=true.

Arguments

  • alg: The solver to be used. Can be any solver from OrdinaryDiffEq.jl.
  • init_on_ddu: If true, the second derivative is also initialized exactly via automatic differentiation with ForwardDiff.jl.

References

  • [4] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
  • [5] Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)
source

References

[4]
[5]
@@ -19,4 +19,4 @@ -
+ diff --git a/dev/priors/index.html b/dev/priors/index.html index 1d2dfa516..3d268061f 100644 --- a/dev/priors/index.html +++ b/dev/priors/index.html @@ -8,12 +8,12 @@ \end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider choices relating to the drift matrix $\textcolor{#389826}{A}$. If you're more interested in the diffusion $\textcolor{#4063D8}{\Gamma}$ check out the Diffusion models and calibration section, and for info on the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$ check out the Initialization section.

If you're unsure which prior to use, just stick to the integrated Wiener process prior IWP! This is also the default choice for all solvers. The other priors are rather experimental / niche at the time of writing.

API

ProbNumDiffEq.IWP — Type
IWP([wiener_process_dimension::Integer,] num_derivatives::Integer)

Integrated Wiener process.

This is the recommended prior! It is the most well-tested prior, both in this package and in the probabilistic numerics literature in general (see the references). It is also the prior that has the most efficient implementation.

The IWP can be created without specifying the dimension of the Wiener process, in which case it will be inferred from the dimension of the ODE during the solve. This is typically the preferred usage.

In math

\[\begin{aligned} \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\ \text{d} Y^{(q)}(t) &= \Gamma \ \text{d}W(t). -\end{aligned}\]

Examples

julia> solve(prob, EK1(prior=IWP(2)))
source
ProbNumDiffEq.IOUP — Type
IOUP([wiener_process_dimension::Integer,]
+\end{aligned}\]
Examples
julia> solve(prob, EK1(prior=IWP(2)))
source
ProbNumDiffEq.IOUP — Type
IOUP([wiener_process_dimension::Integer,]
      num_derivatives::Integer,
      rate_parameter::Union{Number,AbstractVector,AbstractMatrix})

Integrated Ornstein–Uhlenbeck process.

As with the IWP, the IOUP can be created without specifying its dimension, in which case it will be inferred from the dimension of the ODE during the solve. This is typically the preferred usage. The rate parameter however always needs to be specified.

In math

\[\begin{aligned} \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\ \text{d} Y^{(q)}(t) &= L Y^{(q)}(t) \ \text{d}t + \Gamma \ \text{d}W(t), -\end{aligned}\]

where $L$ is the rate_parameter.

Examples

julia> solve(prob, EK1(prior=IOUP(2, -1)))
source
ProbNumDiffEq.Matern — Type
Matern([wiener_process_dimension::Integer,]
+\end{aligned}\]
where $L$ is the rate_parameter.
Examples
julia> solve(prob, EK1(prior=IOUP(2, -1)))
source
ProbNumDiffEq.Matern — Type
Matern([wiener_process_dimension::Integer,]
        num_derivatives::Integer,
        lengthscale::Number)

Matern process.

As with the IWP, the Matern can be created without specifying its dimension, in which case it will be inferred from the dimension of the ODE during the solve. This is typically the preferred usage. The lengthscale parameter however always needs to be specified.

In math

\[\begin{aligned} \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\ @@ -21,4 +21,4 @@ \begin{pmatrix} q+1 \\ j \end{pmatrix} \left( \frac{\sqrt{2q - 1}}{l} \right)^{q-j} Y^{(j)}(t) \right) \ \text{d}t + \Gamma \ \text{d}W(t). -\end{aligned}\]

where $l$ is the lengthscale.

Examples

julia> solve(prob, EK1(prior=Matern(2, 1)))
source
+\end{aligned}\]

where $l$ is the lengthscale.

Examples

julia> solve(prob, EK1(prior=Matern(2, 1)))
source diff --git a/dev/references/index.html b/dev/references/index.html index 05618e37e..0f617627e 100644 --- a/dev/references/index.html +++ b/dev/references/index.html @@ -45,4 +45,4 @@
- + diff --git a/dev/solvers/index.html b/dev/solvers/index.html index be48bde62..13af02986 100644 --- a/dev/solvers/index.html +++ b/dev/solvers/index.html @@ -4,14 +4,14 @@ prior=IWP(order), diffusionmodel=DynamicDiffusion(), initialization=TaylorModeInit(), - kwargs...)

Gaussian ODE filter with first-order vector field linearization.

Arguments

Some additional kwargs relating to implicit solvers are supported; check out DifferentialEquations.jl's Extra Options page. Right now, we support autodiff, chunk_size, and diff_type. In particular, autodiff=false can come in handy to use finite differences instead of ForwardDiff.jl to compute Jacobians.

Examples

julia> solve(prob, EK1())

References

source
ProbNumDiffEq.EK0 — Type
EK0(; order=3,
+      kwargs...)

Gaussian ODE filter with first-order vector field linearization.

Arguments

  • order::Integer: Order of the integrated Wiener process (IWP) prior.
  • smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
  • prior::AbstractODEFilterPrior: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
  • diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
  • initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Some additional kwargs relating to implicit solvers are supported; check out DifferentialEquations.jl's Extra Options page. Right now, we support autodiff, chunk_size, and diff_type. In particular, autodiff=false can come in handy to use finite differences instead of ForwardDiff.jl to compute Jacobians.

Examples

julia> solve(prob, EK1())

References

source
ProbNumDiffEq.EK0 — Type
EK0(; order=3,
       smooth=true,
       prior=IWP(order),
       diffusionmodel=DynamicDiffusion(),
-      initialization=TaylorModeInit())

Gaussian ODE filter with zeroth-order vector field linearization.

Arguments

  • order::Integer: Order of the integrated Wiener process (IWP) prior.
  • smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
  • prior::AbstractODEFilterPrior: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
  • diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
  • initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Examples

julia> solve(prob, EK0())

References

source

Probabilistic Exponential Integrators

ProbNumDiffEq.ExpEK — Function
ExpEK(; L, order=3, kwargs...)

Probabilistic exponential integrator

Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs that provide improved stability by essentially solving the linear part of the ODE exactly. In probabilistic numerics, this amounts to including the linear part into the prior model of the solver.

ExpEK is therefore just a short-hand for EK0 with IOUP prior:

ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)

See also RosenbrockExpEK, EK0, EK1.

Arguments

See EK0 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
-julia> solve(prob, ExpEK(L=-1))

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source
ProbNumDiffEq.RosenbrockExpEK — Function
RosenbrockExpEK(; order=3, kwargs...)

Probabilistic Rosenbrock-type exponential integrator

A probabilistic exponential integrator similar to ExpEK, but with automatic linearization along the mean numerical solution. This brings the advantage that the linearity does not need to be specified manually, and the more accurate local linearization can sometimes also improve stability; but since the "prior" is adjusted at each step the probabilistic interpretation becomes more complicated.

RosenbrockExpEK is just a short-hand for EK1 with appropriete IOUP prior:

RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)

See also ExpEK, EK0, EK1.

Arguments

See EK1 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
-julia> solve(prob, RosenbrockExpEK())

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source

References

[2]
+ initialization=TaylorModeInit())

Gaussian ODE filter with zeroth-order vector field linearization.

Arguments

  • order::Integer: Order of the integrated Wiener process (IWP) prior.
  • smooth::Bool: Turn smoothing on/off; smoothing is required for dense output.
  • prior::AbstractODEFilterPrior: Prior to be used by the ODE filter. By default, uses a 3-times integrated Wiener process prior IWP(3). See also: Priors.
  • diffusionmodel::ProbNumDiffEq.AbstractDiffusion: See Diffusion models and calibration.
  • initialization::ProbNumDiffEq.InitializationScheme: See Initialization.

Examples

julia> solve(prob, EK0())

References

source

Probabilistic Exponential Integrators

ProbNumDiffEq.ExpEK — Function
ExpEK(; L, order=3, kwargs...)

Probabilistic exponential integrator

Probabilistic exponential integrators are a class of integrators for semi-linear stiff ODEs that provide improved stability by essentially solving the linear part of the ODE exactly. In probabilistic numerics, this amounts to including the linear part into the prior model of the solver.

ExpEK is therefore just a short-hand for EK0 with IOUP prior:

ExpEK(; order=3, L, kwargs...) = EK0(; prior=IOUP(order, L), kwargs...)

See also RosenbrockExpEK, EK0, EK1.

Arguments

See EK0 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, ExpEK(L=-1))

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source
ProbNumDiffEq.RosenbrockExpEK — Function
RosenbrockExpEK(; order=3, kwargs...)

Probabilistic Rosenbrock-type exponential integrator

A probabilistic exponential integrator similar to ExpEK, but with automatic linearization along the mean numerical solution. This brings the advantage that the linearity does not need to be specified manually, and the more accurate local linearization can sometimes also improve stability; but since the "prior" is adjusted at each step the probabilistic interpretation becomes more complicated.

RosenbrockExpEK is just a short-hand for EK1 with appropriete IOUP prior:

RosenbrockExpEK(; order=3, kwargs...) = EK1(; prior=IOUP(order, update_rate_parameter=true), kwargs...)

See also ExpEK, EK0, EK1.

Arguments

See EK1 for available keyword arguments.

Examples

julia> prob = ODEProblem((du, u, p, t) -> (@. du = - u + sin(u)), [1.0], (0.0, 10.0))
+julia> solve(prob, RosenbrockExpEK())

Reference

  • [2] Bosch et al, "Probabilistic Exponential Integrators", arXiv (2021)
source

References

[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators (2023), arXiv:2305.14978 [math.NA].
-
+ diff --git a/dev/tutorials/dae/81b0c82b.svg b/dev/tutorials/dae/23ce2a2c.svg similarity index 97% rename from dev/tutorials/dae/81b0c82b.svg rename to dev/tutorials/dae/23ce2a2c.svg index 9bd03ac67..7b3ac3161 100644 --- a/dev/tutorials/dae/81b0c82b.svg +++ b/dev/tutorials/dae/23ce2a2c.svg @@ -1,59 +1,59 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dae/12e69d49.svg b/dev/tutorials/dae/fb9c5960.svg similarity index 99% rename from dev/tutorials/dae/12e69d49.svg rename to dev/tutorials/dae/fb9c5960.svg index eb2db0a8b..12910f676 100644 --- a/dev/tutorials/dae/12e69d49.svg +++ b/dev/tutorials/dae/fb9c5960.svg @@ -1,90 +1,90 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dae/index.html b/dev/tutorials/dae/index.html index cb53a7f11..83a113c85 100644 --- a/dev/tutorials/dae/index.html +++ b/dev/tutorials/dae/index.html @@ -28,7 +28,7 @@ ylabel=["u₁(t)" "u₂(t)" "u₃(t)"], xlabel=["" "" "t"], denseplot=false, -)Example block output

Looks good!

Solving an Index-3 DAE directly

The following is based on the "Automatic Index Reduction of DAEs" tutorial by ModelingToolkit.jl, which demonstrates how the classic Rodas4 solver fails to solve a DAE due to the fact that it is of index 3; which is why ModelingToolkit's automatic index reduction is so useful.

It turns out that our probabilistic numerical solvers can directly solve the index-3 DAE!

First, define the pendulum problem as in the tutorial:

function pendulum!(du, u, p, t)
+)
Example block output

Looks good!

Solving an Index-3 DAE directly

The following is based on the "Automatic Index Reduction of DAEs" tutorial by ModelingToolkit.jl, which demonstrates how the classic Rodas4 solver fails to solve a DAE due to the fact that it is of index 3; which is why ModelingToolkit's automatic index reduction is so useful.

It turns out that our probabilistic numerical solvers can directly solve the index-3 DAE!

First, define the pendulum problem as in the tutorial:

function pendulum!(du, u, p, t)
     x, dx, y, dy, T = u
     g, L = p
     du[1] = dx
@@ -104,7 +104,7 @@
  [0.9606564704752983, -0.6488912965111842, -0.2777393773684324, -2.241904545274168, -8.20649809026962]
  [0.953580703532598, -0.7329362795843499, -0.301137606051972, -2.317049716858693, -8.911119223489587]
  [0.9472620899569895, -0.8049576742312687, -0.32045990833782334, -2.3741938112672893, -9.49621461341327]
- [0.9445849315579167, -0.8347208413797078, -0.32826712908215605, -2.3960584356227934, -9.733564093090328]

Nope! The EK1 is able to solve the index-3 DAE directly. Pretty cool!

plot(sol)
Example block output

Is index-reduction still worth it?

The point of the "Automatic Index Reduction of DAEs" tutorial is to demonstrate ModelingToolkit's utility for automatic index reduction, which enables the classic implicit Runge-Kutta solvers such as Rodas5 to solve this DAE. Let's see if that still helps in this context here.

First, modelingtoolkitize the problem:

traced_sys = modelingtoolkitize(pendulum_prob)

\[ \begin{align} + [0.9445849315579167, -0.8347208413797078, -0.32826712908215605, -2.3960584356227934, -9.733564093090328]

Nope! The EK1 is able to solve the index-3 DAE directly. Pretty cool!

plot(sol)
Example block output

Is index-reduction still worth it?

The point of the "Automatic Index Reduction of DAEs" tutorial is to demonstrate ModelingToolkit's utility for automatic index reduction, which enables the classic implicit Runge-Kutta solvers such as Rodas5 to solve this DAE. Let's see if that still helps in this context here.

First, modelingtoolkitize the problem:

traced_sys = modelingtoolkitize(pendulum_prob)

\[ \begin{align} \frac{\mathrm{d} x_1\left( t \right)}{\mathrm{d}t} =& x_2\left( t \right) \\ \frac{\mathrm{d} x_2\left( t \right)}{\mathrm{d}t} =& x_1\left( t \right) x_5\left( t \right) \\ \frac{\mathrm{d} x_3\left( t \right)}{\mathrm{d}t} =& x_4\left( t \right) \\ @@ -140,4 +140,4 @@ Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, 10015–10027, PMLR (2022). -

+ diff --git a/dev/tutorials/dynamical_odes/e22b23c3.svg b/dev/tutorials/dynamical_odes/136579f1.svg similarity index 98% rename from dev/tutorials/dynamical_odes/e22b23c3.svg rename to dev/tutorials/dynamical_odes/136579f1.svg index cbcabf373..3e29c1bdd 100644 --- a/dev/tutorials/dynamical_odes/e22b23c3.svg +++ b/dev/tutorials/dynamical_odes/136579f1.svg @@ -1,50 +1,50 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/5fcbcdb1.svg b/dev/tutorials/dynamical_odes/52e65c07.svg similarity index 99% rename from dev/tutorials/dynamical_odes/5fcbcdb1.svg rename to dev/tutorials/dynamical_odes/52e65c07.svg index 99338afc5..16b13beaf 100644 --- a/dev/tutorials/dynamical_odes/5fcbcdb1.svg +++ b/dev/tutorials/dynamical_odes/52e65c07.svg @@ -1,52 +1,52 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/27923a4c.svg b/dev/tutorials/dynamical_odes/ac97650f.svg similarity index 88% rename from dev/tutorials/dynamical_odes/27923a4c.svg rename to dev/tutorials/dynamical_odes/ac97650f.svg index 16bf07c0a..827da79d4 100644 --- a/dev/tutorials/dynamical_odes/27923a4c.svg +++ b/dev/tutorials/dynamical_odes/ac97650f.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/373f8431.svg b/dev/tutorials/dynamical_odes/e1d09d06.svg similarity index 88% rename from dev/tutorials/dynamical_odes/373f8431.svg rename to dev/tutorials/dynamical_odes/e1d09d06.svg index 6e17b7d38..a17a5f392 100644 --- a/dev/tutorials/dynamical_odes/373f8431.svg +++ b/dev/tutorials/dynamical_odes/e1d09d06.svg @@ -1,46 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/dynamical_odes/index.html b/dev/tutorials/dynamical_odes/index.html index 62c6079a8..31548a9aa 100644 --- a/dev/tutorials/dynamical_odes/index.html +++ b/dev/tutorials/dynamical_odes/index.html @@ -19,13 +19,13 @@ tspan = (0.0, 100.0) prob = ODEProblem(Hénon_Heiles, [du0; u0], tspan) sol = solve(prob, EK1()); -plot(sol, idxs=(3, 4)) # where `idxs=(3,4)` is used to plot x agains yExample block output

Solving the second-order ODE directly

Instead of first transforming the problem, we can also solve it directly as a second-order ODE, by defining it as a SecondOrderODEProblem.

Note

The SecondOrderODEProblem type is not defined in ProbNumDiffEq.jl but is provided by SciMLBase.jl. For more information, check out the DifferentialEquations.jl documentation on Dynamical, Hamiltonian and 2nd Order ODE Problems.

function Hénon_Heiles2(ddu, du, u, p, t)
+plot(sol, idxs=(3, 4)) # where `idxs=(3,4)` is used to plot x agains y
Example block output

Solving the second-order ODE directly

Instead of first transforming the problem, we can also solve it directly as a second-order ODE, by defining it as a SecondOrderODEProblem.

Note

The SecondOrderODEProblem type is not defined in ProbNumDiffEq.jl but is provided by SciMLBase.jl. For more information, check out the DifferentialEquations.jl documentation on Dynamical, Hamiltonian and 2nd Order ODE Problems.

function Hénon_Heiles2(ddu, du, u, p, t)
     ddu[1] = -u[1] - 2 * u[1] * u[2]
     ddu[2] = u[2]^2 - u[2] - u[1]^2
 end
 prob2 = SecondOrderODEProblem(Hénon_Heiles2, du0, u0, tspan)
 sol2 = solve(prob2, EK1());
-plot(sol2, idxs=(3, 4))
Example block output

Benchmark: Solving second order ODEs is faster

Solving second-order ODEs is not just a matter of convenience - in fact, SciMLBase's SecondOrderODEProblem is neatly designed in such a way that all the classic solvers from OrdinaryDiffEq.jl can handle it by solving the corresponding first-order ODE. But, transforming the ODE to first order increases the dimensionality of the problem, and comes therefore at increased computational cost; this also motivates classic specialized solvers for second-order ODEs.

The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the measurement model [1]. As a result, we can use the EK1 both for first and second order ODEs, but it automatically specializes on the latter to provide a 2x performance boost:

julia> @btime solve(prob, EK1(order=3), adaptive=false, dt=1e-2);
+plot(sol2, idxs=(3, 4))
Example block output

Benchmark: Solving second order ODEs is faster

Solving second-order ODEs is not just a matter of convenience - in fact, SciMLBase's SecondOrderODEProblem is neatly designed in such a way that all the classic solvers from OrdinaryDiffEq.jl can handle it by solving the corresponding first-order ODE. But, transforming the ODE to first order increases the dimensionality of the problem, and comes therefore at increased computational cost; this also motivates classic specialized solvers for second-order ODEs.

The probabilistic numerical solvers from ProbNumDiffEq.jl have the same internal state representation for first and second order ODEs; all that changes is the measurement model [1]. As a result, we can use the EK1 both for first and second order ODEs, but it automatically specializes on the latter to provide a 2x performance boost:

julia> @btime solve(prob, EK1(order=3), adaptive=false, dt=1e-2);
   766.312 ms (400362 allocations: 173.38 MiB)
 
 julia> @btime solve(prob2, EK1(order=4), adaptive=false, dt=1e-2);
@@ -37,13 +37,13 @@
 E(dx, dy, x, y) = PotentialEnergy(x, y) + KineticEnergy(dx, dy)
 E(u) = E(u...); # convenient shorthand
E (generic function with 2 methods)

So, let's have a look at how the total energy changes over time when we numerically simulate the Hénon-Heiles model over a long period of time: Standard solve

longprob = remake(prob2, tspan=(0.0, 1e3))
 longsol = solve(longprob, EK1(smooth=false), dense=false)
-plot(longsol.t, E.(longsol.u))
Example block output

It visibly loses energy over time, from an initial 0.12967 to a final 0.12899. Let's fix this to get a physically more meaningful solution.

Energy preservation with the ManifoldUpdate callback

In the language of ODE filters, preserving energy over time amounts to just another measurement model [1]. The most convenient way of updating on this additional zero measurement with ProbNumDiffEq.jl is with the ManifoldUpdate callback.

Note

The ManifoldUpdate callback can be thought of a probabilistic counterpart to the ManifoldProjection callback provided by DiffEqCallbacks.jl.

To do so, first define a (vector-valued) residual function, here chosen to be the difference between the current energy and the initial energy, and build a ManifoldUpdate callback

residual(u) = [E(u) - E(du0..., u0...)]
+plot(longsol.t, E.(longsol.u))
Example block output

It visibly loses energy over time, from an initial 0.12967 to a final 0.12899. Let's fix this to get a physically more meaningful solution.

Energy preservation with the ManifoldUpdate callback

In the language of ODE filters, preserving energy over time amounts to just another measurement model [1]. The most convenient way of updating on this additional zero measurement with ProbNumDiffEq.jl is with the ManifoldUpdate callback.

Note

The ManifoldUpdate callback can be thought of a probabilistic counterpart to the ManifoldProjection callback provided by DiffEqCallbacks.jl.

To do so, first define a (vector-valued) residual function, here chosen to be the difference between the current energy and the initial energy, and build a ManifoldUpdate callback

residual(u) = [E(u) - E(du0..., u0...)]
 cb = ManifoldUpdate(residual)
DiscreteCallback{ProbNumDiffEq.var"#condition#45", ProbNumDiffEq.var"#affect!#46"{Int64, Float64, Float64, typeof(Main.residual)}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(ProbNumDiffEq.var"#condition#45"(), ProbNumDiffEq.var"#affect!#46"{Int64, Float64, Float64, typeof(Main.residual)}(100, 1.0e-25, 1.0e-15, Main.residual), SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])

Then, solve the ODE with this callback

longsol_preserving = solve(longprob, EK1(smooth=false), dense=false, callback=cb)
 plot(longsol.t, E.(longsol.u))
-plot!(longsol_preserving.t, E.(longsol_preserving.u))
Example block output

Voilà! With the ManifoldUpdate callback we could preserve the energy over time and obtain a more truthful probabilistic numerical long-term simulation of the Hénon-Heiles model.

References

[1]
+plot!(longsol_preserving.t, E.(longsol_preserving.u))Example block output

Voilà! With the ManifoldUpdate callback we could preserve the energy over time and obtain a more truthful probabilistic numerical long-term simulation of the Hénon-Heiles model.

References

[1]
N. Bosch, F. Tronarp and P. Hennig. Pick-and-Mix Information Operators for Probabilistic ODE Solvers. In: Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, 10015–10027, PMLR (2022).
-
+ diff --git a/dev/tutorials/exponential_integrators/d844aea3.svg b/dev/tutorials/exponential_integrators/1c4f7b27.svg similarity index 93% rename from dev/tutorials/exponential_integrators/d844aea3.svg rename to dev/tutorials/exponential_integrators/1c4f7b27.svg index cf3336237..289160d11 100644 --- a/dev/tutorials/exponential_integrators/d844aea3.svg +++ b/dev/tutorials/exponential_integrators/1c4f7b27.svg @@ -1,51 +1,51 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/exponential_integrators/d7a376a7.svg b/dev/tutorials/exponential_integrators/9b15c58f.svg similarity index 93% rename from dev/tutorials/exponential_integrators/d7a376a7.svg rename to dev/tutorials/exponential_integrators/9b15c58f.svg index 36fb17dec..c52192d0a 100644 --- a/dev/tutorials/exponential_integrators/d7a376a7.svg +++ b/dev/tutorials/exponential_integrators/9b15c58f.svg @@ -1,63 +1,63 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/exponential_integrators/5e6b2b1f.svg b/dev/tutorials/exponential_integrators/df738690.svg similarity index 92% rename from dev/tutorials/exponential_integrators/5e6b2b1f.svg rename to dev/tutorials/exponential_integrators/df738690.svg index 198d619bb..96b2f335b 100644 --- a/dev/tutorials/exponential_integrators/5e6b2b1f.svg +++ b/dev/tutorials/exponential_integrators/df738690.svg @@ -1,65 +1,65 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/exponential_integrators/5f074e4e.svg b/dev/tutorials/exponential_integrators/e9c80a4b.svg similarity index 93% rename from dev/tutorials/exponential_integrators/5f074e4e.svg rename to dev/tutorials/exponential_integrators/e9c80a4b.svg index 421b015a2..bb25a6a9e 100644 --- a/dev/tutorials/exponential_integrators/5f074e4e.svg +++ b/dev/tutorials/exponential_integrators/e9c80a4b.svg @@ -1,55 +1,55 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/exponential_integrators/index.html b/dev/tutorials/exponential_integrators/index.html index 1f2bb259e..9f4c96254 100644 --- a/dev/tutorials/exponential_integrators/index.html +++ b/dev/tutorials/exponential_integrators/index.html @@ -12,7 +12,7 @@ prob = ODEProblem(f, u0, tspan) ref = solve(prob, EK1(), abstol=1e-10, reltol=1e-10) -plot(ref, color=:black, linestyle=:dash, label="Reference")Example block output

But for fixed (large) step sizes this ODE is more challenging: The explicit EK0 method oscillates and diverges due to the stiffness of the ODE, and the semi-implicit EK1 method is stable but the solution is not very accurate.

STEPSIZE = 4
+plot(ref, color=:black, linestyle=:dash, label="Reference")
Example block output

But for fixed (large) step sizes this ODE is more challenging: The explicit EK0 method oscillates and diverges due to the stiffness of the ODE, and the semi-implicit EK1 method is stable but the solution is not very accurate.

STEPSIZE = 4
 DM = FixedDiffusion() # recommended for fixed steps
 
 # we don't smooth the EK0 here to show the oscillations more clearly
@@ -22,18 +22,18 @@
 plot(ylims=(0.3, 1.05))
 plot!(ref, color=:black, linestyle=:dash, label="Reference")
 plot!(sol0, denseplot=false, marker=:o, markersize=2, label="EK0", color=1)
-plot!(sol1, denseplot=false, marker=:o, markersize=2, label="EK1", color=2)
Example block output

Probabilistic exponential integrators leverage the semi-linearity of the ODE to compute more accurate solutions for the same fixed step size. You can use either the ExpEK method and provide the linear part (with the keyword argument L), or the RosenbrockExpEK to automatically linearize along the mean of the numerical solution:

sol_exp = solve(prob, ExpEK(L=-1, diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+plot!(sol1, denseplot=false, marker=:o, markersize=2, label="EK1", color=2)
Example block output

Probabilistic exponential integrators leverage the semi-linearity of the ODE to compute more accurate solutions for the same fixed step size. You can use either the ExpEK method and provide the linear part (with the keyword argument L), or the RosenbrockExpEK to automatically linearize along the mean of the numerical solution:

sol_exp = solve(prob, ExpEK(L=-1, diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
 sol_ros = solve(prob, RosenbrockExpEK(diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
 
 plot(ylims=(0.3, 1.05))
 plot!(ref, color=:black, linestyle=:dash, label="Reference")
 plot!(sol_exp, denseplot=false, marker=:o, markersize=2, label="ExpEK", color=3)
-plot!(sol_ros, denseplot=false, marker=:o, markersize=2, label="RosenbrockExpEK", color=4)
Example block output

The solutions are indeed much more accurate than those of the standard EK1, for the same fixed step size!

Background: Integrated Ornstein-Uhlenbeck priors

Probabilistic exponential integrators "solve the linear part exactly" by including it into the prior model of the solver. Namely, the solver chooses a (q-times) integrated Ornstein-Uhlenbeck prior with rate parameter equal to the linearity. The ExpEK solver is just a short-hand for an EK0 with appropriate prior:

julia> ExpEK(order=3, L=-1) == EK0(prior=IOUP(3, -1))true

Similarly, the RosenbrockExpEK solver is also just a short-hand:

julia> RosenbrockExpEK(order=3) == EK1(prior=IOUP(3, update_rate_parameter=true))true

This means that you can also construct other probabilistic exponential integrators by hand! In this example the EK1 with IOUP prior with rate parameter -1 performs extremely well:

sol_expek1 = solve(prob, EK1(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
+plot!(sol_ros, denseplot=false, marker=:o, markersize=2, label="RosenbrockExpEK", color=4)
Example block output

The solutions are indeed much more accurate than those of the standard EK1, for the same fixed step size!

Background: Integrated Ornstein-Uhlenbeck priors

Probabilistic exponential integrators "solve the linear part exactly" by including it into the prior model of the solver. Namely, the solver chooses a (q-times) integrated Ornstein-Uhlenbeck prior with rate parameter equal to the linearity. The ExpEK solver is just a short-hand for an EK0 with appropriate prior:

julia> ExpEK(order=3, L=-1) == EK0(prior=IOUP(3, -1))true

Similarly, the RosenbrockExpEK solver is also just a short-hand:

julia> RosenbrockExpEK(order=3) == EK1(prior=IOUP(3, update_rate_parameter=true))true

This means that you can also construct other probabilistic exponential integrators by hand! In this example the EK1 with IOUP prior with rate parameter -1 performs extremely well:

sol_expek1 = solve(prob, EK1(prior=IOUP(3, -1), diffusionmodel=DM), adaptive=false, dt=STEPSIZE)
 
 plot(ylims=(0.3, 1.05))
 plot!(ref, color=:black, linestyle=:dash, label="Reference")
-plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")
Example block output

References

[2]
+plot!(sol_expek1, denseplot=false, marker=:o, markersize=2, label="EK1 + IOUP")Example block output

References

[2]
N. Bosch, P. Hennig and F. Tronarp. Probabilistic Exponential Integrators (2023), arXiv:2305.14978 [math.NA].
-
+ diff --git a/dev/tutorials/fenrir/17dc9379.svg b/dev/tutorials/fenrir/17dc9379.svg deleted file mode 100644 index c8e60b8cb..000000000 --- a/dev/tutorials/fenrir/17dc9379.svg +++ /dev/null @@ -1,126 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/fenrir/3f2e3f62.svg b/dev/tutorials/fenrir/3f2e3f62.svg new file mode 100644 index 000000000..ecc47a2a1 --- /dev/null +++ b/dev/tutorials/fenrir/3f2e3f62.svg @@ -0,0 +1,126 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/fenrir/7974861f.svg b/dev/tutorials/fenrir/7974861f.svg new file mode 100644 index 000000000..9961bb308 --- /dev/null +++ b/dev/tutorials/fenrir/7974861f.svg @@ -0,0 +1,131 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/fenrir/e5275bab.svg b/dev/tutorials/fenrir/e5275bab.svg deleted file mode 100644 index 2aba80ba2..000000000 --- a/dev/tutorials/fenrir/e5275bab.svg +++ /dev/null @@ -1,131 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/dev/tutorials/fenrir/index.html b/dev/tutorials/fenrir/index.html index 62f30656f..ee7e81b72 100644 --- a/dev/tutorials/fenrir/index.html +++ b/dev/tutorials/fenrir/index.html @@ -26,12 +26,12 @@ odedata = [true_sol(t) .+ sqrt(observation_noise_var) * randn(length(u0)) for t in times] plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""]) -scatter!(times, stack(odedata), markersize=2, markerstrokewidth=0.1, color=1, label=["Noisy Data" ""])Example block output

Our goal is then to recover the true parameter p (and thus also the true trajectoy plotted above) the noisy data.

Computing the negative log-likelihood

To do parameter inference - be it maximum-likelihod, maximum a posteriori, or full Bayesian inference with MCMC - we need to evaluate the likelihood of given a parameter estimate $\theta_\text{est}$. This is exactly what Fenrir.jl's fenrir_nll provides:

p_est = (0.1, 0.1, 2.0)
+scatter!(times, stack(odedata), markersize=2, markerstrokewidth=0.1, color=1, label=["Noisy Data" ""])
Example block output

Our goal is then to recover the true parameter p (and thus also the true trajectoy plotted above) the noisy data.

Computing the negative log-likelihood

To do parameter inference - be it maximum-likelihod, maximum a posteriori, or full Bayesian inference with MCMC - we need to evaluate the likelihood of given a parameter estimate $\theta_\text{est}$. This is exactly what Fenrir.jl's fenrir_nll provides:

p_est = (0.1, 0.1, 2.0)
 prob = remake(true_prob, p=p_est)
 data = (t=times, u=odedata)
 κ² = 1e10
 nll, _, _ = fenrir_nll(prob, data, observation_noise_var, κ²; dt=1e-1)
-nll
274.2589058328008

This is the negative marginal log-likelihood of the parameter p_est. You can use it as any other NLL: Optimize it to compute maximum-likelihood estimates or MAPs, or plug it into MCMC to sample from the posterior. In our paper [3] we compute MLEs by pairing Fenrir with Optimization.jl and ForwardDiff.jl. Let's quickly explore how to do this next.

Maximum-likelihood parameter inference

To compute a maximum-likelihood estimate (MLE), we just need to maximize $\theta \to p(\mathcal{D} \mid \theta)$ - that is, minimize the nll from above. We use Optimization.jl for this. First, define a loss function and create an OptimizationProblem

function loss(x, _)
+nll
253.4660556614911

This is the negative marginal log-likelihood of the parameter p_est. You can use it as any other NLL: Optimize it to compute maximum-likelihood estimates or MAPs, or plug it into MCMC to sample from the posterior. In our paper [3] we compute MLEs by pairing Fenrir with Optimization.jl and ForwardDiff.jl. Let's quickly explore how to do this next.

Maximum-likelihood parameter inference

To compute a maximum-likelihood estimate (MLE), we just need to maximize $\theta \to p(\mathcal{D} \mid \theta)$ - that is, minimize the nll from above. We use Optimization.jl for this. First, define a loss function and create an OptimizationProblem

function loss(x, _)
     ode_params = x[begin:end-1]
     prob = remake(true_prob, p=ode_params)
     κ² = exp(x[end]) # the diffusion parameter of the EK1
@@ -49,15 +49,15 @@
  2.0
  1.0

Then, just solve it! Here we use LBFGS:

optsol = solve(optprob, LBFGS())
 p_mle = optsol.u[1:3]
3-element Vector{Float64}:
- 0.18670836355256565
- 0.36196938701472103
- 2.8837585725081145

Success! The computed MLE is quite close to the true parameter which we used to generate the data. As a final step, let's plot the true solution, the data, and the result of the MLE:

plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""])
+ 0.19306758835549193
+ 0.2695184060347543
+ 2.9167155953642308

Success! The computed MLE is quite close to the true parameter which we used to generate the data. As a final step, let's plot the true solution, the data, and the result of the MLE:

plot(true_sol, color=:black, linestyle=:dash, label=["True Solution" ""])
 scatter!(times, stack(odedata), markersize=2, markerstrokewidth=0.1, color=1, label=["Noisy Data" ""])
 mle_sol = solve(remake(true_prob, p=p_mle), EK1())
-plot!(mle_sol, color=3, label=["MLE-parameter Solution" ""])
Example block output

Looks good!

Reference

[3]
+plot!(mle_sol, color=3, label=["MLE-parameter Solution" ""])Example block output

Looks good!

Reference

[3]
F. Tronarp, N. Bosch and P. Hennig. Fenrir: Physics-Enhanced Regression for Initial Value Problems. In: Proceedings of the 39th International Conference on Machine Learning, Proceedings of Machine Learning Research, 21776–21794, PMLR (2022).
-
+ diff --git a/dev/tutorials/getting_started/5e51c0bd.svg b/dev/tutorials/getting_started/63d2c42c.svg similarity index 96% rename from dev/tutorials/getting_started/5e51c0bd.svg rename to dev/tutorials/getting_started/63d2c42c.svg index 650f42d3e..132e8b8bf 100644 --- a/dev/tutorials/getting_started/5e51c0bd.svg +++ b/dev/tutorials/getting_started/63d2c42c.svg @@ -1,50 +1,50 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/eeead470.svg b/dev/tutorials/getting_started/97855bee.svg similarity index 91% rename from dev/tutorials/getting_started/eeead470.svg rename to dev/tutorials/getting_started/97855bee.svg index 56d0fb34e..bc4af05bb 100644 --- a/dev/tutorials/getting_started/eeead470.svg +++ b/dev/tutorials/getting_started/97855bee.svg @@ -1,54 +1,54 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/08155c4a.svg b/dev/tutorials/getting_started/ddf784da.svg similarity index 96% rename from dev/tutorials/getting_started/08155c4a.svg rename to dev/tutorials/getting_started/ddf784da.svg index 050015353..2d7818ef1 100644 --- a/dev/tutorials/getting_started/08155c4a.svg +++ b/dev/tutorials/getting_started/ddf784da.svg @@ -1,50 +1,50 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/tutorials/getting_started/index.html b/dev/tutorials/getting_started/index.html index f7413a569..edcae59d3 100644 --- a/dev/tutorials/getting_started/index.html +++ b/dev/tutorials/getting_started/index.html @@ -15,7 +15,7 @@ prob = ODEProblem(fitz, u0, tspan, p) sol = solve(prob, EK1()) -plot(sol)Example block output

Step 1: Define the problem

First, import ProbNumDiffEq.jl

using ProbNumDiffEq

Then, set up the ODEProblem exactly as you would in DifferentialEquations.jl. Define the vector field

function fitz(du, u, p, t)
+plot(sol)
Example block output

Step 1: Define the problem

First, import ProbNumDiffEq.jl

using ProbNumDiffEq

Then, set up the ODEProblem exactly as you would in DifferentialEquations.jl. Define the vector field

function fitz(du, u, p, t)
     a, b, c = p
     du[1] = c * (u[1] - u[1]^3 / 3 + u[2])
     du[2] = -(1 / c) * (u[1] - a - b * u[2])
@@ -66,13 +66,13 @@
  [2.035086153414764, 0.7169394922557271]
  [2.0210436519484176, 0.6717007174965183]
  [2.0104405118668915, 0.6383145073764349]

That's it! we just computed a probabilistic numerical ODE solution!

Step 3: Analyze the solution

Let's plot the result with Plots.jl.

using Plots
-plot(sol)
Example block output

Looks good! Looks like the EK1 managed to solve the Fitzhugh-Nagumo problem quite well.

Tip

To learn more about plotting ODE solutions, check out the plotting tutorial for DifferentialEquations.jl + Plots.jl provided here. Most of that works exactly as expected with ProbNumDiffEq.jl.

Plot the probabilistic error estimates

The plot above looks like a standard ODE solution – but it's not! The numerical errors are just so small that we can't see them in the plot, and the probabilistic error estimates are too. We can visualize them by plotting the errors and error estimates directly:

using OrdinaryDiffEq, Statistics
+plot(sol)
Example block output

Looks good! Looks like the EK1 managed to solve the Fitzhugh-Nagumo problem quite well.

Tip

To learn more about plotting ODE solutions, check out the plotting tutorial for DifferentialEquations.jl + Plots.jl provided here. Most of that works exactly as expected with ProbNumDiffEq.jl.

Plot the probabilistic error estimates

The plot above looks like a standard ODE solution – but it's not! The numerical errors are just so small that we can't see them in the plot, and the probabilistic error estimates are too. We can visualize them by plotting the errors and error estimates directly:

using OrdinaryDiffEq, Statistics
 reference = solve(prob, Vern9(), abstol=1e-9, reltol=1e-9, saveat=sol.t)
 errors = reduce(hcat, mean.(sol.pu) .- reference.u)'
 error_estimates = reduce(hcat, std.(sol.pu))'
 plot(sol.t, errors, label="error", color=[1 2], xlabel="t", ylabel="err")
 plot!(sol.t, zero(errors), ribbon=3error_estimates, label="error estimate",
-      color=[1 2], alpha=0.2)
Example block output

More about the ProbabilisticODESolution

The solution object returned by ProbNumDiffEq.jl mostly behaves just like any other ODESolution in DifferentialEquations.jl – with some added uncertainties and related functionality on top. So, sol can be indexed

julia> sol[1]2-element Vector{Float64}:
+      color=[1 2], alpha=0.2)
Example block output

More about the ProbabilisticODESolution

The solution object returned by ProbNumDiffEq.jl mostly behaves just like any other ODESolution in DifferentialEquations.jl – with some added uncertainties and related functionality on top. So, sol can be indexed

julia> sol[1]2-element Vector{Float64}:
  -1.0
   1.0
julia> sol[end]2-element Vector{Float64}:
  2.0104405118668915
@@ -93,4 +93,4 @@
  2.8014003609811094e-6
  2.7505533924662915e-6

Dense output

Probabilistic numerical ODE solvers approximate the posterior distribution

\[p \Big( y(t) ~\big|~ y(0) = y_0, \{ \dot{y}(t_i) = f_\theta(y(t_i), t_i) \} \Big),\]

which describes a posterior not just for the discrete steps but for any $t$ in the continuous space $t \in [0, T]$; in classic ODE solvers, this is also known as "interpolation" or "dense output". The probabilistic solutions returned by our solvers can be interpolated as usual by treating them as functions, but they return Gaussian distributions

julia> sol(0.45)Gaussian{Vector{Float64},PSDMatrix{Float64, Matrix{Float64}}}([-0.2773821283086793, 1.1675659430627094], 2x2 PSDMatrix{Float64, Matrix{Float64}}; R=[-3.208354407822495e-5 -4.78996984535617e-6; 0.0 1.926478845667932e-5; 0.0 1.4605790175623882e-5; 0.0 1.0839739026001365e-6; 0.0 5.6069326352075e-7; 0.0 0.0; 0.0 0.0; 0.0 0.0])
julia> mean(sol(0.45))2-element Vector{Float64}:
  -0.2773821283086793
-  1.1675659430627094

Next steps

Check out one of the other tutorials:

+ 1.1675659430627094

Next steps

Check out one of the other tutorials: