In preparation for a course I came across a model that required that the points where the ODE is evaluated be dependent on parameters, necessitating their gradients. Currently we require that the ODE time points be data, but it’s straightforward to promote them to parameters.
In other words, the derivative of the states with respect to a given time point is just ODE RHS evaluated at that time point.
Operationally this would require
promoting the time input from vector<double> to vector<T>
using meta programs to grab the value of the times where they are currently used in the code
passing the vector of ODE RHS values at each time point (n_times * n_states values) into decouple_states so that they can be included in the precomputed_gradients call
Does this all sounds good? Charles and I would like to get this in within the next week so that we can use it as a demonstration in the aforementioned course so prompt comments are appreciated.
Could well work… but shouldn’t it be possible to work around it with what we have by multiplying the ODE with
dI/dt
on both sides? So you would have to multiply the ODE RHS with the (total) derivative of you function I and then you are done. Maybe this is nonsense.
However, what certainly worries me is how the ODE integrator reacts to this. Now the time differential changes through that function. If you make it part of the ODE RHS, then all should be fine, but if you put it like that then I am not sure if the integrator is able to figure out good step sizes. Of course, you can try.
… and last: It’s going to be a hell to implement this for all the 9 cases we then end up having …
No, you cannot multiply by a Jacobian as the transformation would also involve evaluating the states at a nonlocal time.
But more importantly there are no additional states that need to be integrated or anything like that. The additional sensitivities for variations in t_{i} are literally just the evaluation of the ODE system at z(t_{i}). This includes the initial time.
This is definitely worthwhile. With pharmacometrics applications we regularly encounter situations where some of the times at which we evaluate the ODE solution are parameters.
Also the proposed approach is more attractive than one we were thinking of implementing in Torsten. It would use a change of variable like one I believe Sebastian proposed in the past, e.g., u = (t - t_i) / (t_{i+1} - t_i) where t_i and t_{i+1} are the lower and upper bounds of an integration interval, i.e., times at which ODE solutions are desired. The integration intervals are then fixed at (0, 1) and the t_i’s have to be passed to the ODE as part of the parameter vector. Not pretty…
Hi @betanalpha, is this feature still under development? I was checking the integrate_ode_rk45’s signature in develop branch and notice that “ts” is still data.
No reason we can’t generalize if you know how to calculate the derivatives and are willing to crack open the C++ code.
I don’t think anyone else has immediate plans to work on time derivatives, so I don’t think it’ll need to be coordinated with other work on the ODE solvers.