Ode_rk45 signature

Hi,

I’m trying to model activity propagation in a network with N nodes for N_TS time points that spreads as described by the following evolution equation:

\frac{dx_i}{dt}\tau_i = -x_i(t) + s\phi\big(x_i(t)\big) + g\Bigg(\sum_{j\neq i}^{N} W_{ij}\phi\big(x_j(t)\big)\Bigg) + I_i(t)

This is the first time I’m trying to solve an ODE in a Stan model so it’s very possible that I’ve made other mistakes as well. I haven’t even gotten to the point of testing the actual model yet because I can’t compile it due to the following error:

Compiling Stan program...
/
Semantic error in '/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmphjhKZn/model-e51f176820ee.stan', line 120, column 22 to column 82:
   -------------------------------------------------
   118:    // states x are estimated as parameters z with some uncertainty. 
   119:    // these estimated parameters are used in the model description to relate them to measured data
   120:    vector[N] z[N_TS] = ode_rk45(dx_dt, z_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
                               ^
   121:  }
   122:  
   -------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45'. Available signatures:
(real, vector, ...) => vector, vector, real, real[],  ...
.

Instead supplied arguments of incompatible type:
(real, vector, int, vector[], vector, vector, real, real, real, real) => vector, vector, real, real[], int, vector[], vector, real[], real, real, real, real.

The supplied arguments of type real and vector look like they should agree with the required signature for the ode_rk45 function so I’m not sure how to debug this. As I said, it is entirely possible that the error lies somewhere else but I’m not sure where to start. Any help would be appreciated.

Here’s the full model:



functions {
  
  real phi(real a){
    real out = (exp(2*a)-1)/(exp(2*a)+1);
    return out;
  }
  
  // From
  // https://aws1.discourse-cdn.com/standard14/uploads/mc_stan/original/1X/25b0af2e6aba491a267c45fa090c6374ad57c4b4.pdf
  // https://discourse.mc-stan.org/t/zero-chain-movement-mixing-in-forced-ode-model/921
  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);
  }
  
  vector dx_dt(real t,       // time
               vector x,     // system state {1 x N - network activity for each time point}
               int N,
               vector[] N_t, //vector array
               vector I,
               vector ts,
               real s,
               real g,
               real b, 
               real tau) {
  
    vector[N] res;
    
    int d = find_interval_elem(t, ts, 1);
    
    for (i in 1:N){
      res[i] = (-x[i] + s * phi(x[i]) + g *N_t[d, i] + b * I[d])/tau;
    }
    
    // returns change for one time point for all nodes
    return  res;
  }
}

data {
  int<lower = 0> N_TS;           // number of measurement times
  int<lower = 0> N;             // number of nodes
  real t0;
  real ts[N_TS];                 // measurement times > 0
  vector[N] y_init;             // initial measured activity level
  
  // Arrays are declared by enclosing the dimensions in square brackets following the name of the variable.
  // this is an array of length N_TS consisting of (column) vectors of length N
  // Each array element is like a row
  vector[N] y[N_TS]; // measured activity level with nodes in cols and timepoints in rows
  matrix[N, N] W;// adjacency matrix
  
  vector[N_TS] I; // task stimulation
}

transformed data{
  vector[N] N_t[N_TS];
  
  for(i in 1:N){
    for(t in 1:N_TS){
      //row of incoming connections to node i of length n
      //multiplied with activity (column) vector of length n for all nodes at time t
      N_t[t, i] = dot_product(W[i,], y[t]);
    }
  }
}

parameters {
  real<lower = 0> s;   // self coupling
  real<lower = 0> g;  // global coupling
  real<lower = 0> b;  // task modulation
  real<lower = 0> tau; // time constant
  real<lower = 0> sigma;   // measurement error
  vector[N] x_init; 
}


transformed parameters {
  // states x are estimated as parameters z with some uncertainty. 
  // these estimated parameters are used in the model description to relate them to measured data
  vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
}

model {
  s ~ lognormal(-1, 1);
  g ~ lognormal(-1, 1);
  beta ~ lognormal(-1, 1);
  tau ~ normal(1, 0.5);
  sigma ~ lognormal(-1, 1);

  
  for (k in 1:N) {
    y_init[k] ~ lognormal(log(x_init[k]), sigma); 
    y[ , k] ~ lognormal(log(x[, k]), sigma[k]);
  }
}


Looks like ts is an array but dx_dt expects a vector.

vector[N] z[N_TS] = ode_rk45(dx_dt, z_init, t0, ts, N, N_t, I, to_vector(ts), s, g, b, tau);

It still doesn’t compile but with this change the errors are unrelated to the ODE.

1 Like

With this change I get the error:

Compiling Stan program...
-
Semantic error in '/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmphjhKZn/model-e51f15233314.stan', line 120, column 22 to column 93:
   -------------------------------------------------
   118:    // states x are estimated as parameters z with some uncertainty. 
   119:    // these estimated parameters are used in the model description to relate them to measured data
   120:    vector[N] z[N_TS] = ode_rk45(dx_dt, z_init, t0, to_vector(ts), N, N_t, I, ts, s, g, b, tau);
                               ^
   121:  }
   122:  
   -------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45'. Expected arguments:
