Update deprecated syntax for future rstan compatibility

DESCRIPTION

@@ -39,7 +39,7 @@ Imports:
     methods,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.18.1),
+    rstan (>= 2.26.0),
     rstantools (>= 2.1.1),
     data.table,
     tidyr,
@@ -55,8 +55,8 @@ LinkingTo:
     Rcpp (>= 0.12.0),
     RcppEigen (>= 0.3.3.3.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.18.1),
-    StanHeaders (>= 2.18.0)
+    rstan (>= 2.26.0),
+    StanHeaders (>= 2.26.0)
 SystemRequirements: GNU make
 URL: https://github.com/bgoussen/BeeGUTS
 BugReports: https://github.com/bgoussen/BeeGUTS/issues

inst/stan/GUTS_IT.stan

@@ -7,17 +7,17 @@ functions {
 
 #include /include/common_stan_functions.stan
 
-  real[] TKTD_varIT( real t,
-                     real[] y,
-                     real[] theta,
-                     real[] x_r,
-                     int[]  x_i) {
+  vector TKTD_varIT( real t,
+                     vector y,
+                     array[] real theta,
+                     array[] real x_r,
+                     array[] int  x_i) {
 
     // - parameters
     real kd = theta[1];
 
     // - new variables
-    real dy_dt[1]; //
+    vector[1] dy_dt; //
 
       // - latent variables
     int Nconc = x_i[1]; // Number of point measuring concentration
@@ -36,18 +36,21 @@ functions {
     return(dy_dt);
   }
 
-  matrix solve_TKTD_varIT(real[] y0, real t0, real[] ts, real[] theta, data real[] tconc, data real[] conc, data real[] odeParam){
+  matrix solve_TKTD_varIT(array[] real y0, real t0, array[] real ts, array[] real theta, data array[] real tconc, data array[] real conc, data real relTol, data real absTol, int maxSteps){
 
-    int x_i[1];
+    array[1] int x_i;
     x_i[1] = size(tconc);
-
-    return(to_matrix(
-      integrate_ode_rk45(TKTD_varIT, y0, t0, ts, theta,
-                         to_array_1d(append_row(to_vector(tconc), to_vector(conc))),
-                         x_i,
+    array[size(ts)] vector[1] ode_res
+      = ode_rk45_tol(TKTD_varIT, to_vector(y0), t0, ts,
                          // additional control parameters for the solver: real rel_tol, real abs_tol, int max_num_steps
-                         odeParam[1], odeParam[2], odeParam[3])));
-
+                         relTol, absTol, maxSteps, theta,
+                         to_array_1d(append_row(to_vector(tconc), to_vector(conc))),
+                         x_i);
+    matrix[size(ts), 1] rtn;
+    for(i in 1:size(ts)) {
+      rtn[i] = transpose(ode_res[i]);
+    }
+    return rtn;
   }
 }
 
@@ -65,19 +68,13 @@ data {
 }
 transformed data{
 
-  real<lower=0> y0[1];
-  real odeParam[3];
+  array[1] real<lower=0> y0;
 
-  real tNsurv_ode[nData_Nsurv]; // time of Nbr survival to include in the ode !
-  real tconc_ode[nData_conc]; // time of Nbr survival to include in the ode !
+  array[nData_Nsurv] real tNsurv_ode; // time of Nbr survival to include in the ode !
+  array[nData_conc] real tconc_ode; // time of Nbr survival to include in the ode !
 
   y0[1] = 1e-20; // cannot start at 0 for log(0)
 
-  // Add odeSolveParameters
-  odeParam[1] = relTol;
-  odeParam[2] = absTol;
-  odeParam[3] = maxSteps;
-
   for(gr in 1:nGroup){
     tNsurv_ode[idS_lw[gr]:idS_up[gr]] = tNsurv[idS_lw[gr]:idS_up[gr]];
     tNsurv_ode[idS_lw[gr]] = tNsurv[idS_lw[gr]] + 1e-9 ; // to start ode integrator at 0
@@ -90,17 +87,17 @@ parameters {
 
   real beta_log10  ;
 
-  real sigma[2 + nDatasets] ;
+  array[2 + nDatasets] real sigma ;
 
 
 }
 transformed parameters{
 
   real kd_log10 = kdMean_log10 + kdSD_log10 * sigma[1] ;
   real mw_log10 = mwMean_log10 + mwSD_log10 * sigma[2] ;
-  real hb_log10[nDatasets];
+  array[nDatasets] real hb_log10;
 
-  real<lower=0> param[1]; //
+  array[1] real<lower=0> param; //
 
   matrix[nData_Nsurv, 1] y_hat;
   vector<lower=0, upper=1>[nData_Nsurv] Psurv_hat;
@@ -118,7 +115,7 @@ transformed parameters{
   for(gr in 1:nGroup){
      real hb  = 10^hb_log10[groupDataset[gr]]; // hb
     /* initial time must be less than t0 = 0, so we use a very small close small number -1e-9 */
-      y_hat[idS_lw[gr]:idS_up[gr], 1] = solve_TKTD_varIT(y0, 0, tNsurv_ode[idS_lw[gr]:idS_up[gr]], param, tconc_ode[idC_lw[gr]:idC_up[gr]], conc[idC_lw[gr]:idC_up[gr]], odeParam)[,1];
+      y_hat[idS_lw[gr]:idS_up[gr], 1] = solve_TKTD_varIT(y0, 0, tNsurv_ode[idS_lw[gr]:idS_up[gr]], param, tconc_ode[idC_lw[gr]:idC_up[gr]], conc[idC_lw[gr]:idC_up[gr]], relTol, absTol, maxSteps)[,1];
 
     for(i in idS_lw[gr]:idS_up[gr]){
 

inst/stan/GUTS_SD.stan

@@ -7,11 +7,11 @@ functions {
 
 #include /include/common_stan_functions.stan
 
-  real[] TKTD_varSD( real t,
-                     real[] y,
-                     real[] theta,
-                     real[] x_r,
-                     int[]  x_i) {
+  vector TKTD_varSD( real t,
+                     vector y,
+                     array[] real theta,
+                     array[] real x_r,
+                     array[] int  x_i) {
 
     // - parameters
     real kd = theta[1];
@@ -20,8 +20,8 @@ functions {
     real hb = theta[4];
 
     // - new variables
-    real max_zw[2]; //
-    real dy_dt[2]; //
+    array[2] real max_zw; //
+    vector[2] dy_dt; //
 
     // - latent variables
     int Nconc = x_i[1]; // Number of point measuring concentration
@@ -46,17 +46,22 @@ functions {
     return(dy_dt);
   }
 
-  matrix solve_TKTD_varSD(real[] y0, real t0, real[] ts, real[] theta, data real[] tconc, data real[] conc, data real[] odeParam){
+  matrix solve_TKTD_varSD(array[] real y0, real t0, array[] real ts, array[] real theta, data array[] real tconc, data array[] real conc, data real relTol, data real absTol, int maxSteps){
 
-    int x_i[1];
+    array[1] int x_i;
     x_i[1] = size(tconc);
 
-    return(to_matrix(
-      integrate_ode_rk45(TKTD_varSD, y0, t0, ts, theta,
-                         to_array_1d(append_row(to_vector(tconc), to_vector(conc))),
-                         x_i,
+    array[size(ts)] vector[2] ode_res
+      = ode_rk45_tol(TKTD_varSD, to_vector(y0), t0, ts,
                          // additional control parameters for the solver: real rel_tol, real abs_tol, int max_num_steps
-                          odeParam[1], odeParam[2], odeParam[3])));
+                         relTol, absTol, maxSteps, theta,
+                         to_array_1d(append_row(to_vector(tconc), to_vector(conc))),
+                         x_i);
+    matrix[size(ts), 2] rtn;
+    for(i in 1:size(ts)) {
+      rtn[i] = transpose(ode_res[i]);
+    }
+    return rtn;
   }
 }
 
@@ -73,20 +78,14 @@ data {
 }
 transformed data{
 
-  real<lower=0> y0[2];
-  real odeParam[3];
+  array[2] real<lower=0> y0;
 
-  real tNsurv_ode[nData_Nsurv]; // time of Nbr survival to include in the ode !
-  real tconc_ode[nData_conc]; // time of Nbr survival to include in the ode !
+  array[nData_Nsurv] real tNsurv_ode; // time of Nbr survival to include in the ode !
+  array[nData_conc] real tconc_ode; // time of Nbr survival to include in the ode !
 
   y0[1] = 0;
   y0[2] = 0;
 
-  // Add odeSolveParameters
-  odeParam[1] = relTol;
-  odeParam[2] = absTol;
-  odeParam[3] = maxSteps;
-
   for(gr in 1:nGroup){
     tNsurv_ode[idS_lw[gr]:idS_up[gr]] = tNsurv[idS_lw[gr]:idS_up[gr]];
     tNsurv_ode[idS_lw[gr]] = tNsurv[idS_lw[gr]] + 1e-9 ; // to start ode integrator at 0
@@ -98,17 +97,17 @@ transformed data{
 }
 parameters {
 
-  real sigma[3 + nDatasets];
+  array[3 + nDatasets] real sigma;
 
 }
 transformed parameters{
 
   real kd_log10 = kdMean_log10 + kdSD_log10 * sigma[1];
   real zw_log10 = zwMean_log10 + zwSD_log10 * sigma[2];
   real bw_log10 = bwMean_log10 + bwSD_log10 * sigma[3];
-  real hb_log10[nDatasets];
+  array[nDatasets] real hb_log10;
 
-  real<lower=0> param[4]; //
+  array[4] real<lower=0> param; //
 
   matrix[nData_Nsurv,2] y_hat;
   vector<lower=0, upper=1>[nData_Nsurv] Psurv_hat;
@@ -128,7 +127,7 @@ transformed parameters{
     param[4] = 10^hb_log10[groupDataset[gr]]; // hb
 
   /* initial time must be less than t0 = 0, so we use a very small close small number 1e-9 to at at time tNsurv and tconc */
-    y_hat[idS_lw[gr]:idS_up[gr],1:2] = solve_TKTD_varSD(y0, 0, tNsurv_ode[idS_lw[gr]:idS_up[gr]], param, tconc_ode[idC_lw[gr]:idC_up[gr]], conc[idC_lw[gr]:idC_up[gr]], odeParam);
+    y_hat[idS_lw[gr]:idS_up[gr],1:2] = solve_TKTD_varSD(y0, 0, tNsurv_ode[idS_lw[gr]:idS_up[gr]], param, tconc_ode[idC_lw[gr]:idC_up[gr]], conc[idC_lw[gr]:idC_up[gr]], relTol, absTol, maxSteps);
 
     Psurv_hat[idS_lw[gr]:idS_up[gr]] = exp( - y_hat[idS_lw[gr]:idS_up[gr], 2]);
 

inst/stan/include/common_stan_functions.stan

@@ -33,7 +33,7 @@ int find_interval_elem(real x, vector sorted, int start_ind){
         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_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; }

inst/stan/include/data_guts.stan

@@ -6,24 +6,24 @@ int <lower=1> nDatasets;
 
 // Number of groups
 int<lower=1> nGroup; // Number of groups (one group is combination of one dataset and one treatment)
-int groupDataset[nGroup]; // Corresponding dataset for each group
+array[nGroup] int groupDataset; // Corresponding dataset for each group
 
 // Concentration
 int<lower=1> nData_conc; // length of data for concentration
-real conc[nData_conc]; // concentration
-real tconc[nData_conc]; // time of concentration
+array[nData_conc] real conc; // concentration
+array[nData_conc] real tconc; // time of concentration
 
-int<lower=1> idC_lw[nGroup]; // e.g. 1 6 12 18
-int<lower=1> idC_up[nGroup]; // e.g. 6 12 18 24
+array[nGroup] int<lower=1> idC_lw; // e.g. 1 6 12 18
+array[nGroup] int<lower=1> idC_up; // e.g. 6 12 18 24
 
 // Survivors
 int<lower=1> nData_Nsurv; // number of group: 4
-int Nsurv[nData_Nsurv];
-int Nprec[nData_Nsurv];
-real tNsurv[nData_Nsurv]; // time of Nbr survival
+array[nData_Nsurv] int Nsurv;
+array[nData_Nsurv] int Nprec;
+array[nData_Nsurv] real tNsurv; // time of Nbr survival
 
-int<lower=1> idS_lw[nGroup]; // e.g. 1 6 12 18
-int<lower=1> idS_up[nGroup]; // e.g. 6 12 18 24
+array[nGroup] int<lower=1> idS_lw; // e.g. 1 6 12 18
+array[nGroup] int<lower=1> idS_up; // e.g. 6 12 18 24
 
 // PRIORS
 real hbMean_log10;

inst/stan/include/gen_quantities_guts.stan

@@ -1,9 +1,9 @@
 /* Code adapted from Virgile Baudrot
 https://github.com/virgile-baudrot/gutsRstan */
 
-int Nsurv_ppc[nData_Nsurv];
-int Nsurv_sim[nData_Nsurv];
-int Nsurv_sim_prec[nData_Nsurv];
+array[nData_Nsurv] int Nsurv_ppc;
+array[nData_Nsurv] int Nsurv_sim;
+array[nData_Nsurv] int Nsurv_sim_prec;
 
 vector[nData_Nsurv] log_lik;