CoDatMo Liverpool & UNINOVE models (slow ODE implementation and Trapezoidal Solver)

@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.

Thank you.

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 SEIR-model/results/: deaths/ is the trapezoidal and deaths_rk4 is the RK4. Feel free to run analyses in those CSVs.

1 Like

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 deaths_new_rk45-v2.stan with 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

Looks good

1 Like

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:

  • daily_deaths underestimates new_deaths by 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_rk4 returns a vector of size num_steps+1, ie the initial conditions followed by the states at the time points. state_estimate should thus be of length T+1, but it is of length T. The latter is appropriate for the stan ode interface, but not for odeint_rk4.
1 Like

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.

Regarding the 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[1]==initial_state, which is not what you want, even though this probably has very little effect. You concat eg

  S = append_row(initial_state[1], to_vector(state_estimate[, 1]));

such that

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 state_estimate, e.g.

  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.