Variadic argument lists for ODEs

[bbbales2]

Summary

See issue #1642.

See design here: stan-dev/design-docs#20

Edit: Updating description to latest changes (April 13th, 2020)

There are four new functions here:

ode_rk45, ode_rk45_tol
ode_bdf, ode_bdf_tol

I switched to shorter names because typing integrator every time seems unnecessary. It has the bonus effect of the old interfaces can remain like they were.

The _tol versions of the functions accept the relative tolerance, absolute tolerance, and max number of iterations arguments.

The interface is meant to looks like the example in the issue:

real[ , ] ode_bdf(function ode, real[] initial_state,
                  real initial_time, real[] times,
                  T1 arg1, T2 arg2, ...)

and

real[ , ] ode_bdf_tol(function ode, real[] initial_state,
                  real initial_time, real[] times,
                  real rel_tol, real abs_tol, int max_num_steps,
                  T1 arg1, T2 arg2, ...)

where the ode right hand side function looks like:

real[] ode(real time, real[] state, T1 arg1, T2 arg2, ...)

The old interfaces (intergrate_ode_bdf and integrate_ode_rk45) are meant to work as they did previously (they just wrap ode_bdf_tol and ode_rk45_tol now).

C++ Notes

The variadic arguments need to come last in the ode rhs functor. Something like:

operator()(const T0& t_in, const std::vector<T1>& y_in, std::ostream* msgs, const Args&... args)

This means all the existing ode rhs functors need changed even though the stan interface to integrate_ode_rk45 and integrate_ode_bdf aren't changing. Something like this would work:

operator()(const T0& t_in, const std::vector<T1>& y_in, std::ostream* msgs,
             const std::vector<T2>& theta, const std::vector<double>& x,
             const std::vector<int>& x_int)

The old interface looked like:

operator()(const T0& t_in, const std::vector<T1>& y_in,
             const std::vector<T2>& theta, const std::vector<double>& x,
             const std::vector<int>& x_int, std::ostream* msgs)

If this is a problem (which I assume it is for Math library versioning cuz this would break code), I believe we can add an adapter interface inside the new placeholder integrate_ode_* functions. This would mean doing code generation in Stan conditional on what function is using the functor which isn't ideal.

I removed the coupled_ode_observer class. That functionality is now split between fun/ode_store_sensitivities.hpp and functor/ode_add_time_gradients (both of these files are used in the boost and cvodes integrators).

I removed the coupled_ode_data class. Most of that just lives in the class in rev/cvodes_integrator.hpp now.

functor/coupled_ode_system.hpp has lots of changes. Generally I was trying to simplify that class.

Tests

This passes the existing tests, but will need more cause it supports a lot more stuff.

Side Effects

Did a lot of reorganizing the ODE stuff. All sorts of things could be broken.

Checklist

• Math issue Variadic argument lists for ODEs #1642

• Copyright holder: Columbia University

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • docs build, (make doxygen)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested

[wds15]

I just ported the SIR example to use variadic and I am running it. Observations:

• I am getting Exception: integrate_ode_rk45: Failed to integrate to next output time (7) in less than max_num_steps steps (in '../performance-tests-cmdstan/stat_comp_benchmarks/benchmarks/sir_variadic/sir.stan', line 49, column 4 to column 16) as an error message, which is odd since it should read "ode_rk45" for the exception to be thrown.

• Speedwise the variadic ODE seems to be of the same speed as the old codes (when using SIR with develop) => so we are all good here - I am getting for 5 replications 88 for the old and 90 for the new variadic, which looks fine to me.

We really need to motivate people to switch with 2.24+ to use variadic ODE... or they pay 30% performance penalty. However, rewriting the models like this is almost a joy to do since it is so much more natural to write beta instead of theta[1], etc.

For reference, here is the SIR variadic model:

// Simple SIR model inspired by the presentation in
// http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3380087/pdf/nihms372789.pdf

functions {
  
  // beta, water contact rate
  // kappa, C_{50}
  // gamma, recovery rate
  // xi, bacteria production rate
  // delta, bacteria removal rate
  vector simple_SIR(real t,
                    vector y,
                    real beta,
                    real kappa,
                    real gamma,
                    real xi,
                    real delta) {

    vector[4] dydt;

    dydt[1] = - beta * y[4] / (y[4] + kappa) * y[1];
    dydt[2] = beta * y[4] / (y[4] + kappa) * y[1] - gamma * y[2];
    dydt[3] = gamma * y[2];
    dydt[4] = xi * y[2] - delta * y[4];

    return dydt;
  }
}

data {
  int<lower=0> N_t;
  real t[N_t];
  vector[4] y0;
  int stoi_hat[N_t];
  real B_hat[N_t];
}

transformed data {
  real t0 = 0;
  real<lower=0> known_kappa = 1000000.0;
}

parameters {
  real<lower=0> beta;
  real<lower=0> gamma;
  real<lower=0> xi;
  real<lower=0> delta;
}

transformed parameters {
  real kappa = known_kappa;
  vector<lower=0>[4] y[N_t] = ode_rk45(simple_SIR, y0, t0, t, beta, kappa, gamma, xi, delta);
}
  
