Update deprecated syntax for future rstan compatibility

DESCRIPTION

@@ -25,7 +25,7 @@ Imports:
     survival,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.19.2),
+    rstan (>= 2.26.0),
     rstantools (>= 2.0.0),
     Rdpack (>= 0.7),
     tibble (>= 2.1.3),
@@ -50,8 +50,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
 RoxygenNote: 7.2.3
 Roxygen: list(markdown = TRUE)

inst/stan/binomial_1par.stan

@@ -6,9 +6,9 @@ data {
 #include /include/data_common.stan
 
   // Outcomes
-  int<lower=0, upper=1> ipd_r[ni_ipd];
-  int<lower=0> agd_arm_n[ni_agd_arm];
-  int<lower=0> agd_arm_r[ni_agd_arm];
+  array[ni_ipd] int<lower=0, upper=1> ipd_r;
+  array[ni_agd_arm] int<lower=0> agd_arm_n;
+  array[ni_agd_arm] int<lower=0> agd_arm_r;
 }
 transformed data {
 #include /include/transformed_data_common.stan

inst/stan/binomial_2par.stan

@@ -6,9 +6,9 @@ data {
 #include /include/data_common.stan
 
   // Outcomes
-  int<lower=0, upper=1> ipd_r[ni_ipd];
-  int<lower=0> agd_arm_n[ni_agd_arm];
-  int<lower=0> agd_arm_r[ni_agd_arm];
+  array[ni_ipd] int<lower=0, upper=1> ipd_r;
+  array[ni_agd_arm] int<lower=0> agd_arm_n;
+  array[ni_agd_arm] int<lower=0> agd_arm_r;
 }
 transformed data {
 #include /include/transformed_data_common.stan
@@ -18,8 +18,8 @@ parameters {
 }
 transformed parameters {
   vector[ni_agd_arm] theta2_agd_arm_bar;
-  real<lower=0> nprime[ni_agd_arm];
-  real<lower=0, upper=1> pprime[ni_agd_arm];
+  array[ni_agd_arm] real<lower=0> nprime;
+  array[ni_agd_arm] real<lower=0, upper=1> pprime;
 #include /include/transformed_parameters_theta.stan
 #include /include/transformed_parameters_common.stan
 

inst/stan/include/data_common.stan

@@ -10,25 +10,25 @@ int<lower=0> ni_agd_contrast; // total number of AgD (contrast-based) data point
 
 // Treatment IDs
 int<lower=0> narm_ipd; // Number of IPD arms
-int<lower=1> ipd_arm[ni_ipd]; // Arm indicator for IPD (i.e. picking element of which_RE)
-int<lower=1> ipd_trt[narm_ipd];
+array[ni_ipd] int<lower=1> ipd_arm; // Arm indicator for IPD (i.e. picking element of which_RE)
+array[narm_ipd] int<lower=1> ipd_trt;
 int<lower=0> narm_agd_arm;
-int<lower=1> agd_arm_trt[narm_agd_arm];
-int<lower=1> agd_contrast_trt[ni_agd_contrast];
-int<lower=1> agd_contrast_trt_b[ni_agd_contrast];
+array[narm_agd_arm] int<lower=1> agd_arm_trt;
+array[ni_agd_contrast] int<lower=1> agd_contrast_trt;
+array[ni_agd_contrast] int<lower=1> agd_contrast_trt_b;
 
 // Study IDs
-// int<lower=1> ipd_study[max(ipd_arm)];
-// int<lower=1> agd_arm_study[ni_agd_arm];
-// int<lower=1> agd_contrast_study[ni_agd_contrast];
+// array[max(ipd_arm)] int<lower=1> ipd_study;
+// array[ni_agd_arm] int<lower=1> agd_arm_study;
+// array[ni_agd_contrast] int<lower=1> agd_contrast_study;
 
 int<lower=1> nt; // number of treatments
 int<lower=0> nX; // number of columns of design matrix
 
 // Integration
 int<lower=1> nchains;
 int<lower=1,upper=nchains> CHAIN_ID;
-int<lower=1> nint_vec[nchains]; // number of samples for numerical integration (1 = no integration)
+array[nchains] int<lower=1> nint_vec; // number of samples for numerical integration (1 = no integration)
 int<lower=1> nint_max;
 int<lower=0> int_thin; // thinning factor for saved integration points
 
@@ -49,7 +49,7 @@ vector[has_offset ? ni_ipd + nint_max * (ni_agd_arm + ni_agd_contrast) : 0] offs
 
 // -- Random effects --
 int<lower=0, upper=1> RE; // Random effects flag (yes = 1)
-int<lower=0> which_RE[RE ? narm_ipd + narm_agd_arm + ni_agd_contrast : 0]; // ID of RE delta for each arm (0 for no RE delta)
+array[RE ? narm_ipd + narm_agd_arm + ni_agd_contrast : 0] int<lower=0> which_RE; // ID of RE delta for each arm (0 for no RE delta)
 corr_matrix[RE ? max(which_RE) : 1] RE_cor; // RE correlation matrix
 
 // -- Node-splitting --

inst/stan/include/transformed_data_common.stan

@@ -14,12 +14,12 @@ matrix[0, 0] REdummy; // Use dummy zero-dim matrix to set RE_L to [0x0] when n_d
 cholesky_factor_corr[n_delta] RE_L = n_delta ? cholesky_decompose(RE_cor) : REdummy;
 // Sparse representation
 vector[0] wdummy;
-int vudummy[0];
+array[0] int vudummy;
 int RE_L_nz = count_nonzero(RE_L); // Number of non-zero entries
 int RE_sparse = RE_L_nz * 1.0 / num_elements(RE_L) <= 0.1; // Use sparse representation? (yes = 1)
 vector[RE_sparse ? RE_L_nz : 0] RE_L_w = RE_sparse ? csr_extract_w(RE_L): wdummy; // Non-zero entries
-int RE_L_v[RE_sparse ? RE_L_nz : 0] = RE_sparse ? csr_extract_v(RE_L): vudummy; // v sparse component
-int RE_L_u[RE_sparse ? n_delta + 1 : 0] = RE_sparse ? csr_extract_u(RE_L) : vudummy; // u sparse component
+array[RE_sparse ? RE_L_nz : 0] int RE_L_v = RE_sparse ? csr_extract_v(RE_L): vudummy; // v sparse component
+array[RE_sparse ? n_delta + 1 : 0] int RE_L_u = RE_sparse ? csr_extract_u(RE_L) : vudummy; // u sparse component
 
 // Total number of data points
 // int totni = ni_ipd + nint * (ni_agd_arm + ni_agd_contrast);
@@ -31,7 +31,7 @@ int totns = ns_ipd + ns_agd_arm; // + ns_agd_contrast;
 // int<lower=0> narm_ipd = ni_ipd ? max(ipd_arm) : 0;
 
 // All treatments vector
-int<lower=1> trt[narm_ipd + narm_agd_arm + ni_agd_contrast] = append_array(append_array(ipd_trt, agd_arm_trt), agd_contrast_trt);
+array[narm_ipd + narm_agd_arm + ni_agd_contrast] int<lower=1> trt = append_array(append_array(ipd_trt, agd_arm_trt), agd_contrast_trt);
 
 // Split Q matrix or X matrix into IPD and AgD rows
 matrix[0, nX] Xdummy;
@@ -46,14 +46,14 @@ vector[has_offset && ni_agd_arm ? nint_max * ni_agd_arm : 0] offset_agd_arm = ha
 vector[has_offset && ni_agd_contrast ? nint_max * ni_agd_contrast : 0] offset_agd_contrast = has_offset && ni_agd_contrast ? offsets[(ni_ipd + nint_max * ni_agd_arm + 1):(ni_ipd + nint_max * (ni_agd_arm + ni_agd_contrast))] : odummy;
 
 // nint/int_thin for numerical integration checks
-int n_int_thin = (nint_max > 1 && int_thin > 0) ? nint / int_thin : 0;
+int n_int_thin = (nint_max > 1 && int_thin > 0) ? nint %/% int_thin : 0;
 
 // Inverse covariance matrix for contrasts
 matrix[ni_agd_contrast ? ni_agd_contrast : 1, ni_agd_contrast ? ni_agd_contrast : 1] inv_Sigma = inverse_spd(agd_contrast_Sigma);
 
 // Construct number of contrasts in each study for contrast-based AgD by looking at Sigma covariance matrix
 // NOTE: Sigma must be block diagonal (i.e. all contrasts for a single study together)
-int nc_agd_contrast[ns_agd_contrast];
+array[ns_agd_contrast] int nc_agd_contrast;
 if (ns_agd_contrast) {
   int s = 1;
   int c = 1;

inst/stan/include/vector_functions.stan

@@ -1,9 +1,9 @@
 // -- Vector functions --
 
 // Which entries are equal to a given value
-int[] which(int[] x, int y) {
+array[] int which(array[] int x, int y) {
   int n = num_elements(x);
-  int w[n]; // Over-allocate w and then truncate later
+  array[n] int w; // Over-allocate w and then truncate later
   int c = 1;
   for (i in 1:n) {
     if (x[i] == y) {
@@ -15,9 +15,9 @@ int[] which(int[] x, int y) {
 }
 
 // Which entries are greater than 0
-int[] which_gt0(vector x) {
+array[] int which_gt0(vector x) {
   int n = num_elements(x);
-  int w[n]; // Over-allocate w and then truncate later
+  array[n] int w; // Over-allocate w and then truncate later
   int c = 1;
   for (i in 1:n) {
     if (x[i] > 0) {

inst/stan/normal.stan

@@ -12,9 +12,9 @@ data {
   real<lower=0> prior_aux_df;
 
   // Outcomes
-  real ipd_y[ni_ipd];
-  real agd_arm_y[ni_agd_arm];
-  real<lower=0> agd_arm_se[ni_agd_arm];
+  array[ni_ipd] real ipd_y;
+  array[ni_agd_arm] real agd_arm_y;
+  array[ni_agd_arm] real<lower=0> agd_arm_se;
 }
 transformed data {
 #include /include/transformed_data_common.stan

inst/stan/ordered_multinomial.stan

@@ -8,14 +8,14 @@ data {
   // Outcomes
   int<lower=2> ncat;
 
-  int<lower=1, upper=ncat> ipd_r[ni_ipd];
-  int<lower=0, upper=ncat> ipd_cat[ni_ipd, ncat];  // Category specs, left-aligned, padded with zeros
-  int<lower=2, upper=ncat> ipd_ncat[ni_ipd];  // Number of categories observed
+  array[ni_ipd] int<lower=1, upper=ncat> ipd_r;
+  array[ni_ipd, ncat] int<lower=0, upper=ncat> ipd_cat;  // Category specs, left-aligned, padded with zeros
+  array[ni_ipd] int<lower=2, upper=ncat> ipd_ncat;  // Number of categories observed
 
-  int<lower=0> agd_arm_r[ni_agd_arm, ncat];
+  array[ni_agd_arm, ncat] int<lower=0> agd_arm_r;
   vector[ni_agd_arm] agd_arm_n; // AgD arm sample sizes
-  int<lower=0, upper=ncat> agd_arm_cat[ni_agd_arm, ncat];  // Category specs, left-aligned, padded with zeros
-  int<lower=2, upper=ncat> agd_arm_ncat[ni_agd_arm];  // Number of categories observed
+  array[ni_agd_arm, ncat] int<lower=0, upper=ncat> agd_arm_cat;  // Category specs, left-aligned, padded with zeros
+  array[ni_agd_arm] int<lower=2, upper=ncat> agd_arm_ncat;  // Number of categories observed
 
   // Prior on differences between cutpoints
   int<lower=0,upper=6> prior_aux_dist;
@@ -25,8 +25,8 @@ data {
 
 }
 transformed data {
-  vector[ncat] theta_ipd0[ni_ipd] = rep_array(rep_vector(0, ncat), ni_ipd);
-  vector[ncat] theta_agd_arm_bar0[ni_agd_arm] = rep_array(rep_vector(0, ncat), ni_agd_arm);
+  array[ni_ipd] vector[ncat] theta_ipd0 = rep_array(rep_vector(0, ncat), ni_ipd);
+  array[ni_agd_arm] vector[ncat] theta_agd_arm_bar0 = rep_array(rep_vector(0, ncat), ni_agd_arm);
   // matrix[ni_agd_arm * n_int_thin, ncat] theta_bar_cum_agd_arm0 = rep_matrix(0, ni_agd_arm * n_int_thin, ncat);
 #include /include/transformed_data_common.stan
 }
@@ -40,11 +40,11 @@ parameters {
 transformed parameters {
   vector[ncat - 1] cc;
 
-  vector[ncat] theta_ipd[ni_ipd]; // IPD transformed predictor
+  array[ni_ipd] vector[ncat] theta_ipd; // IPD transformed predictor
 
   matrix[nint_max > 1 ? nint * ni_agd_arm : 0, nint_max > 1 ? ncat - 1 : 0] theta_agd_arm_ii; // Use these as q_ii intermediates
-  vector[ncat] q_agd_arm_bar[ni_agd_arm]; // AgD arm transformed predictor
-  vector[ncat] theta_agd_arm_bar[ni_agd_arm]; // AgD arm transformed predictor
+  array[ni_agd_arm] vector[ncat] q_agd_arm_bar; // AgD arm transformed predictor
+  array[ni_agd_arm] vector[ncat] theta_agd_arm_bar; // AgD arm transformed predictor
 
 #include /include/transformed_parameters_common.stan
 
@@ -271,8 +271,8 @@ model {
 }
 generated quantities {
   // Note: fitted values and theta_bar_cum_agd_arm will be 0 for missing categories
-  vector[ncat] fitted_ipd[ni_ipd];
-  vector[ncat] fitted_agd_arm[ni_agd_arm];
+  array[ni_ipd] vector[ncat] fitted_ipd;
+  array[ni_agd_arm] vector[ncat] fitted_agd_arm;
   matrix[ni_agd_arm * n_int_thin, ncat] theta_bar_cum_agd_arm = rep_matrix(0, ni_agd_arm * n_int_thin, ncat);
 
 #include /include/generated_quantities_common.stan
@@ -286,7 +286,7 @@ generated quantities {
 
   // AgD (arm-based) log likelihood and residual deviance
   {
-    vector[ncat] dv[ni_agd_arm];
+    array[ni_agd_arm] vector[ncat] dv;
     for (i in 1:ni_agd_arm) {
       log_lik[ni_ipd + i] = multinomial_lpmf(agd_arm_r[i] | theta_agd_arm_bar[i]);
       fitted_agd_arm[i] = agd_arm_n[i] * theta_agd_arm_bar[i];

inst/stan/poisson.stan

@@ -6,9 +6,9 @@ data {
 #include /include/data_common.stan
 
   // Outcomes
-  int<lower=0> ipd_r[ni_ipd];
+  array[ni_ipd] int<lower=0> ipd_r;
   vector<lower=0>[ni_ipd] ipd_E;
-  int<lower=0> agd_arm_r[ni_agd_arm];
+  array[ni_agd_arm] int<lower=0> agd_arm_r;
   vector<lower=0>[ni_agd_arm] agd_arm_E;
 }
 transformed data {

inst/stan/survival_mspline.stan

@@ -102,26 +102,26 @@ data {
 
   // AgD arm-based
   // int<lower=0> narm_agd_arm; // Number of arm-based AgD arms
-  int<lower=1> agd_arm_arm[ni_agd_arm]; // Arm indicator for AgD (arm-based) (i.e. picking element of which_RE)
+  array[ni_agd_arm] int<lower=1> agd_arm_arm; // Arm indicator for AgD (arm-based) (i.e. picking element of which_RE)
 
   // Outcomes
-  row_vector[n_scoef] ipd_time[ni_ipd];
-  row_vector[n_scoef] ipd_itime[ni_ipd];
-  row_vector[n_scoef] ipd_start_itime[ni_ipd];
-  row_vector[n_scoef] ipd_delay_itime[ni_ipd];
-  int<lower=0, upper=1> ipd_delayed[ni_ipd];
-  int<lower=0, upper=3> ipd_status[ni_ipd];
-
-  row_vector[n_scoef] agd_arm_time[ni_agd_arm];
-  row_vector[n_scoef] agd_arm_itime[ni_agd_arm];
-  row_vector[n_scoef] agd_arm_start_itime[ni_agd_arm];
-  row_vector[n_scoef] agd_arm_delay_itime[ni_agd_arm];
-  int<lower=0, upper=1> agd_arm_delayed[ni_agd_arm];
-  int<lower=0, upper=3> agd_arm_status[ni_agd_arm];
+  array[ni_ipd] row_vector[n_scoef] ipd_time;
+  array[ni_ipd] row_vector[n_scoef] ipd_itime;
+  array[ni_ipd] row_vector[n_scoef] ipd_start_itime;
+  array[ni_ipd] row_vector[n_scoef] ipd_delay_itime;
+  array[ni_ipd] int<lower=0, upper=1> ipd_delayed;
+  array[ni_ipd] int<lower=0, upper=3> ipd_status;
+
+  array[ni_agd_arm] row_vector[n_scoef] agd_arm_time;
+  array[ni_agd_arm] row_vector[n_scoef] agd_arm_itime;
+  array[ni_agd_arm] row_vector[n_scoef] agd_arm_start_itime;
+  array[ni_agd_arm] row_vector[n_scoef] agd_arm_delay_itime;
+  array[ni_agd_arm] int<lower=0, upper=1> agd_arm_delayed;
+  array[ni_agd_arm] int<lower=0, upper=3> agd_arm_status;
 
   // aux IDs for independent spline coefficients parameters
   int<lower=0, upper=1> aux_by; // Flag subgroup aux parameters within each arm (1 = yes)
-  int<lower=1> aux_id[ni_ipd + ni_agd_arm*(aux_by ? nint_max : 1)];
+  array[ni_ipd + ni_agd_arm*(aux_by ? nint_max : 1)] int<lower=1> aux_id;
 }
 transformed data {
   // Dirichlet prior vector
@@ -136,7 +136,7 @@ parameters {
 #include /include/parameters_common.stan
 
   // Spline coefficients
-  simplex[n_scoef] scoef[n_aux];
+  array[n_aux] simplex[n_scoef] scoef;
 }
 transformed parameters {
   // Log likelihood contributions

inst/stan/survival_param.stan

@@ -152,7 +152,7 @@ functions {
   }
 
   // -- Log likelihood with censoring and truncation --
-  vector loglik(int dist, vector time, vector start_time, vector delay_time, int[] status, vector eta, vector aux, vector aux2) {
+  vector loglik(int dist, vector time, vector start_time, vector delay_time, array[] int status, vector eta, vector aux, vector aux2) {
     int n = num_elements(eta);
     vector[n] l;
 
@@ -162,11 +162,11 @@ functions {
     int nwhich2 = num_elements(which(status, 2));
     int nwhich3 = num_elements(which(status, 3));
     int nwhichd = num_elements(which_gt0(delay_time));
-    int which0[nwhich0];
-    int which1[nwhich1];
-    int which2[nwhich2];
-    int which3[nwhich3];
-    int whichd[nwhichd];
+    array[nwhich0] int which0;
+    array[nwhich1] int which1;
+    array[nwhich2] int which2;
+    array[nwhich3] int which3;
+    array[nwhichd] int whichd;
     if (nwhich0) which0 = which(status, 0);
     if (nwhich1) which1 = which(status, 1);
     if (nwhich2) which2 = which(status, 2);
@@ -267,22 +267,22 @@ data {
 
   // AgD arm-based
   // int<lower=0> narm_agd_arm; // Number of arm-based AgD arms
-  int<lower=1> agd_arm_arm[ni_agd_arm]; // Arm indicator for AgD (arm-based) (i.e. picking element of which_RE)
+  array[ni_agd_arm] int<lower=1> agd_arm_arm; // Arm indicator for AgD (arm-based) (i.e. picking element of which_RE)
 
   // Outcomes
   vector[ni_ipd] ipd_time;
   vector[ni_ipd] ipd_start_time;
   vector[ni_ipd] ipd_delay_time;
-  int<lower=0, upper=3> ipd_status[ni_ipd];
+  array[ni_ipd] int<lower=0, upper=3> ipd_status;
 
   vector[ni_agd_arm] agd_arm_time;
   vector[ni_agd_arm] agd_arm_start_time;
   vector[ni_agd_arm] agd_arm_delay_time;
-  int<lower=0, upper=3> agd_arm_status[ni_agd_arm];
+  array[ni_agd_arm] int<lower=0, upper=3> agd_arm_status;
 
   // Aux IDs for independent shape parameters
   int<lower=0, upper=1> aux_by; // Flag subgroup aux parameters within each arm (1 = yes)
-  int<lower=1> aux_id[ni_ipd + ni_agd_arm * (aux_by ? nint_max : 1)];
+  array[ni_ipd + ni_agd_arm * (aux_by ? nint_max : 1)] int<lower=1> aux_id;
 }
 transformed data {
   // Exponential model indicator, 0 = exponential
@@ -292,7 +292,7 @@ transformed data {
   // Number of auxiliary parameters
   int n_aux = nonexp ? max(aux_id) : 0;
   // AgD arm indicator - shifted
-  int agd_arm_arm2[ni_agd_arm];
+  array[ni_agd_arm] int agd_arm_arm2;
   // Number of integration points
 
 #include /include/transformed_data_common.stan