# Compling error with ODE functions

#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);
}

``````

#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.