(real, vector) => vector, vector, real, real[]
Instead supplied arguments of incompatible type:
(real, vector, int, vector[], vector, vector, real, real, real, real) => vector, vector, real, vector, int, vector[], vector, real[], real, real, real, real

which still appears related to ode_rk45.

I had specified the type of ts based on https://mc-stan.org/docs/2_25/stan-users-guide/measurement-error-models.html so I’m not sure why this would be the problem.

A few things I noticed that are hopefully helpful:

Looking at your original post it seems like this should be x on the LHS and x_init on the RHS (instead of z and z_init). Or did I misunderstand?

Did you mean b instead of beta? It looks like you have a parameter b in the parameters block but no parameter beta. This will result in an error (once the ode error is fixed).

In the parameters block sigma is a scalar so indexing it here will give you an error (once the ode error goes away).

If I make all those changes then I can get this Stan program to compile successfully!

I’m sorry about the sloppy coding. I’ve been borrowing bits and pieces from other examples and have made a bit of a mess.

I made the corrections you suggested but unfortunately still got:

Compiling Stan program...
/
Semantic error in '/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmpDfdFlL/model-11b2c325c2e68.stan', line 120, column 22 to column 82:
   -------------------------------------------------
   118:    // states x are estimated as parameters z with some uncertainty. 
   119:    // these estimated parameters are used in the model description to relate them to measured data
   120:    vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
                               ^
   121:  }
   122:  
   -------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45'. Available signatures:
(real, vector, ...) => vector, vector, real, real[],  ...
.

Instead supplied arguments of incompatible type:
(real, vector, int, vector[], vector, vector, real, real, real, real) => vector, vector, real, real[], int, vector[], vector, real[], real, real, real, real.

The corrected stan file now looks like:


functions {
  
  real phi(real a){
    real out = (exp(2*a)-1)/(exp(2*a)+1);
    return out;
  }
  
  // From
  // https://discourse.mc-stan.org/t/zero-chain-movement-mixing-in-forced-ode-model/921
  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);
  }
  
  vector dx_dt(real t,       // time
              vector x,     // system state {1 x N - network activity for each time point}
              int N,
              vector[] N_t, //vector array
              vector I,
              vector ts,
              real s,
              real g,
              real b, 
              real tau) {
    
    vector[N] res;
    
    int d = find_interval_elem(t, ts, 1);
    
    for (i in 1:N){
      res[i] = (-x[i] + s * phi(x[i]) + g * N_t[d, i] + b * I[d])/tau;
    }
    
    // returns change for one time point for all nodes
    return  res;
  }
  
}

data {
  int<lower = 0> N_TS;           // number of measurement times
  int<lower = 0> N;             // number of nodes
  real t0;
  real ts[N_TS];                 // measurement times > 0
  vector[N] y_init;             // initial measured activity level

  // Arrays are declared by enclosing the dimensions in square brackets following the name of the variable.
  // this is an array of length N_TS consisting of (column) vectors of length N
  // Each array element is like a row
  vector[N] y[N_TS]; // measured activity level with nodes in cols and timepoints in rows
  matrix[N, N] W;// adjacency matrix

  vector[N_TS] I; // task stimulation
}

transformed data{
  vector[N] N_t[N_TS];

  for(i in 1:N){
    for(t in 1:N_TS){
      //row of incoming connections to node i of length n (W[i,] - 1xn)
      //multiplied with activity (column) vector of length n for all nodes at time t (y[t] - nx1)
      N_t[t, i] = dot_product(W[i,], y[t]);
    }
  }
}

parameters {
  real<lower = 0> s;   // self coupling
  real<lower = 0> g;  // global coupling
  real<lower = 0> b;  // task modulation
  real<lower = 0> tau; // time constant
  real<lower = 0> sigma;   // measurement error
  vector[N] x_init;
}


transformed parameters {
  // states x are estimated as parameters z with some uncertainty. 
  // these estimated parameters are used in the model description to relate them to measured data
  vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, ts, s, g, b, tau);
}

model {
  s ~ lognormal(-1, 1);
  g ~ lognormal(-1, 1);
  b ~ lognormal(-1, 1);
  tau ~ normal(1, 0.5);
  sigma ~ lognormal(-1, 1);


  for (k in 1:N) {
    y_init[k] ~ lognormal(log(x_init[k]), sigma);
    y[ , k] ~ lognormal(log(x[, k]), sigma);
  }
}

No worries, that can easily happen!

I think you also still need @nhuurre’s suggestion of using to_vector:

vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, to_vector(ts), s, g, b, tau);

Does that work?

Unfortunately not. I get a very similar but slightly different looking error

