Different solutions for seeded SDEProblems from ModelingToolkit and StochasticDiffEq #3971
Description
Describe the bug 🐞
It seems like this happens because MTK constructs the WienerProcess that gets passed into the SDEProblem (in calculate_noise_and_rate_prototype). Then, when the integrator is initialized, the RNG for the noise process is reseeded when the SDEProblem is constructed. However, in the non-MTK case, the noise process is constructed inside the integrator initialization (and the RNG is initialized). But as it turns out reseeding the RNG that WienerProcess uses does not actually result in an identical object as initializing it with a given seed (JuliaRandom/RandomNumbers.jl#89). So you end up having a different RNG for the noise process in the MTK vs non-MTK case. There are a few ways to make the solutions consistent (assuming the reseeding behavior is intended), but I am not sure which one would be right:
- ModelingToolkit can not pass in a noise process, so that the RNG is initialized when the integrator initialized (The noise process that MTK also seems to assume that the problem is out-of-place, which might be a bug?)
- StochasticDiffEq can reseed the RNG of the noise process that is passed in, so that the results are guaranteed to be the same as in the case when the noise process is constructed in the integrator initialization
Minimal Reproducible Example 👇
using ModelingToolkit, StochasticDiffEq
f(u, p, t) = u
g(u, p, t) = [1]
tspan = (0.0, 1.0)
prob = SDEProblem(f, g, [0.5], tspan, seed = 99)
sol = solve(prob, RKMil(), dt = 0.01, adaptive = false)
@variables x(t)
@brownians ξ
@mtkcompile sys = System([D(x) ~ x + ξ], t)
prob2 = SDEProblem(sys, [x => 0.5], tspan, seed = 99)
sol2 = solve(prob2, RKMil(), dt = 0.01, adaptive = false)
isequal(sol, sol2) # false