DiscreteCallbacks for ODEs generated from function closures behaving in a way I can't understand

General Usage
,
1

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

  1. a condition, which returns a boolean dependent upon the state of the integrator
  2. 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

Notice that 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 do_stuff with 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 do_stuff ?

Thanks a lot in advance, and sorry for the essay!

2

The former assigns the result of the righthand side to the variable on the left, i.e. it just rebinds what a is referring to. This is a different array than a. The latter writes the result into the existing array, modifying it. Since you don’t return a, the former version gets the old a passed in again, while the latter version has the updated a.

3

Indeed. but what the function returns shouldn’t depend on the initial value of a, only the updated value (which, if you look at line2 of do_stuff) is independent of the initial value. So do_stuff should have the same behaviour in either case…? a[:] was only to prevent allocations for storing a each time i do_stuff.

4

I’m not terribly familiar with how all these pieces fit together and lacking a MWE, I can only guess, sorry! Your code is not very explicit about what’s a matrix and what’s a scalar, so I don’t quite know why you need a temporary array at all.

However, looking through the documentation of the integrator, why not use SciMLBase.get_tmp_cache? According to its documentation, it should be a perfect fit for your usecase:

Returns a tuple of internal cache vectors which are safe to use as temporary arrays. This should be used for integrator interface and callbacks which need arrays to write into in order to be non-allocating. The length of the tuple is dependent on the method.

2 Likes
5

That’s absolutely fine, sorry for the lack of MWE.

Indeed, get_tmp_cache is exactly what I need! And its existence implies that it’s not safe to make your own temporary arrays. I can look into the code and figure out why so I’m all set. Thanks for the tip!

1 Like