Compiling Stan program...
-
Semantic error in '/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmpDfdFlL/model-11b2c6b1a58d.stan', line 120, column 22 to column 93:
   -------------------------------------------------
   118:    // states x are estimated as parameters z with some uncertainty. 
   119:    // these estimated parameters are used in the model description to relate them to measured data
   120:    vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, to_vector(ts), N, N_t, I, ts, s, g, b, tau);
                               ^
   121:  }
   122:  
   -------------------------------------------------

Ill-typed arguments supplied to function 'ode_rk45'. Expected arguments:
(real, vector) => vector, vector, real, real[]
Instead supplied arguments of incompatible type:
(real, vector, int, vector[], vector, vector, real, real, real, real) => vector, vector, real, vector, int, vector[], vector, real[], real, real, real, real

Thank you for your patience with me.

No problem, thanks for your patience with Stan!

Which version of CmdStan do you have? You can check with cmdstanr::cmdstan_version(). I don’t get an error with CmdStan 2.25 and this Stan program:

functions {

  real phi(real a){
    real out = (exp(2*a)-1)/(exp(2*a)+1);
    return out;
  }

  // From
  // https://discourse.mc-stan.org/t/zero-chain-movement-mixing-in-forced-ode-model/921
  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);
  }

  vector dx_dt(real t,       // time
              vector x,     // system state {1 x N - network activity for each time point}
              int N,
              vector[] N_t, //vector array
              vector I,
              vector ts,
              real s,
              real g,
              real b,
              real tau) {

    vector[N] res;

    int d = find_interval_elem(t, ts, 1);

    for (i in 1:N){
      res[i] = (-x[i] + s * phi(x[i]) + g * N_t[d, i] + b * I[d])/tau;
    }

    // returns change for one time point for all nodes
    return  res;
  }

}

data {
  int<lower = 0> N_TS;           // number of measurement times
  int<lower = 0> N;             // number of nodes
  real t0;
  real ts[N_TS];                 // measurement times > 0
  vector[N] y_init;             // initial measured activity level

  // Arrays are declared by enclosing the dimensions in square brackets following the name of the variable.
  // this is an array of length N_TS consisting of (column) vectors of length N
  // Each array element is like a row
  vector[N] y[N_TS]; // measured activity level with nodes in cols and timepoints in rows
  matrix[N, N] W;// adjacency matrix

  vector[N_TS] I; // task stimulation
}

transformed data{
  vector[N] N_t[N_TS];

  for(i in 1:N){
    for(t in 1:N_TS){
      //row of incoming connections to node i of length n (W[i,] - 1xn)
      //multiplied with activity (column) vector of length n for all nodes at time t (y[t] - nx1)
      N_t[t, i] = dot_product(W[i,], y[t]);
    }
  }
}

parameters {
  real<lower = 0> s;   // self coupling
  real<lower = 0> g;  // global coupling
  real<lower = 0> b;  // task modulation
  real<lower = 0> tau; // time constant
  real<lower = 0> sigma;   // measurement error
  vector[N] x_init;
}


transformed parameters {
  // states x are estimated as parameters z with some uncertainty.
  // these estimated parameters are used in the model description to relate them to measured data
  vector[N] x[N_TS] = ode_rk45(dx_dt, x_init, t0, ts, N, N_t, I, to_vector(ts), s, g, b, tau);
}

model {
  s ~ lognormal(-1, 1);
  g ~ lognormal(-1, 1);
  b ~ lognormal(-1, 1);
  tau ~ normal(1, 0.5);
  sigma ~ lognormal(-1, 1);


  for (k in 1:N) {
    y_init[k] ~ lognormal(log(x_init[k]), sigma);
    y[ , k] ~ lognormal(log(x[, k]), sigma);
  }
}

I have 2.25.0 as well. I restarted my session, copied the Stan program you posted and ran

mod = cmdstanr::cmdstan_model("nmm_ode_forum.stan")

yet this time I got:

Compiling Stan program...
error: PCH file built from a different branch ((clang-1200.0.32.21)) than the compiler ((clang-1200.0.32.27))
1 error generated.
make: *** [/var/folders/10/zpzbg__x3hlc56gwq9zf6r600000gn/T/RtmpejWVOO/model-11d61773c7560] Error 1
Error: An error occured during compilation! See the message above for more information.

Oh that’s annoying. That’s a totally unrelated error. I think/hope

cmdstanr::rebuild_cmdstan()

will fix that PCH file error.

To be clear, I meant annoying for you the Stan user, not that your post was annoying!

No worries. I wasn’t offended ;) And more importantly: That worked! Thank you again for the suggestions and patience.

Ok great, and happy to help.

Do you mean just the “PCH file” error or that plus the Stan program also compiled successfully?

Both the PCH file error has disappeared and the Stan program has compiled.

I have other issues with the model but those aren’t relevant for this topic.

Ok great. At least it compiles now and I hope you’re able to sort out the other issues. Definitely feel free to start new topics on the forum if you get stuck. There are a bunch of people here who are good at troubleshooting ODE models (although I’m not one of them).