Bug in expose_stan_functions for integrate_ode_bdf

Hi,

I found that changing:

integrate_ode_rk45

to

integrate_ode_bdf

Means that the following code prevents,

expose_stan_functions

from compiling/working. (However the Stan model itself can be compiled and run.) The top of the error messages are as follows,

In file included from file775e09fb39.cpp:5:
In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rstan/include/boost_not_in_BH/boost/fusion/support/unused.hpp:10:
In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/fusion/support/config.hpp:11:
In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config.hpp:39:
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include/boost/config/compiler/clang.hpp:196:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
#  define BOOST_NO_CXX11_RVALUE_REFERENCES
          ^
<command line>:6:9: note: previous definition is here
#define BOOST_NO_CXX11_RVALUE_REFERENCES 1
        ^
In file included from file775e09fb39.cpp:8:
In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:56:
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat/functor/ode_system.hpp:79:35: error: no matching function for call to object of type 'const deriv_aslanidi_functor__'
          vector<var> dy_dt_var = f_(t, y_var, theta_, x_, x_int_, msgs_);
                                  ^~
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat/functor/cvodes_ode_data.hpp:104:21: note: in instantiation of function template specialization 'stan::math::ode_system<deriv_aslanidi_functor__>::jacobian<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Map<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0, Eigen::Stride<0, 0> > >' requested here
        ode_system_.jacobian(t, y_vec, fy, Jy_map);
                    ^
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat/functor/cvodes_ode_data.hpp:89:30: note: in instantiation of member function 'stan::math::cvodes_ode_data<deriv_aslanidi_functor__, double, double>::dense_jacobian' requested here
        return explicit_ode->dense_jacobian(NV_DATA_S(y), J, t);
                             ^
/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/stan/math/rev/mat/functor/integrate_ode_bdf.hpp:157:57: note: in instantiation of member function 'stan::math::cvodes_ode_data<deriv_aslanidi_functor__, double, double>::dense_jacobian' requested here
                                             &ode_data::dense_jacobian),
                                                        ^
file775e09fb39.cpp:259:35: note: in instantiation of function template specialization 'stan::math::integrate_ode_bdf<deriv_aslanidi_functor__, double, double>' requested here
            stan::math::assign(I, integrate_ode_bdf(deriv_aslanidi_functor__(), rep_array(X0,1), t0, ts, theta, to_array_1d(append_row(to_vector(ts),to_vector(V))), x_i, pstream__));
                                  ^
file775e09fb39.cpp:219:5: note: candidate function not viable: no known conversion from 'vector<stan::math::var>' to 'const vector<double>' for 2nd argument
    operator()(const double& t,
    ^
1 warning and 1 error generated.
make: *** [file775e09fb39.o] Error 1
clang++ -I/Library/Frameworks/R.framework/Resources/include -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rstan/include/boost_not_in_BH" -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DBOOST_NO_CXX11_RVALUE_REFERENCES -I/usr/local/include -I/usr/local/include/freetype2 -I/opt/X11/include  -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/Rcpp/include" -I"/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rstan/include" -I"/private/var/folders/f2/b6912q052s38c1tl7ymgk_s40000gn/T/Rtmp4B22mo/sourceCpp-x86_64-apple-darwin13.4.0-0.12.8"   -fPIC  -Wall -mtune=core2 -g -O2  -c file775e09fb39.cpp -o file775e09fb39.o
Error in Rcpp::sourceCpp(code = paste(lines, collapse = "\n"), env = env) : 
  Error 1 occurred building shared library.
In addition: Warning message:
In readLines(file, warn = TRUE) :
  incomplete final line found on 'aslanidi_1.stan'
Here is the C++ code that does not compile. Please report bug.

The Stan model is as follows,

functions{
  int find_interval_elem(real x, vector sorted, int start_ind){
    int res;
    int N;
    int max_iter;
    real left;
    real right;
    int left_ind;
    int right_ind;
    int iter;

    N = num_elements(sorted);

    if(N == 0) return(0);

    left_ind  = start_ind;
    right_ind = N;

    max_iter = 100 * N;
    left  = sorted[left_ind ] - x;
    right = sorted[right_ind] - x;

    if(0 <= left)  return(left_ind-1);
    if(0 == right) return(N-1);
    if(0 >  right) return(N);

    iter = 1;
    while((right_ind - left_ind) > 1  && iter != max_iter) {
      int mid_ind;
      real mid;
      // is there a controlled way without being yelled at with a
      // warning?
      mid_ind = (left_ind + right_ind) / 2;
      mid = sorted[mid_ind] - x;
      if (mid == 0) return(mid_ind-1);
      if (left  * mid < 0) { right = mid; right_ind = mid_ind; }
      if (right * mid < 0) { left  = mid; left_ind  = mid_ind; }
      iter = iter + 1;
    }
    if(iter == max_iter)
      print("Maximum number of iterations reached.");
    return(left_ind);
  }
  
  real[] deriv_aslanidi(real t, real[] I, real[] theta, real[] x_r, int[] x_i){
    
    int aLen = x_i[1];
    vector[aLen] ts = to_vector(x_r[1:aLen]);
    vector[aLen] V = to_vector(x_r[(aLen+1):(2*aLen)]);
    int aT = find_interval_elem(t, ts, 1);
    real aV = (aT==0) ? V[1] : V[aT];
    
    real xtau = theta[1] / (1 + exp(aV/ theta[2])) + theta[3];
    real xinf = 1 / (1 + exp(-(aV + theta[4]) / theta[5]));
    real rinf = 1 / (1 + exp((aV + theta[6]) / theta[7]));
    real dydt[1];
    dydt[1] = (xinf - I[1]) / xtau;
    return dydt;
  }
  
  vector solve_aslanidi_forced_ode(real[] ts, real X0, real[] theta, real[] V, real t0){
    int x_i[1];
    real I[size(V),1];
    x_i[1] = size(V);

    I = integrate_ode_bdf(deriv_aslanidi, rep_array(X0, 1), t0, ts, theta, to_array_1d(append_row(to_vector(ts), to_vector(V))), x_i);
    return(to_vector(I[,1]));
  }
}

data{
  int N;
  real V[N];
  real I[N];
  real ts[N];
  real t0;
}

transformed data {
  int x_i[0];
}


parameters{
  real<lower=0> p1;     // ms
  real<lower=0> p2;     // mV
  real<lower=0> p3;     // ms
  real<lower=0> p4;     // mV
  real<lower=0> p5;     // mV
  real p6;              // mV
  real<lower=0> p7;     // mV
  real<lower=0,upper=1> X0;
  real<lower=0> sigma;
}

transformed parameters{
  real theta[7];
  theta[1] = p1;
  theta[2] = p2;
  theta[3] = p3;
  theta[4] = p4;
  theta[5] = p5;
  theta[6] = p6;
  theta[7] = p7;
}

model{
  // solve ODE using stiff solver
  vector[N] I_int;
  I_int = solve_aslanidi_forced_ode(ts, X0, theta, V,-0.1);
  
  // likelihood
  for(i in 1:N){
    I[i] ~ normal(I_int[i],sigma);
  }
  
  //priors
  p1 ~ normal(900,500);
  p2 ~ normal(5,1);
  p3 ~ normal(100,10);
  p4 ~ normal(0.1,0.02);
  p5 ~ normal(12.25,3);
  p6 ~ normal(-5.6,1);
  p7 ~ normal(20.4,3);
  sigma ~ normal(1,0.1);
}

I think this is one for @bgoodri

We should have this sorted soon with standalone function compilation.