@wds15: I suspect you know and I don’t have an easy way to have a go myself quickly. So: how do the runtimes compare?
No, I was too impatient. Probably one should decrease the number of pieces and also look into model details like initialization and priors. At least I still had to use bdf solvers to make any progress, but you seem to expect a non-stiff solution. For this rather detailed work one should better have more familiarity with the model and the data as compared to myself. In any case, I hope to have put you in a helpful direction.
Just one warning: my code assumes that the initial time is equal to beta_left_t…which it was from the data I saw.
Fair enough. We can do that analysis.
The number of pieces needs to be high, unfortunately. Vaguely similarly, the initialisation and prior have been a focus for us over recent months. So, I’m not sure we can make much headway there.
I do think your pointers have been really helpful though: we want to extend this model to consider much larger numbers of coupled ODEs (about 100x as many) and need to understand how to best configure the ODE aspects (while we build a modified version of Stan that might be able to do the inference in sub-glacial timescales): Our bespoke trapezoidal solver was used here with that future application in mind. It’ll be really interesting to see how each of the candidate solvers work in the current and future context.
I’ve run the model using Brazil’s data (62 weeks)
- Bespoke Trapezoidal Solver: 1h33
- @jtimonen RK4 Solver: 2h56
- new Adjoint Solver by @wds15 : still running but from an estimate it might take around 14h
The old solution using just the RK45 and 900 ODEs would take 12 days. So, I guess this is a nice advance from 12 days to 14 hours. I guess that better analytical solutions to reduce the number of ODEs would reduce sampling time even more.
How do the different solvers compare in estimates / lp__ / diagnostics?
The Adjoint with RK45 is still running. But the trapezoidal vs RK4 CmdStan CSVs files are already in the
CoDatMo/UNINOVE_Sao_Paulo repo under folder
deaths/ is the trapezoidal and
deaths_rk4 is the RK4. Feel free to run analyses in those CSVs.
Here is another variant where I am fully exploiting the new adjoint ODE solver. The cool thing is that the many parameters being passed into the adjoint ODE solver really do not hurt that much as it looks. The sampler does make progress on my screen and this is being compiled with the latest adjoint ode cmdstan:
deaths_new_rk45-v2.stan (9.8 KB)
I am running both the
deaths_new_rk45-v3.stan with the
cmdstan-ode-adjoint-v2 in a 2018 Mac Mini i5-8500B and the
cmdstan-ode-adjoint-v3 in a 2020 Dell G5 with i7-9750H .
Will report back on times for both runs.
You need to put into make/local this line:
for good performance.
echo "CXXFLAGS+=-DSTAN_NO_RANGE_CHECKS" > make/local make build -j4
Ok so the Dell G5 is running now with the flags
This new flag is needed in Stan 2.27 onwards to turn off index checking which eats ~20% performance, but you don‘t want this flag while developing a model.
BTW, what did you mean above by „adjoint RK45“ approach of mine?
Sorry I went over the files and there was one or a proposal that you said to numerical solve the betas inside the
seeiittd ODE system using some other solver then feed it into a rk_45 or adjoint solver. I will correct the statements above… I made a whole mess of confusing names…
I think there are (at least) two bugs in your code, but I might have misunderstood something:
new_deathsby a factor of 7 due to the R script taking a mean where the name of the variables suggests a sum would be correct.
odeint_rk4returns a vector of size
num_steps+1, ie the initial conditions followed by the states at the time points.
state_estimateshould thus be of length
T+1, but it is of length
T. The latter is appropriate for the
stanode interface, but not for
For the mean vs sum of weekly deaths, you are absolutely right. Due to the nature (horrific nature) of the Brazil’s weekly deaths I thought that I was getting cumulative sums but it’s just how the weekly deaths are horrible in Brazil. Thank you for pointing it outt.
odeint_rk4 you are right, but take a look in the
transformed parameters block:
int num_steps = size(times) - 1;
I am already doing a
-1 there so the array of vectors are of length
T -1 +1 which is length
T. Does this clear things out?
Yes, I have seen this. This still leaves
state_estimate==initial_state, which is not what you want, even though this probably has very little effect. You concat eg
S = append_row(initial_state, to_vector(state_estimate[, 1]));
S = [S(0),S(0), S(1),...,S(T-2),S(T-1)]
while you actually want
S = [S(0),S(1),...,S(T-1),S(T)]
I understand now. How would I correct this?
int num_steps = size(times);
Is this appropriate?
Yes I think so. then make
state_estimate of length
T+1 and the extracted states are just the respective rows in
S = to_vector(state_estimate[, 1]);
That didn’t work:
Chain 4 Unrecoverable error evaluating the log probability at the initial value. Chain 4 Exception: subtract: m1 ((435, 1)436, 1) must match in size (in '/var/folders/5_/9mb9ggq50bl0cypgy7sd4d8r0000gp/T/RtmpK3rJQA/model-9ade69dfaaee.stan', line 211, column 2 to column 57)
for all chains.
EDIT: Maybe in line 211-212
daily_infections = S[:T+1] - S[2:] + machine_precision(); daily_deaths = D[2:] - D[:T+1] + machine_precision();
and line 216
effective_reproduction_number= (daily_infections ./ I[:T+1])*dI;
Sampling now lets see how the model will behave with Brazil’s data.