model {
  real y_hat[N_t];
  beta ~ cauchy(0, 2.5);
  gamma ~ cauchy(0, 1);
  xi ~ cauchy(0, 25);
  delta ~ cauchy(0, 1);

  stoi_hat[1] ~ poisson(y0[1] - y[1, 1]);
  for (n in 2:N_t)
    stoi_hat[n] ~ poisson(y[n - 1, 1] - y[n, 1]);

  for (n in 1:N_t)
    y_hat[n] = y[n,4];

  B_hat ~ lognormal(log(y_hat), 0.15);
  
}

I made kappa explicitly a parameter to get really the same problem to solve, but in practice one should use the known kappa directly as data, of course.

[bbbales2]

I am getting Exception: integrate_ode_rk45: Failed to integrate

Good catch. Fixed and changed the tests to catch this here

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.15	3.96	1.05	4.56% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.99	-0.52% slower
eight_schools/eight_schools.stan	0.09	0.09	1.0	0.43% faster
gp_regr/gp_regr.stan	0.19	0.19	0.99	-0.97% slower
irt_2pl/irt_2pl.stan	5.45	5.34	1.02	2.04% faster
performance.compilation	86.38	85.09	1.02	1.49% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.71	7.71	1.0	0.04% faster
pkpd/one_comp_mm_elim_abs.stan	20.19	26.66	0.76	-32.05% slower
sir/sir.stan	91.94	113.02	0.81	-22.93% slower
gp_regr/gen_gp_data.stan	0.04	0.04	0.99	-1.39% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.05	3.04	1.0	0.46% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.32	0.4	0.81	-23.47% slower
arK/arK.stan	1.78	1.77	1.0	0.49% faster
arma/arma.stan	0.72	0.73	0.99	-0.6% slower
garch/garch.stan	0.51	0.52	0.99	-0.99% slower
Mean result: 0.962306926494

---

Jenkins Console Log
Blue Ocean
Commit hash: a3f438b

---

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

[wds15]

Do we need a "Release Notes" section for this PR? I would say no, since we will anyway feature the variadic ODE boldly in the release notes.

[rok-cesnovar]

We can do that. We won't forget to mention this one. Should we just merge this? There will be a conflict with #1798 The stanc3 one is also close.

[wds15]

I think we are ready to merge. I would merge once we are confident that develop will stay "green" given that stanc3 must go in and the stan repo may need to switch to stanc3.

[rok-cesnovar]

The old interface is the same as it was and there are currently no models with the new interface in the pipeline so no reason anything would fail.

Stan repo now uses stanc3 though that is also not required now that the old interface uses the old functor. The tests that went green here have no special branches set. So merge away.

[wds15]

Great... if that's the case then we can merge. Maybe we let @bbbales2 the pleasure of doing that?

[bbbales2]

I'm gonna do #1798 first since that'll let me simplify this. Should be able to get this green again by the end of the day.

[bbbales2]

That was an easier merge than I expected. This'll need approved again (here). Should be good to go once things are green.

[rok-cesnovar]

@serban-nicusor-toptal the asnyc OpenCL are commented out now until we figure out what is going on with the gelman-group-linux setup. Just a heads up that you are aware.

[wds15]

Should we merge without that async GPU tests? Or what is the plan?

[rok-cesnovar]

Yes, if possible that we do that here. I can also open a separate PR to do only the comment-out change and then once that is merged we re-run this. I can do that if you would prefer.

That async tests should not run anyways on this PR, as this PR does not touch the /opencl folder. But I think due to the merge of recent develop Jenkins now thinks that there were changes there. I double checked that there are no changes there.

The problem is that for some reason gelman-group-linux is acting weirdly, cant find the GPU and also some header file problems with compiling TBB on distribution tests.

[wds15]

I am fine with commenting out the OpenCL stuff here on this PR... maybe create an issue for sorting out the issues so that we track this?

[rok-cesnovar]

Good call. Reopened #1940

[wds15]

LGTM

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.17	4.12	1.01	1.24% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.95	-4.76% slower
eight_schools/eight_schools.stan	0.09	0.09	1.02	1.64% faster
gp_regr/gp_regr.stan	0.19	0.19	0.99	-0.61% slower
irt_2pl/irt_2pl.stan	5.41	5.43	1.0	-0.43% slower
performance.compilation	90.79	86.01	1.06	5.26% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.38	8.39	1.0	-0.16% slower
pkpd/one_comp_mm_elim_abs.stan	21.64	26.31	0.82	-21.57% slower
sir/sir.stan	93.83	112.2	0.84	-19.58% slower
gp_regr/gen_gp_data.stan	0.04	0.04	1.0	-0.02% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.27	3.29	0.99	-0.57% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.31	0.4	0.79	-26.93% slower
arK/arK.stan	1.85	1.82	1.02	1.99% faster
arma/arma.stan	0.67	0.68	0.99	-1.14% slower
garch/garch.stan	0.61	0.6	1.02	2.37% faster
Mean result: 0.966767708464

---

Jenkins Console Log
Blue Ocean
Commit hash: bec1c2a

---

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

[wds15]

bump... @bbbales2 merge ?!