I suggest that we move beyond speculation and actually look at the code, in particular https://github.com/stan-dev/math/blob/develop/stan/math/rev/mat/functor/integrate_ode_bdf.hpp.
(Incidentally, why are the integrate_ode_pdf and CVODES utilities in mat? There are no mat dependencies and the last time I was working on the code it was in arr. Took me forever to find the files.)
Let's go to line 191. When the CVODES solver is called it tries to integrate the system to the next time point, and in doing that it will evaluate the ODE equations many times at a variety of input points. So at this point we can't guarantee the inputs.
How are the ODE equations called? For that we go up to line 140 where the function pointer to the static method is passed into CVODES' memory sandbox.
That method is defined in https://github.com/stan-dev/math/blob/develop/stan/math/rev/mat/functor/cvodes_ode_data.hpp on line 65, which then calls the
rhs method defined on line 93.
Here is the last place that we can check the output of the ODE system. The failure mode that we have been discussing is if the
ode_system_ call on line 96 resizes the
dy_dt_vec which has already been set to the correct size.
One option is to prevent the
ode_system_ functor from resizing
dy_dt_vec altogether which honestly I prefer. Still, that would require some work in the functor code to check for conflicts that would require a resize and err in an appropriate way.
Alternatively we could check the size of
dy_dt_vec immediately after line 96 and throw as necessary.
The relevant question is what is the effective API for
ode_system? This is a pattern common across lots of our code where functions don't expect the right size vector and instead of checking they just resize. I never liked this because it makes it hard for me to reason about the code, but if people prefer it then we'd have to go with the latter alternative.
And, of course, when dealing with checks that should be flagged rarely then the
unlikely macro is your friend.