@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 nonstiff 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 subglacial 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 SEIRmodel/results/
: 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:
https://github.com/standev/cmdstan/suites/2733052299/artifacts/60343510
deaths_new_rk45v2.stan (9.8 KB)
I am running both the deaths_new_rk45v3.stan
with the cmdstanodeadjointv2
in a 2018 Mac Mini i58500B and the deaths_new_rk45v2.stan
with cmdstanodeadjointv3
in a 2020 Dell G5 with i79750H .
Will report back on times for both runs.
You need to put into make/local this line:
CXXFLAGS+=DSTAN_NO_RANGE_CHECKS
for good performance.
Like?
echo "CXXFLAGS+=DSTAN_NO_RANGE_CHECKS" > make/local
make build j4
Looks good
Ok so the Dell G5 is running now with the flags
CXXFLAGS+=DSTAN_NO_RANGE_CHECKS
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
underestimatesnew_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 sizenum_steps+1
, ie the initial conditions followed by the states at the time points.state_estimate
should thus be of lengthT+1
, but it is of lengthT
. The latter is appropriate for thestan
ode interface, but not forodeint_rk4
.
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(T2),S(T1)]
while you actually want
S = [S(0),S(1),...,S(T1),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/model9ade69dfaaee.stan', line 211, column 2 to column 57)
for all chains.
EDIT: Maybe in line 211212
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.