Compling error with ODE functions

General
1

I’m trying use Stan’s numerical integrator and replicating the stochastic differential equation model from http://ctsm.info/ (example building 2). Here’s the Stan code:

functions {
  // the TiTe state space model 
  real[] TiTe_ODE(real t, real[] X, real[] parms, real[] rdata,	int[] idata) {
    // the state variables:
    // X[1]: internal temperature
    // X[2]: envelope (the thermal mass) temperature
    real Ci = parms[1];   // capacitance internal
    real Ce = parms[2];   // capacitance envelope
    real Rie = parms[3];  // resistance internal-envelope 
    real Rea = parms[4];  // resistance envelope-ambient 
    real Aw = parms[5];   // effective glazing area (aka gA)
    
    real Ta = rdata[1];   // ambient air temperature
    real Ps = rdata[2];   // thermal power from solar
    real Ph = rdata[3];   // thermal power from heaters
    
    real dXdt[2]; 
    
    // the differential equations
    dXdt[1] =  1/(Ci*Rie)*(X[2]-X[1]) + Aw/Ci*Ps + 1/Ci*Ph;
    dXdt[2] =  1/(Ce*Rie)*(X[1]-X[2]) + 1/(Ce*Rea)*(Ta-X[2]);
    
    return dXdt;
    
  }
  
  // needed to not get syntax error on inital time for integrate_ode_rk45, some referencing/data issue
  real[] TiTe_integrate(real[] init, real t0, real[] time, real[] parms, real[] rdata,int[] idummy) {
    real temp[1, 2] = integrate_ode_rk45(TiTe_ODE, init, t0, time, parms, rdata, idummy);
    real res[2] = to_array_1d(temp); 
    return res;
  }
  
  
  
  matrix TiTe(real[] time, real[] Ta, real[] Ps, real[] Ph,
              real Ci, real Ce, real Rie, real Rea, real Aw) {
    int nX = 2;            // number of state variables
    int nt = size(time);   // number of time steps
    matrix[nt, nX] result;
    real parms[5] = to_array_1d([Ci, Ce, Rie, Rea, Aw]); // parameters have to be in array
    int idata[0];     // integer data, not used
    //init = rep_array(0, nT);
    
    real t0 = time[1];
    real init[nX];
    init[1] = 15;  // init internal temperature
    init[2] = 10;  // init same as CTSM-model
    // init[2] = (init[1] +  Ta[1]) / 2;  // envelope temp as average of internal & ambient

    for(i in 2:nt){
      real rdata[3] = to_array_1d([Ta[i], Ps[i], Ph[i]]);  // real data, have to be in array
      init = TiTe_integrate(init, t0, time[i:i], parms, rdata, idata);
      for(j in 1:nX) result[i, j] = init[j]; 
      t0 = time[i];
    }
    return result;
  }  
}

data {
  int nt;
  real time[nt];
  real Ta[nt];
  real Ps[nt];
  real Ph[nt];
  
  // observation
  real yTi[nt];
}

parameters {
  // R-C network parameters
  real<lower=0.0, upper=20> Ci;
  real<lower=0.0, upper=20> Ce;
  real<lower=0.0, upper=50> Rie;
  real<lower=0.0, upper=50> Rea;
  real<lower=0.0, upper=40> Aw;
  real<lower=0> sigmaTi;
  //real<lower=0> sigmaTe;
}
transformed parameters {
  // the state variables
  matrix[nt, 2] X;  // X[, 1]: Ti (internal temperature), X[, 2]: Te (envelope temperature)
  X = TiTe(time, Ta, Ps, Ph, Ci, Ce, Rie, Rea, Aw);
}
model {
  // observed data likelihood
  yTi ~ normal(X[, 1], sigmaTi); 
}

Exposing the functions with rstan::expose_stan_functions works fine and calling the TiTe-function from R with data and parameters return expected results.

But, trying to compile with stan_model returns following error: Compilation ERROR, function(s)/method(s) not created!

Any suggestions on what is going wrong are welcome!

Regards,
Lukas Lundström

2

The compilation error is ultimately

