Refactor DAE solver

[yizhang-yiz]

Addressing #2643.

Summary

The refactored DAE implementation consists of

• A new IDAS interface based on RAII
• support of variadic args based on existing ODE design (wrapping/unwrapping of tuples, utilizing stan::math::for_each & stan::math::apply, etc)
• typed tests based on ODE's test fixtures. Three tests problems are
  • chemical reaction (Robertson's model)
  • predator-prey-harvesting model
  • a simple DAE that has analytical solution
• new user-facing DAE solver that follows ODE naming:
  • dae(...) with default tolerance & max num steps
  • dae_tol(...) with user-supplied controls

Tests

Typed tests:

• analytical_dae_typed_test: a simple DAE with closed-form solution
• chem_dae_typed_test: Robertson's model
• pph_dae_typed_test: predator-prey-harvesting model

The tests are done with type combinations of yy (state), yp (state derivative), and theta (parameters). Each of the three is tested with double and stan::math::var_value<double>. The validity of the results are tested against

• analytical solutions (analytical_dae_typed_test)
• finite difference gradients (chem_dae_typed_test & pph_dae_typed_test)

Side Effects

n/a

Release notes

Refactor DAE solver.

Checklist

• Math issue Refactor DAE implementation #2643

• Copyright holder: Metrum Research Group. 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)
  • dependencies checks pass, (make test-math-dependencies)
  • 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

[yizhang-yiz]

@charlesm93 interest in reviewing this?

[charlesm93]

Yes, I can take a stab at reviewing this by the end of the week.

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.57	3.56	1.0	0.2% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.96	-3.66% slower
eight_schools/eight_schools.stan	0.09	0.09	0.98	-2.0% slower
gp_regr/gp_regr.stan	0.14	0.14	0.99	-0.62% slower
irt_2pl/irt_2pl.stan	5.74	5.72	1.0	0.26% faster
performance.compilation	93.91	90.96	1.03	3.14% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.07	8.06	1.0	0.16% faster
pkpd/one_comp_mm_elim_abs.stan	31.52	30.69	1.03	2.64% faster
sir/sir.stan	121.64	119.55	1.02	1.71% faster
gp_regr/gen_gp_data.stan	0.04	0.03	1.02	2.16% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.01	3.16	0.95	-5.24% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.4	0.39	1.05	4.5% faster
arK/arK.stan	2.07	2.07	1.0	-0.25% slower
arma/arma.stan	0.23	0.24	0.97	-3.36% slower
garch/garch.stan	0.6	0.59	1.03	2.74% faster
Mean result: 1.00231834572

---

Jenkins Console Log
Blue Ocean
Commit hash: ff7815b

---

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

[charlesm93]

Overall this looks good. I'm not familiar with some of the C++ code, so I'll need to have another read and I'll likely ask more questions.

The API suggests a more general solver than I anticipated, because f doesn't need to depend linearly on y'. This means we can specify ODEs which are non-linear in y'. By contrast, functors we pass to the ODE solver have the format y' = f. Does this mean the DAE solver expands the class of ODEs we can solve?

Just checking: how are the sensitivities computed? I'm assuming we're using a forward method which computes the full Jacobian of sensitivities and that the adjoint method is currently not an option for DAEs.

I would propose to add two unit tests: one in which the solver solves a regular ODE and one in which it solves an algebraic equation. These seem to be important limiting cases. That said, I don't think these tests are required for the PR.

[yizhang-yiz]

The API suggests a more general solver than I anticipated, because f doesn't need to depend linearly on y'. This means we can specify ODEs which are non-linear in y'. By contrast, functors we pass to the ODE solver have the format y' = f. Does this mean the DAE solver expands the class of ODEs we can solve?

For the DAE the functor specifies the residual, not the right-hand-side. The residual f is supposed to be zero, which is a more general form than ODE:

f(t, y, y', ts) = 0

Just checking: how are the sensitivities computed? I'm assuming we're using a forward method which computes the full Jacobian of sensitivities and that the adjoint method is currently not an option for DAEs.

Adjoint is provided by IDAS but I choose to only implement the forward method. Since this is a refactor of things I did a few years ago, IMO it's better to leave major features to the future.

I would propose to add two unit tests: one in which the solver solves a regular ODE and one in which it solves an algebraic equation. These seem to be important limiting cases. That said, I don't think these tests are required for the PR.

Good point, I should add them to unit tests.

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.63	3.55	1.02	2.16% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.99	-1.02% slower
eight_schools/eight_schools.stan	0.08	0.09	0.97	-2.73% slower
gp_regr/gp_regr.stan	0.14	0.15	0.99	-1.2% slower
irt_2pl/irt_2pl.stan	5.68	5.73	0.99	-0.73% slower
performance.compilation	93.36	91.7	1.02	1.78% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.15	8.14	1.0	0.17% faster
pkpd/one_comp_mm_elim_abs.stan	31.82	32.02	0.99	-0.61% slower
sir/sir.stan	120.08	120.52	1.0	-0.36% slower
gp_regr/gen_gp_data.stan	0.03	0.04	0.94	-6.49% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.97	2.97	1.0	-0.17% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.42	0.92	-8.74% slower
arK/arK.stan	2.09	2.07	1.01	0.76% faster
arma/arma.stan	0.29	0.28	1.03	2.89% faster
garch/garch.stan	0.61	0.61	1.0	0.19% faster
Mean result: 0.99151693635

---

Jenkins Console Log
Blue Ocean
Commit hash: 1ebc06b

---

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

[charlesm93]

At this point my main concerns have been addressed. @yizhang-yiz you can do additional changes as you described in your comments and then merge the PR.

[yizhang-yiz]

Thanks a lot @charlesm93

[stan-buildbot]

---

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	3.55	3.65	0.97	-2.89% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.01	0.52% faster
eight_schools/eight_schools.stan	0.09	0.09	0.98	-1.95% slower
gp_regr/gp_regr.stan	0.15	0.14	1.04	3.72% faster
irt_2pl/irt_2pl.stan	5.79	5.81	1.0	-0.46% slower
performance.compilation	93.1	90.6	1.03	2.68% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	8.67	8.3	1.04	4.18% faster
pkpd/one_comp_mm_elim_abs.stan	30.7	31.72	0.97	-3.31% slower
sir/sir.stan	123.95	121.72	1.02	1.8% faster
gp_regr/gen_gp_data.stan	0.04	0.03	1.06	5.55% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan	3.01	3.08	0.98	-2.32% slower
pkpd/sim_one_comp_mm_elim_abs.stan	0.4	0.4	1.01	1.32% faster
arK/arK.stan	2.08	2.04	1.02	1.78% faster
arma/arma.stan	0.28	0.28	1.0	-0.38% slower
garch/garch.stan	0.61	0.63	0.97	-3.1% slower
Mean result: 1.00555815348

---

Jenkins Console Log
Blue Ocean
Commit hash: 56ee22d

---

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