I just realized that I didn’t mention it in my earlier post, but sometimes Stan runs and reports nothing when the ODE is clearly not specified correctly.
@wds15, I just created a new issue: https://github.com/stan-dev/math/issues/552. Please take a look at it. It doesn’t have to do with timing. This has to do with us not stomping on other memory and trying not to produce wrong results.
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.
I was suggesting checking where called, which you located as line 96 in cvodes_ode_data.hpp.
We should check the size of the return and just copy as is if it matches the size of the state vector. If not, I’m not sure how to fail in a way that’s recoverable outside of CVODES. We could try to fish a message through at that point, but I don’t think we can throw. We could fill in with dummy values, but I don’t know what that would do downstream.
P.S. It’s not so much that ode_system_ resizes dy_dt_vec, but that it will just set it to a new value of a different size—everything gets passed to the functor by constant reference.
So then there’s no point in declaring dy_dt_vec beforehand, right? In which case we can just do
void rhs(const double y[], double dy_dt[], double t) const {
const std::vector<double> y_vec(y, y + N_);
std::vector<double> dy_dt_vec = ode_system_(t, y_vec, dy_dt_vec);
if (unlikely(dy_dt_vec.size() != N_))
throw_exception();
std::copy(dy_dt_vec.begin(), dy_dt_vec.end(), dy_dt); // Now safe
}
Even better if the functor does a move!
To repeat myself, why is all the CVODES stuff in mat instead of arr?
First off, thanks @betanalpha and @Bob_Carpenter for hunting this down! I think this will address my concerns.
I think so. It should just copy over the return std::vector.
I think the function you have is close, but won’t work exactly as written. I’m making some minor adjustments:
void rhs(const double y[], double dy_dt[], double t) const {
const std::vector<double> y_vec(y, y + N_);
std::vector<double> dy_dt_vec;
ode_system_(t, y_vec, dy_dt_vec);
if (unlikely(dy_dt_vec.size() != N_))
throw_exception();
std::copy(dy_dt_vec.begin(), dy_dt_vec.end(), dy_dt); // Now safe
}
Unless you change the function signature of ode_system_'s void operator()(t, y_vec, dy_dt_vec) to something like std::vector<T> operator()(t, y_vec).
I took a look – it’s out of convenience, but not necessary. It’s used more for rectangular storage of data, which is quite natural: Jacobians. It looks like there’s at least one matrix multiplication in there.
Hi!
The intention of the ode_system object was to encapsulate the handling of the ode functor. Hence, I would place that code there.
It has not been defined, but the () operator takes a dy_dt which should already have the correct size N. As this object is only internally called within Stan, we can expect that to be correct (you cannot get N with the current design of that object).
The point of ode_system had been to also refactor the coupled_ode_system which has never happened. Hence, I am not sure whats best to do here.
The reason why cvodes stuff is in mat is that
a) I wanted to used matrix-vector multiplication
b) I can use the Eigen::Map functionality which allows me to write stuff into memory allocated from CVODES or std::vectors. So its far more flexible. My intentions for a refactor was to take advantage of that flexibility.
Sebastian
Ah… and not to forget. You also have to check the sizes in the jacobian function inside the ode_system object. There are also calls to the ODE RHS. Hence, the checking size code should go into ode_system to be consistent.
1 Like