R CMD SHLIB error for ODE-model

Hello,

I have written an IRT model which involves a user-defined probability function. To compute the latter, I need to integrate the exponential of a spline. I accomplished this via “integrate_ode_rk45” and it compiled and worked.

However, my first solution involved a lot of loops and I tried to refactor the code to use vectorized expressions instead of loops. Now the new model (see below) does compile when calling “stanc”. However, when I try to fit the model, I get an error message involving the following error message (translated from German):

"Execution of command ‘D:/R-4.0.2/bin/x64/R CMD SHLIB file3588724c109.cpp 2> file3588724c109.cpp.err.txt’ resulted in Status 1".

Via uncommenting and changing some statements in the underlying stan code, I was able to narrow the error down to the line, wherein I call the integrate_ode_rk45 function. However, I have no idea why this line causes the error.

Thank you for any advices/sugesstions on how to resolve this problem!

functions 
{
  vector build_b_spline(real[] t, real[] ext_knots, int ind, int order);
  vector build_b_spline(real[] t, real[] ext_knots, int ind, int order) 
 {
    vector[size(t)] b_spline;
    vector[size(t)] w1 = rep_vector(0, size(t));
    vector[size(t)] w2 = rep_vector(0, size(t));
    if (order==1)
      for (i in 1:size(t)) 
        b_spline[i] = (ext_knots[ind] <= t[i]) && (t[i] < ext_knots[ind+1]);
    else {
      if (ext_knots[ind] != ext_knots[ind+order-1])
        w1 = (to_vector(t) - rep_vector(ext_knots[ind], size(t))) /
             (ext_knots[ind+order-1] - ext_knots[ind]);
      if (ext_knots[ind+1] != ext_knots[ind+order])
        w2 = 1 - (to_vector(t) - rep_vector(ext_knots[ind+1], size(t))) /
                 (ext_knots[ind+order] - ext_knots[ind+1]);
      b_spline = w1 .* build_b_spline(t, ext_knots, ind, order-1) +
                 w2 .* build_b_spline(t, ext_knots, ind+1, order-1);
    }
    return b_spline;
 }

  real[] exp_spline_ode(real time, real[] state, real[] theta, real[] x_r, int[] x_i)
 {
  real dydt[1];
  real first[1];
  first[1] = time;
  dydt[1] = exp(theta[1]*build_b_spline(first, x_r, x_i[1], x_i[2])[1]);
  return dydt;
 }
  
  real espl_lpdf(vector y, real[] theta, real coeff, data real[] ext_knots) 
 {
    real loc[rows(y)];
    real deriv_spline[rows(y)];
    real ruck[rows(y)];
    int x_i[2];
    real times[rows(y)];
    real theta_h[1];
    real init_st[rows(y)] = rep_array(0.0, rows(y));
    //real init_st[1] = {0.0};
    x_i[1] = 1;
    x_i[2] = 2;
    times = to_array_1d(y);
    theta_h[1] = coeff;
    deriv_spline = to_array_1d(build_b_spline(to_array_1d(y), ext_knots, 1, 2));
    //This call is causing the trouble!
    loc = (integrate_ode_rk45(exp_spline_ode, init_st, ext_knots[1], times, theta_h, ext_knots, x_i))[,1]; 
    return sum(std_normal_lpdf(to_vector(loc)-to_vector(theta))+coeff*to_vector(deriv_spline));
 }
}

data {
  int<lower=0> K;         // number of items
  int<lower=1> N;         // number of test takers
  matrix[N, K] Y;         // responses of test takers to the items
  int num_knots;
  real ext_knots[num_knots];
}


parameters {
  vector[N] z;        
  vector[K] intercepts; 
  vector[K] coeffs;  
  real<lower=0> sigma_theta;    
}

transformed parameters {
  vector[N] theta = sigma_theta*z;
}


model { 
  sigma_theta ~ cauchy(0,1);
  intercepts ~ normal(0,1);
  coeffs ~ normal(0,1);
  z ~ normal(0, 1);
  for(j in 1:K)
 {
   Y[,j] ~ espl(to_array_1d(theta-intercepts[j]),coeffs[j],ext_knots); 
 }

}

Not sure whats going on…could you try integrate_1d?

Well, in fact, I started with the “integrate_1d” function - but this did not work. Hence, I switched to the ODE-solver which already worked in the “loop”-setting. But after my attempt to rewrite the model in a more efficient manner, the above error appears.