Compile CVODES functions with each Stan model

In 2017, model compilation takes ~40-50 seconds. Compiling CVODES
functions with each model adds about 7 seconds.

Closes #209.

Author: ariddell

pystan/model.py

@@ -269,6 +269,45 @@ def __init__(self, file=None, charset='utf-8', model_name="anon_model",
             s = template.safe_substitute(model_cppname=self.model_cppname)
             outfile.write(s)
 
+        ## cvodes sources
+
+        # cvodes sources are complied and linked together with the Stan model
+        # extension module. This is not ideal. In theory, build_clib could be
+        # used to build a library once and models would be complied and then
+        # linked with this library. This would save 7 or more seconds from every build.
+        # But such a strategy is frustrated by the
+        # lack of ``install_clib`` functionality in Python's distutils.
+        #
+        # TODO: numpy provides install_clib functionality, use that.
+        cvodes_src_path = os.path.join(pystan_dir, 'stan', 'lib', 'stan_math', 'lib', 'cvodes_2.9.0', 'src')
+        cvodes_sources = [
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodea.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodea_io.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_band.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_bandpre.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_bbdpre.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_dense.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_diag.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_direct.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_io.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_sparse.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_spbcgs.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_spgmr.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_spils.c'),
+            os.path.join(cvodes_src_path, 'cvodes', 'cvodes_sptfqmr.c'),
+            os.path.join(cvodes_src_path, 'nvec_ser', 'nvector_serial.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_band.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_dense.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_direct.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_iterative.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_math.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_nvector.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_spbcgs.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_spgmr.c'),
+            os.path.join(cvodes_src_path, 'sundials', 'sundials_sptfqmr.c'),
+        ]
+
         stan_macros = [
             ('BOOST_RESULT_OF_USE_TR1', None),
             ('BOOST_NO_DECLTYPE', None),
@@ -291,7 +330,7 @@ def __init__(self, file=None, charset='utf-8', model_name="anon_model",
         distutils.log.set_verbosity(verbose)
         extension = Extension(name=self.module_name,
                               language="c++",
-                              sources=[pyx_file],
+                              sources=[pyx_file] + cvodes_sources,
                               define_macros=stan_macros,
                               include_dirs=include_dirs,
                               extra_compile_args=extra_compile_args)

pystan/tests/test_cvodes.py

@@ -0,0 +1,62 @@
+import unittest
+
+import pystan
+
+
+class TestCVODES(unittest.TestCase):
+
+    def test_cvodes_program(self):
+        # from integrate_ode_bdf.stan
+        model_code = """
+            functions {
+              real[] sho(real t,
+                         real[] y,
+                         real[] theta,
+                         real[] x,
+                         int[] x_int) {
+                real dydt[2];
+                dydt[1] = y[2];
+                dydt[2] = -y[1] - theta[1] * y[2];
+                return dydt;
+              }
+            }
+            data {
+              int<lower=1> T;
+              real y0_d[2];
+              real t0;
+              real ts[T];
+              real theta_d[1];
+              real x[0];
+              int x_int[0];
+            }
+            parameters {
+              real y0_p[2];
+              real theta_p[1];
+            }
+            model {
+              real y_hat[T,2];
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_d, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_p, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_d, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_p, x, x_int);
+
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_d, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_p, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_d, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_p, x, x_int, 1e-10, 1e-10, 1e8);
+            }
+            generated quantities {
+              real y_hat[T,2];
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_d, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_p, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_d, x, x_int);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_p, x, x_int);
+
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_d, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_d, t0, ts, theta_p, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_d, x, x_int, 1e-10, 1e-10, 1e8);
+              y_hat = integrate_ode_bdf(sho, y0_p, t0, ts, theta_p, x, x_int, 1e-10, 1e-10, 1e8);
+            }
+        """
+        model = pystan.StanModel(model_code=model_code, verbose=True)
+        self.assertIsNotNone(model)