Hi,
I have been using DifferentialEquations.jl to solve a set of around 250 coupled SDEs that used to be solved with custom Euler-Maruyama. In
the model this is simulating there is a parameter, an adjacency matrix, that is constantly
changing. The changes in the parameter are recorded.
I have been trying to use a SavingCallback together with DiscreteCallback in a
CallbackSet to reproduce a parameter change that used to be implemented as a function
call after every step of the Euler-Maruyama integration procedure. Even though the
DiscreteCallback is implemented to be run after every integration step, it appears like
it acts even before the first step, messing the parameter it changes.
Here’s an example of what I mean. The animation is comparing the DiffEq implementation vs.
the legacy implementation. The parameter of interest determines the color of the dots. The
initial conditions are exactly the same.
The expected behaviour is for all points on both simulations to coincide exactly at the
beginning, and for the switches in color to match as time advances. Another clue that there is something wrong is the fact that points of a color should always be attracted to their “corner” like on the right, but on the diffeq version sometimes red dots move to the “blue” corner, and so on.
The (considerably) simplified code of both solutions look like this:
The drift function and the parameter changing function:
function point_drift(X, B)
return sum(B .* X; dims = 2)
end
# For illustrative purposes
function param_switch(X, B, dt)
avg = sum(B; dims = 2) .* dt
return B ./ avg
end
The legacy, Euler-Maruyama solver taking fixed size steps
function simulate!(X, B; Nt=200, dt=0.01)
sol = zeros(...)
params = zeros(...)
t_points = 1:(Nt - 1)
for i in t_points:
force_x = point_drift(X, B)
sol[i] .= sol + dt * force_x + sigma * sqrt(dt) * randn()
# Act on the parameters
new_B = param_switch(X, B, dt)
params[i] = new_B
B .= B
end
end
The DiffEq implementation uses the same functions
function drift(du, u, p, t)
# Indexing made simpler for illustrative purposes
du[:] .= point_drift(X, p.B)
return nothing
end
function noise(du, u, p, t)
# Indexing made simpler for illustrative purposes
du[:] .= p.σ
return nothing
end
function parameter_switch_affect!(integrator)
dt = get_proposed_dt(integrator)
u = integrator.u
p = integrator.p
new_B = param_switch(u, p.B, dt)
integrator.p.B = new_B
return nothing
end
# Callback functions
true_condition = function (u, t, integrator)
return true
end
param_switching_callback = DiscreteCallback(true_condition,
parameter_switch_affect!;
save_positions=(true, false))
save_B = function (u, t, integrator)
return integrator.p.B
end
The callbacks are combined while preparing the problem
B_cache = SavedValues(Float64, BitMatrix)
saving_callback = SavingCallback(save_B, B_cache; saveat = 0.01, save_end=false, save_start=true)
cbs = callbackset(param_switching_callback, saving_callback)
prob = SDEProblem(drift, noise, u0, (B = BitMatrix(250, 250)))
sol = solve(prob, SRIW1(); callback = cbs, alg_hints=:additive, save_everystep = true)
Without the Saving Callback I have no way of comparing the solutions, so there’s no way of only using the discrete callback. Any hints on why the switching callback is not behaving as the legacy implementation?
