nathanaelbosch/probnumdiffeq.jl

Probabilistic Numerical Differential Equation solvers via Bayesian filtering and smoothing

ProbNumDiffEq.jl

ProbNumDiffEq.jl provides probabilistic numerical ODE solvers to the DifferentialEquations.jl ecosystem. The implemented ODE filters solve differential equations via Bayesian filtering and smoothing. The filters 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 the ProbNumDiffEq.jl poster presentation at JuliaCon2021.

Installation

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

julia> ]
(v1.8) pkg> add ProbNumDiffEq

Example: Solving the FitzHugh-Nagumo ODE

using ProbNumDiffEq

# ODE definition as in DifferentialEquations.jl
function f(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])
end
u0 = [-1.0, 1.0]
tspan = (0.0, 20.0)
p = (0.2, 0.2, 3.0)
prob = ODEProblem(f, u0, tspan, p)

# Solve the ODE with a probabilistic numerical solver: EK1
sol = solve(prob, EK1())

# Plot the solution with Plots.jl
using Plots
plot(sol, color=["#CB3C33" "#389826" "#9558B2"])

In probabilistic numerics, the solution also contains error estimates - it just happens that they are too small to be visible in the plot above. But we can just plot them directly:

using Statistics
stds = std.(sol.pu)
plot(sol.t, hcat(stds...)', color=["#CB3C33" "#389826" "#9558B2"]
     label=["std(u1(t))" "std(u2(t))"], xlabel="t", ylabel="standard-deviation")

Contributing

Contributions are very welcome! Check the existing issues for ideas on how to contribute to the package. If you want to implement a new functionality/algorithm, open an issue to start a discussion.

Please open issues liberally! If there is anything that's unclear or doesn't work, we would very much like to know about it. This includes not just bugs and feature requests but also general questions about the software, feedback and suggestions.

Related packages

• 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.

probnumdiffeq.jl's Issues

Should we use LinearMaps.jl?

Using LinearMaps.jl could have some advantages:

• They provide both kronecker products and a sized Uniform scaling; building on their implementation could be more clean than using my own IsometricKroneckerProduct
• Our projection matrices are currently actual matrices; but their application essentially jsut selects the correct indices. This would be a usecase for a custom linear map implementation.
• We have $H = E_1 - J \cdot E_0 \in \mathbb{R}^{d \times d(q+1)}$. Currently this is implemented as an actual matrix. This makes matrix multiplication $O(d^3 (q+1)^2)$, but it could actually be $O(d^3)$ sinc e projection is cheap. LinearMaps might be a good way to implement this.
• Currently the transition matrices $A(h), Q(h)$ are actually densified for the EK1 as dense matrix matrix multiplication is simply more performant. Maybe a LinearMaps.jl implementation would not have this issue, who knows.

Add an API to specify non-zero observation noise

In our solver formulation we have an observation model of the form
$$z_i \sim \mathcal{N} \left( z_i; \dot{y}(t_i) - f(y(t_i), t_i), R \right),$$
which, together with the zero-data ${z_i = 0}$ define the information operator of the solver. Currently, we always assume $R=0$ as we are concerned with solving ODEs. But, the algorithm would in principle also work with $R>0$, which could be e.g. used when solving discretized PDEs [1]. So let's add this.

[1] "Probabilistic Numerical Method of Lines for Time-Dependent Partial Differential Equations", Krämer et al, AISTATS (2022)

Plotting a `MeanProbODESolution` doesn't work

The README example works fine:

using ProbNumDiffEq

# ODE definition as in DifferentialEquations.jl
function f(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])
end
u0 = [-1.0, 1.0]
tspan = (0.0, 20.0)
p = (0.2, 0.2, 3.0)
prob = ODEProblem(f, u0, tspan, p)

# Solve the ODE with a probabilistic numerical solver: EK1
sol = solve(prob, EK1())

