Use RV at initialization in NoiseTransport

[frankschae]

Addresses https://github.com/SciML/DiffEqNoiseProcess.jl/issues/197

using DiffEqNoiseProcess, JuliaFormatter, StochasticDiffEq, Plots, Random, Distributions

# taken from the first examples about SDE in the doc (https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/#Example-1:-Scalar-SDEs)
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

function σ_lorenz(du, u, p, t)
    du[1] = 5.0
    du[2] = 5.0
    du[3] = 5.0
end

#taken from the doc about abstract noise processes
function f!(out, u, p, t, v)
    @show out # to check if noise is generated correctly
    @show v
    out[1] = sin(v[1] * t)
    out[2] = sin(t + v[2])
    out[3] = cos(t) * v[1] + sin(t) * v[2]
    nothing
end

function RV!(rng, v)
    @show v
    v[1] = randn(rng)
    v[2] = rand(rng)
    @show v
    nothing
end

rv = zeros(2)
t0 = 0.0
W = NoiseTransport(t0, f!, RV!, rv, noise_prototype=zeros(3))

# solving
prob_sde_lorenz = SDEProblem(lorenz, σ_lorenz, [1.0, 0.0, 0.0], (0.0, 1.0), noise=W)
sol = solve(prob_sde_lorenz, EM(), dt=0.01)
plot(sol, idxs=(1, 2, 3))

The values of v are indeed zero for the first run of the SDE solver. This seems to be a bug. For later runs, they do get initialized by reinit! https://github.com/SciML/StochasticDiffEq.jl/pull/502/files#diff-ca1c45d4f0bec64e8340981f6b27f4037bf6ac13c8f1af132656c485c0a8e941R422

An easy fix (?) seems to be to force initialization in the constructor. This is already done for oop/scalar processes but not for inplace/array-like random variables (see changes below). Given that the docs say:

The random variable can be either out-of-place or in-place. It is assumed it is out-of-place when the realization is a subtype of Number, and in-place, when it is a subtype of AbstractArray...

I think the fix below should be fine.

@rmsrosa could you have a look at this change as well, please? I haven't used NoiseTransport myself yet...and I am not fully confident I oversee all potential difficulties.

Checklist

• [ ] Appropriate tests were added
• [x] Any code changes were done in a way that does not break public API
• [ ] All documentation related to code changes were updated
• [x] The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
• [x] Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[rmsrosa]

Hmm, I think I did this on purpose. In the first run, it uses the value provided by rv. In the subsequent runs, it calls RV! and gets a new value. If you want to start with a random value, just call RV! early on or initialize rv with a random value instead of zero or something else.

[frankschae]

Right, I can also see this point -- and I read the docs again, and you wrote it already pretty explicit:

An optional realization rv may be given. The realization rv is used in the first time an AbstractRODEProblem is solved. Subsequent runs of the same problem will draw a different realization from the random variable RV, unless reseed is set to false

We should perhaps add this to the example with some additional comments as well? (That it's not just the rv "prototype" but really the rv for the first run.)

[rmsrosa]

Now that you drew attention to this, I come to realize it might not have been a good design option.

The whole sentence is confusing:

An optional realization rv may be given. The realization rv is used in the first time an AbstractRODEProblem is solved. Subsequent runs of the same problem will draw a different realization from the random variable RV, unless reseed is set to false. In the case of a NoiseProblem, however, a new realization will happen at the first run already, and, in this case, rv can be regarded as a realization prototype, which is necessary in the case of a random vector.

It says it is optional, but at the end it says it is necessary in the case of a random vector.

At the moment, I think the best thing is to consider rv as a prototype and always start with a new realization, and perhaps we add a keyword option to keep the given value of rvas a starting value, if one wishes to do so. What do you think?