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?