# Plot the solution with Plots.jl
using Plots
plot(sol, color=["#CB3C33" "#389826" "#9558B2"])

But this doesn't:

plot(mean(sol))

Write a tutorial on uncertainties returned by our solvers

Right now we have some plots with posterior uncertainties by choosing very large steps and low order solvers, but in most applications they will barely be visible (see e.g. #25). A dedicated tutorial about uncertainties could help clarify these issues.

Topics covered could include:

• Uncertainties in general: When are they large, when are they small? How can they be visualized best?
• EK0 vs EK1 uncertainties
• What is a "diffusion model" and how do they influence the resulting uncertainties

To get started actually writing the tutorial, and for syntax etc, have a look at the current tutorials

• https://github.com/nathanaelbosch/ProbNumDiffEq.jl/blob/main/docs/src/getting_started.md
• https://github.com/nathanaelbosch/ProbNumDiffEq.jl/blob/main/docs/src/dae.md

Improve speed and reduce allocations in the ManifoldUpdate callback

Currently, the manifoldupdate function is implemented very inefficiently. By using the cache matrices already present in integ.cache, it should be possible to implement this in a much more efficient manner.

Calibration with `FixedDiffusion` and `FixedMVDiffusion` is broken with adaptive steps

Currently the calibration is changed in each perform_step! - even in those that end up rejected! This leads to bad parameters. The most likely solution is to save an increment into the cache, and only commit it in DiffEqBase.savevalues!, which we overwrite in ./src/integrator_utils.jl already to similarly handle the way the state x is updated.

This issue should explain some of the behaviour we see in #269;
and fixing this issue should hopefully fix the calibration plots for these diffusion models and make the DynamicDiffusion and FixedDiffusion much more comparable.

workflow: linear_gaussian_filtering_smoothing

While exploring this library I wondered if a workflow like the following is possible within the current implementation:

https://probnum.readthedocs.io/en/latest/tutorials/filtsmooth/discrete_linear_gaussian_filtering_smoothing.html

Basically, we would like to incorporate observations into our Bayesian estimations. Any comments, or hints are appreciated.

cheers!

Make `ClassicSolverInit` work with `SecondOrderODEProblem`s

This test should run:

ProbNumDiffEq.jl/test/secondorderode.jl

Line 46 in c89144f

@test_broken solve(prob_iip, EK1(initialization=ClassicSolverInit())) isa

`std` doesn't work

std no longer works on the solution, and throws the following error:

ERROR: DimensionMismatch: dimensions must match: a has dims (Base.OneTo(2), Base.OneTo(2)), must have singleton at dim 2

This is also in the generated docs.

Make the benchmark plots visible again

The plots in ./benchmarks/multi-language-wrappers.ipynb are not shown anymore. This should be fixed.

Some ideas for how to solve this:

• try other plot formats would help (svg vs pdf vs png)
• save the ipynb to a different file format, supported by github (though html didn't work when I last tried)
• add the benchmark to the documentation; but make sure to upload just the code and plots, and not run the code with Documenter.jl

Try adding JET.jl to the tests

JET is now part of PkgTemplates (JuliaCI/PkgTemplates.jl#432), so maybe their default JET test should be in our test suite too:

JET.test_package(ProbNumDiffEq; target_defined_modules = true)

Make the solvers work for actual `DAEProblem`s

Currently the solvers can handle DAEs, as long as they are given as ODEProblems with a mass-matrix. But in principle DAEProblems should also work, so we should have this.

The main problem is that, as we rely on OrdinaryDiffEq.solve, algorithms can't actually solve both ODEs and DAEs:
https://github.com/SciML/OrdinaryDiffEq.jl/blob/f52a126fa713f1d2d9d5744af63775c2156701ff/src/solve.jl#L74-L80

    if prob isa DiffEqBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
        error("You cannot use an ODE Algorithm with a DAEProblem")
    end


    if prob isa DiffEqBase.AbstractODEProblem && alg isa DAEAlgorithm
        error("You cannot use an DAE Algorithm with a ODEProblem")
    end

I'm not sure what the best way out here would be. A hacky solution would be to transform each ODE problem into a DAE problem like this:

function ode2dae(prob::ODEProblem)
    @assert isinplace(prob)
    function dae!(resid, du, u, p, t)
        prob.f(resid, u, p, t)
        resid .-= du
    end
    du0 = similar(prob.u0)
    prob.f(du0, prob.u0, prob.p, prob.tspan[1])
    return DAEProblem(dae!, prob.du0, prob.u0, prob.tspan, prob.p)
end

and then just make all the methods DAE solvers, but that would also be confusing to users.

The cleanest solution might be to change the OrdinaryDiffEq.jl check from alg isa OrdinaryDiffEqAlgorithm to can_solve(prob, alg) or something similar.

Remove `vars` plotting keyword in favor of `idxs`

Memory leaking when computing gradients

Hey,

I encountered some weird issue, where memory seamingly leaks during ForwardDiff.gradient(ek1_loss). I was able to boil it down to the snippet below:

using ComponentArrays
using ForwardDiff
using ProbNumDiffEq
using BenchmarkTools

function get_lv_prob(tspan=(0.0, 10.0))
    function lotka_volterra(du, u, p, t)
        du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
        du[2] = -p[3] * u[2] + p[4] * u[1] * u[2]
    end
    p = [1.5, 1.0, 3.0, 1.0]
    u0 = [1.0, 1.0]
    prob = ODEProblem(lotka_volterra, u0, tspan, p)
    return prob
end

prob = get_lv_prob((0.0, 17.0)) # no leak
prob = get_lv_prob((0.0, 20.0)) # slowly starts leaking
prob = get_lv_prob((0.0, 100.0)) # continues to leak

function ek1_loss(p, prob)
    # leaks for both EK0 and EK1
    sol = solve(prob, EK1(smooth=false), dense=false, p=p) # leak independet of order or dt
    l = sum(abs2, Array(sol)) / length(sol)
    return l
end

function test(; prob=prob)
    p = ComponentArray(prob.p)
    total_mem = Sys.total_memory() / 2^20
    i = 0
    while true
        i += 1
        ForwardDiff.gradient(p -> ek1_loss(p, prob), p)
        # GC.gc() # avoid memory leak (HOTFIX)
        free_mem = Sys.free_memory() / 2^20
        mem_usage = round((total_mem - free_mem) / total_mem * 100)
        print("\r Iteration $i, memory usage: $mem_usage%")
    end
end

test()

• I am on the latest ProbNumDiffEq v0.12.1.
• Not happening with solvers in OrdinaryDIffEq.
• Not happening with small tspans for some reason.
• Calling the garbage collector "fixes" the issue.
• Leaking happens independent of solver order, stepsize or whether EK0 or EK1 is used

Lemme know in case you need any additional info. Thanks for looking into this.

Make the solvers compatible with general `DynamicalODEProblem`s

Currently, the solvers assume that, if the given vector field f is of type DynamicalODEFunction, it comes from a SecondOrderODEProblem. Instead, we should handle the general case.

Ideally there would still be a distinction somewhere to make sure that secondorder problems are solved in the most efficient way.

Use FastLapackInterface.jl?

https://github.com/DynareJulia/FastLapackInterface.jl
We're already doing something similar in

ProbNumDiffEq.jl/src/fast_linalg.jl

Lines 57 to 91 in 70aee09

	"""
	getupperright!(A)
	Get the upper right part of a matrix `A` without allocating.
	This function is mostly there to make accessing `R` from a `QR` object more efficient, as
	`qr!(A).R` allocates since it does not use views. This function is not exported and is only
	ever used internally by `triangularize!`(@ref).
	"""
	function getupperright!(A)
	m, n = size(A)
	return UpperTriangular(triu!(A[1:min(m, n), 1:n]))
	end
	"""
	triangularize!(A; [cachemat])
	Compute `qr(A).R` in the most efficient and allocation-free way possible.
	The fallback implementation essentially computes `qr!(A).R`. But if `A` is of type
	`StridedMatrix{<:LinearAlgebra.BlasFloat}`, we can make things more efficient by calling
	LAPACK directly and using the preallocated cache `cachemat`.
	"""
	triangularize!(A; cachemat)
	function triangularize!(A; cachemat=nothing)
	QR = qr!(A)
	return getupperright!(getfield(QR, :factors))
	end
	function triangularize!(A::StridedMatrix{<:LinearAlgebra.BlasFloat}; cachemat)
	D = size(A, 2)
	BLOCKSIZE = 36
	R, _ = LinearAlgebra.LAPACK.geqrt!(A, @view cachemat[1:min(BLOCKSIZE, D), :])
	return getupperright!(R)
	end

but not quite. FastLapackInterface.jl might be the cleaner way to do this.

Write a tutorial on second-order ODEs and conserved quantities

Some statements that could be covered in the tutorial:

• "Solving second-order ODEs is easy and convenient"
• "Solving second-order ODEs is more efficient than solving first-order ODEs"; this could use ModelingToolkit.jl's ode_order_lowering functionality (see also https://mtk.sciml.ai/stable/tutorials/higher_order/)
• and just an overall demo on how to use ManifoldUpdate for energy preservation

Content-wise good starting points could be

• https://tutorials.sciml.ai/html/models/01-classical_physics.html
• https://tutorials.sciml.ai/html/models/05-kepler_problem.html

To get started actually writing the tutorial, and for syntax etc, have a look at the current tutorials

• https://github.com/nathanaelbosch/ProbNumDiffEq.jl/blob/main/docs/src/getting_started.md
• https://github.com/nathanaelbosch/ProbNumDiffEq.jl/blob/main/docs/src/dae.md

Solvers do not model uncertainty

Hi there!

Thanks for creating this toolbox!

I just tried the solvers for the first time and I'm running into some problems. My ODE is a Hodgkin-Huxley model (a model in neuroscience). Some specs (full code below):

• the model has 4 states
• it's non-linear
• it's non-homogenous

When simulating it with the EK0 solver and low tolerance, the solution is the same as if I used a simple Euler() solver. The EK0 solution is the light blue curve in the figure below. However, when I increase the tolerance in the EK0 solver, the solution changes, but all samples are the same. These are the other five traces in the figure (you can see only the dark green one because they are almost, but not bitwise, identical).

• Have you tried running the Hodgkin-Huxley model?
• Should I expect different / better results when using https://github.com/probabilistic-numerics/probnum?
• Have you observed this failure to model uncertainty?

Thanks!
Michael

Should `sol.u` be a Gaussian?

This would make it more clear what the solvers are actually doing. But, it might also make the solvers slightly less compatible with any code that just uses sol.u. I think overall we might want to do this, as then it's more consistent as all the solution values are Gaussians, whereas right now there's sol.u which is just the mean, sol.pu the marginals, but interpolations sol(t) actually return a Gaussian again. If everything is Gaussian, it's more consistent.

Question: Should the EK0 by default use a diagonal diffusion?

The issue:
The uncertainties returned by the EK0 can be a bit surprising, as it only returns sacalar covariances (s^2 * I). That is, each output dimension obtains the same amount of uncertainty, independent of its scale. This has previously lead to confusion #25.

Proposal:
Changing the default diffusion model to a DynamicMVDiffusion allows the EK0 to scale each output dimension independently. This could lead to (somewhat) more interpretable output uncertainties.

Related PR: #65 applies the proposed change.

Make partial observations in the data likelihoods work with the DiagonalEK1

This should probably be done after #302: If we handle the solution object better and give it the actual covariance structure it deserves, then implementing this might be much easier.

Add to DiffEqDocs?

It would be good to add these methods to DiffEqDocs in some way.

Make the DiagonalEK1 work for SecondOrderODEs

In principle this is quite straight-forward, but the fact that SecondOrderODEs by default always stack [du; u] is a bit annoying for us as that results in weird covariance structure. For example, if we have isometric Kronecker covariances of the form $I_d \otimes \Sigma$ in our state-space model, then [du; u] has covariances of the form $\tilde{Sigma} \cdot I_d$, i.e. the order of the Kronecker changed. For BlockDiagonals its a bit weirder: Instead of having a BlockDiagonal matrix, which has dense blocks on its diagonal, we have a BlockMatrix densely filled with Diagonals.

So basically, properly resolving this issue might mean implementing a way to efficiently switch between the two orderings. Or in the extreme, it migh even mean changing the ordering we have so far, which should then make this and related things easier.

Reduce compilation times of the package

The package currently takes quite a long time to compile, which is especially frustrating for new Julia users coming from Python. Thus, it would be great to reduce the compilation times.

Bump the OrdinaryDiffEq.jl compat to v6

I removed the v6 compat since I ran into some issues with some changes, I believe to the implicit solver types, so make sure to try out the individual versions and see if the tests run.

Use IntegratedWienerProcesses.jl?

https://github.com/filtron/IntegratedWienerProcesses.jl

It includes functionality to compute the cholesky factor of the transition noise covariance directly. So, we should probably use it!

Make sure to benchmark the change to check that there are no major regressions in performance.

Use Gaussian as Initial Condition

Is it possible to use a Gaussian as an intial condition?

For example something like this?

using LinearAlgebra, Statistics, Distributions
using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots, ODEInterfaceDiffEq
using ModelingToolkit
using ProbNumDiffEq

function osscilator!(ddu, du, u, p,t)
    ddu .= - p .* u
end

# osscilator
u0 = [1]
du0 = [1]
p = 100
T = 3
t = 0:0.1:T
prob = SecondOrderODEProblem(osscilator!, du0, u0, (0,T),p)
@time sol = solve(prob, EK0(;smooth=true), abstol=1e-1, reltol=1e-1)
plot(sol)

Where I would like to use
u0 = [Gaussian(1,1)] or something similar as an intial condition instead.

Currently it does not seem to be suported: ERROR: MethodError: no method matching zero(::Type{Any}) but it would be very useful. For example when the initial condition is a result of a measurement and is not known with infinite precision.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

Add one or more energy conservation benchmarks

Something similar to these SciML Benchmarks would be nice to have:

• Henon-Heiles Energy Conservation: https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DynamicalODE/Henon-Heiles_energy_conservation_benchmark/
• Quadruple Boson Energy Conservation: https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark/#Quadruple-Boson-Energy-Conservation
• Single Pendulum Comparison (currently broken): https://docs.sciml.ai/SciMLBenchmarksOutput/stable/DynamicalODE/single_pendulums/#Solving-single-pendulums-by-DifferentialEquations.jl

This part of the ProbNumDiffEq.jl documentation shows how we can do energy conservation, and that's what we'd want to evaluate and compare to projection methods and symplectic methods.

Use PSDMatrices.jl instead of the small `SquarerootMatrix` implementation in `./src/squarerootmatrix.jl`

Essentially, the functionality implemented in ./src/squarerootmatrix.jl should be moved into a dedicated package.

Some possible modifications / improvements from the current code:

• only storing the square-root
  I'm 80% sure that this would overall improve performance in ProbNumDiffEq.jl, but this might be application dependent.
• arithmetic on square-root matrices

Related examples:

• PDMats.jl, a largely used package to for positive definite matrices
• My try at PSDMatrices.jl
• The PSD type in GaussianDistributions.jl (link)

Update:
PSDMatrices.jl is under active development again and received some updates. So, let's use it here for all our PSD matrix needs!

Solving Forward and Inverse Problems

I was wondering if ProbNumDiffEq can be used to solve forward and inverse problems in one go as described in the following paper?

https://arxiv.org/abs/2103.10153

Or is there a way to explicitly prescribe the latent variables we want to estimate apart from the solution to the ODE?

Remove `smooth!` and replace it's usage with backward kernels in the classic solver init

The only place where smooth! is still used is src/init/classicsolverinit.jl. Use backward kernels instead and remove smooth! completely from the package.

Add unit tests for the functions implemented in `scr/fast_linalg.jl`

Try using TaylorDiff.jl instead of TaylorIntegration.jl

Not quite sure what to gain here, but at the very least the dependency could be more lightweight.

The part that would be changed is this one:

ProbNumDiffEq.jl/src/initialization/taylormode.jl

Lines 35 to 49 in 70b6d57

	"""
	Compute initial derivatives of an IIP ODEProblem with TaylorIntegration.jl
	"""
	function taylormode_get_derivatives(u, f::SciMLBase.AbstractODEFunction{true}, p, t, q)
	tT = Taylor1(typeof(t), q)
	tT[0] = t
	uT = similar(u, Taylor1{eltype(u)})
	@inbounds @simd ivdep for i in eachindex(u)
	uT[i] = Taylor1(u[i], q)
	end
	duT = zero(uT)
	uauxT = similar(uT)
	TaylorIntegration.jetcoeffs!(f, tT, uT, duT, uauxT, p)
	return [evaluate.(differentiate.(uT, i)) for i in 0:q]
	end

Currently I'm not sure how to best get exactly this functionality from TaylorDiff.jl. If you do, please comment!

Fix the usage of `idxs` in plot recipes

Not all of the following work the way they should

plot(sol; idxs=(2,1))
plot(sol; idxs=[(0,1), (1,2)])
plot(sol; idxs=(0,1,2))
plot(sol; idxs=(1,1,2))
plot(sol; idxs=[(1,1,2),(0,1,2)])

We should make sure they all work.

Use FiniteHorizonGramians.jl

• github: https://github.com/filtron/FiniteHorizonGramians.jl
• arXiv paper: https://arxiv.org/abs/2310.13462

TL;DR: Instead of doing what we do in src/priors/..., we could compute A and Q_chol with FiniteHorizonGramians.jl. This is probably particularly relevant for the exponential integrator case, so whenever doing this, make sure to benchmark.

Fix type piracies

Current output:

julia> Aqua.test_piracies(ProbNumDiffEq)
Possible type-piracy detected:
[1] *(M::AbstractMatrix, g::Gaussian{VM, PSDMatrix{T, S}} where {T, S, VM<:AbstractVecOrMat{T}}) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:25
[2] X_A_Xt(A, X) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/ProbNumDiffEq.jl:37
[3] copy(K::Kronecker.KroneckerProduct) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/kronecker.jl:1
[4] copy(P::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:4
[5] copy!(dst::Gaussian, src::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:6
[6] isapprox(M1::PSDMatrix, M2::PSDMatrix; kwargs...) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/filtering/markov_kernel.jl:48
[7] kwcall(::Any, ::typeof(isapprox), M1::PSDMatrix, M2::PSDMatrix) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/filtering/markov_kernel.jl:48
[8] mean(s::StructArray{Gaussian{VM, PSDMatrix{T, S}} where VM<:AbstractVecOrMat{T}} where {T, S}) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:35
[9] ndims(g::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:17
[10] recursivecopy(P::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:11
[11] recursivecopy!(dst::Gaussian, src::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:12
[12] show(io::IO, g::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:13
[13] similar(P::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:5
[14] size(g::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:16
[15] std(g::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:19
[16] var(g::Gaussian) @ ProbNumDiffEq ~/.julia/dev/ProbNumDiffEq/src/gaussians.jl:18