Thank you for the pointer, @Ronis_BR. The DifferentialEquations.jl page even provides an example.
1 Like
For reference, the deprecated feature was:
The problem with DEDataArrays is that they store the discrete values in the array u. At first, this makes sense: u[1], u[2], ... are continuous, and you add u.a = .... Then you update the u.a in the integrator interface. Cool.
The issue is that it ends up being very bug prone, and it’s not necessarily possible to fix all of the bugs. The bug deals with broadcast. Let’s say you have a DEDataArray with one discrete value, so:
u = DEDataArray([1.0,2.0], 1.0)
v = DEDataArray([3.0,4.0], 1.0)
u .+ dt .* v
How is u .+ dt .* v defined? Well this is happening as a continuous update step of the ODE solver, so if the 1.0 is a discrete-time updating quantity, then it shouldn’t be changing. So the definition is to keep the discrete value and do the broadcast on the continuous parts, i.e. DEDataArray([4.0,6.0], 1.0) is the solution. Great!
But what about
u = DEDataArray([1.0,2.0], 1.0)
v = DEDataArray([3.0,4.0], 2.0)
u .+ dt .* v
This can happen after you had a callback with an affect!(integrator) = integrator.u.a = 2.0, so you updated the value. But there’s other cache array in the integrator. What you want do to is update the value of the discrete time to the “new value”… but how do you know that 2.0 is the new value?
The implementation trick that was done for a bit was to have it dependent on broadcast ordering, so here if we define broadcast to take the discrete value of the last DEDataArray in a broadcast expression, then the answer is DEDataArray([4.0,6.0], 2.0) and perfect, you propagated the updated value!
But… this should sound… scary. What if someone does a PR to the ODE solver and reorders a few arrays? Yeah… the exact ordering of the arrays in the broadcast expressions was required and if that was messed up, oops the test now fails. If that sounds brittle that’s because it is.
Ultimately, we couldn’t find a solution to this that would make it robust. And “just make sure we test it well enough so that we put the arrays in the right order” isn’t a great answer, because it can mean the ODE solver works well, but then what if they put it in a machine learning library, or an optimizer… how do we teach the whole Julia ecosystem to do this correctly? So ultimately, it was just a bug-prone way to keep and propagate the information.
The real answer from this is the following phrase:
Don’t mix implementation with representation.
You want to make it simple to represent a mixed discrete-continuous simulation where the discrete is part of the state that is being evolved, but that doesn’t mean you actually want to put the discrete values into the array that is being evolved! You want to make the representation easy, but the implementation can be different from the representation.
This of course leads to ModelingToolkit, which symbolically models the whole system to give a representation, but then it does codegen to do an optimal implementation, which may not match 1-1 the representation. For MTK, discrete values are stored inside of the parameters, because that’s the good way to implement, but not necessarily the best user-interface to the solver. So MTK hides the complexity.
Anyways, the answer of course then is that it was deprecated in favor of ModelingToolkit which can have interfaces here that are bug-free, or at least, any bugs are at least fixable in principle and so we keep being able to improve it 
. Unlike DEDataArray which is fundamentally not able to solve that issue.
2 Likes