Either I’m fundamentally misunderstanding function closures, or DiscreteCallbacks is doing something strange. Either way, input appreciated!
I’m using OrdinaryDiffEq and DiffEqCallbacks to simulate an ODE with (multiple) discrete callbacks. The problem is that the ODE has a slightly different solution each time I run it, if i do x. (it runs reliably and correctly if I don’t do x).
Let me explain x, which has to do with function closures and captured variables, and may not require an understanding of DiffEqCallbacks.
A discrete callback has
- a condition, which returns a boolean dependent upon the state of the integrator
- an affect, which changes the integrator state in a particular way if the condition returns true.
For example, a simulation of a bouncing ball might have a callback with a condition that returns true if the ball’s height is <= 0 (ie hits the floor). And an affect that changes the dynamics appropriately in this context. See Event Handling and Callback Functions · DifferentialEquations.jl for more info if you need.
My discrete callback looks like
function get_callback(c::MyCustomType) condition = build_condition(c::MyCustomType) affect! = build_affect(c) cb = DiscreteCallback(cond, affect!) end
Now I’ll show you the programming pattern for
build_condition. The condition has to be a function dependent upon (u,t,integ), where u is the current state of the system, t is the time, and integ is the integrator interface
function build_condition(c::MyCustomType) a = template_to_avoid_allocation(c) function condition(u,t,integ) b = do_stuff(u,t,integ,a) b > 1e-3 ? return true : return false end end
do_stuff depends upon
a, a captured variable. Actually, the output should be independent of
a. I only pass
a to avoid some allocations, by holding a temporary array in
a[:] as you can see below (and here I have my actual code):
function do_stuff(u,t,integ,a) λ, μ1, θ, θ₀ = get_from(integ) a[:] = (-λ + μ1 * (θ - θ₀)) / (2 * μ2) a /= (sqrt(sum((a).^2))) integ.u[N + 1:end] = -2 * μ2 * dθ return integ end
And herein lies the problem. Everything works fine if I replace the second line of
a = ... instead of
a[:] = .... I don’t understand why the ODE behaviour changes every time if I use the code as is, above, with
a[:] = ....
I have exactly the same programming pattern and the same problem with my
build_affect function as well.
Is there something about the behaviour of function closures I’m fundamentally misunderstanding? Or is there an issue with DiffEqCallbacks, that means this ‘dirty state’ a is not allowable, even if it doesn’t affect function output for
Thanks a lot in advance, and sorry for the essay!