/usr/local/lib/R/site-library/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:70:5: note: candidate function not viable: no known conversion from 'const stan::math::var' to 'double' for 3rd argument
    integrate_ode_rk45(const F& f,

By defining real t0 = time[1] on line 45, you make it a var. Changing line 53, so that the second argument is time[1] helps but then you run into a similar error with the 6th argument.

3

Thanks, that helped.

Seems like the 6th argument (the real data) to integrate_ode_rk45 have to be declared in data or transformed data block. While the 7th argument (the integer data) can be declared later.

So to get real data to the function to be integrated, all of it has to be passed in one argument, eg.:

transformed data {
  real rdata[nt*3] = to_array_1d(append_col(append_col(Ta, Ps), Ph));
}

And then the integer data argument can be used to pass indexes to the real data in the rdata. Is this how it is supposed to be done or is there a more elegant solution? Anyway, here’s a working example on passing data to the Stan integrator:

functions {
  // the TiTe state space model 
  real[] TiTe_ODE(real t, real[] X, real[] parms, real[] rdata, int[] idata) {
    // the state variables:
    // X[1]: internal temperature
    // X[2]: envelope (the thermal mass) temperature
    real Ci = parms[1];   // capacitance internal
    real Ce = parms[2];   // capacitance envelope
    real Rie = parms[3];  // resistance internal-envelope 
    real Rea = parms[4];  // resistance envelope-ambient 
    real Aw = parms[5];   // effective glazing area (aka gA)
    
    real Ta = rdata[idata[2] + 0];          // ambient air temperature
    real Ps = rdata[idata[2] + idata[1]];   // thermal power from solar
    real Ph = rdata[idata[2] + idata[1]*2]; // thermal power from heaters
    
    real dXdt[2]; 

    // the differential equations
    dXdt[1] =  1/(Ci*Rie)*(X[2]-X[1]) + Aw/Ci*Ps + 1/Ci*Ph;
    dXdt[2] =  1/(Ce*Rie)*(X[1]-X[2]) + 1/(Ce*Rea)*(Ta-X[2]);
    
    return dXdt;
    
  }

  
  matrix TiTe(real[] time, real[] rdata, 
              //real[] Ta, real[] Ps, real[] Ph,
              real Ci, real Ce, real Rie, real Rea, real Aw) {
    int nX = 2;            // number of state variables
    int nt = size(time);   // number of time steps
    matrix[nt, nX] result;
    real parms[5] = to_array_1d([Ci, Ce, Rie, Rea, Aw]); // parameters have to be in array
    int idata[2]; 
    
    real temp[1, 2];
    real init[nX];
    init[1] = 20;  
    init[2] = 7;  
    
    idata[1] = nt;   // used for indexing rdata
    
    result[1, 1] = init[1];
    result[1, 2] = init[2];
    for (i in 2:nt) {
      idata[2] = i;   // used for indexing rdata
      temp = integrate_ode_rk45(TiTe_ODE, init, time[i-1], time[i:i], parms, rdata, idata);
      init = temp[1, ];
      result[i, ] = [init[1], init[2]];
    
    }
    return result;
  }  
}

data {
  int nt;
  real time[nt];
  row_vector[nt] Ta;
  row_vector[nt] Ps;
  row_vector[nt] Ph;
  // observation
  real yTi[nt];
}

transformed data {
  real rdata[nt*3] = to_array_1d(append_col(append_col(Ta, Ps), Ph));
}

parameters {
  // R-C network parameters
  real<lower=0.1, upper=20> Ci;
  real<lower=0.1, upper=20> Ce;
  real<lower=0.1, upper=50> Rie;
  real<lower=0.1, upper=50> Rea;
  real<lower=0.1, upper=40> Aw;
  real<lower=0.0> sigmaTi;
  //real<lower=0> sigmaTe;
}
transformed parameters {
  // the state variables
  matrix[nt, 2] X;  // X[, 1]: Ti (internal temperature), X[, 2]: Te (envelope temperature)
  X = TiTe(time, rdata, Ci, Ce, Rie, Rea, Aw);
  
}
model {
  // observed data likelihood
  yTi ~ normal(X[, 1], sigmaTi); 
}
1 Like
4

the Stan reference manual section “Ordinary Differential Equation Solvers” does spell out the argument restrictions quite clearly.
@bgoodri was this problem not caught by the parser because rstan::expose_stan_functions was used?

5

The opposite: expose_stan_functions forces all arguments to be data so it worked. The original model would parse but not compile when trying to estimate the parameters.

6

I see. https://github.com/stan-dev/stan/issues/2369 will be in 2.17 and that should take care of this problem.

7

Ok, thanks. I’ll try it when v 2.17 is out.