Merge pull request #16 from andrjohns/array-syntax

Update deprecated syntax for future rstan compatibility

DESCRIPTION

@@ -33,15 +33,15 @@ Imports:
     Matrix (>= 1.3),
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.18.1),
+    rstan (>= 2.26.0),
     rstantools (>= 2.1.1)    
 LinkingTo: 
     BH (>= 1.66.0),
     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)
 Suggests:
     testthat,
     knitr,

R/stan_car.R

@@ -340,7 +340,7 @@ stan_car <- function(formula,
         y_int = y_int,
         trials = rep(0, length(y)),
         n = n,
-        offset = offset,
+        input_offset = offset,
         has_re = has_re,
         n_ids = n_ids,
         id = id_index$idx,

R/stan_esf.R

@@ -292,7 +292,7 @@ stan_esf <- function(formula,
         y_int = y_int,
         trials = rep(0, length(y)),
         n = n,
-        offset = offset,
+        input_offset = offset,
         has_re = has_re,
         n_ids = n_ids,
         id = id_index$idx,

R/stan_glm.R

@@ -316,7 +316,7 @@ stan_glm <- function(formula,
     y_int = y_int,
     trials = rep(0, length(y)),
     n = n,
-    offset = offset,
+    input_offset = offset,
     has_re = has_re,
     n_ids = n_ids,
     id = id_index$idx,

R/stan_icar.R

@@ -348,7 +348,7 @@ stan_icar <- function(formula,
         y_int = y_int,
         trials = rep(0, length(y)),
         n = n,
-        offset = offset,
+        input_offset = offset,
         has_re = has_re,
         n_ids = n_ids,
         id = id_index$idx,

R/stan_sar.R

@@ -333,7 +333,7 @@ stan_sar <- function(formula,
         y_int = y_int,
         trials = rep(0, length(y)),
         n = n,
-        offset = offset,
+        input_offset = offset,
         has_re = has_re,
         n_ids = n_ids,
         id = id_index$idx,

inst/stan/parts/data.stan

@@ -9,42 +9,42 @@
 // censored counts
   int<lower=0> n_mis;
   int<lower=0> n_obs;
-  int y_mis_idx[n_mis];
-  int y_obs_idx[n_obs];
+  array[n_mis] int y_mis_idx;
+  array[n_obs] int y_obs_idx;
   int censor_point;
 
 // outcome 
   vector[n] y;
-  int<lower=0> y_int[n];
-  int<lower=0> trials[n];
+  array[n] int<lower=0> y_int;
+  array[n] int<lower=0> trials;
   int<lower=0,upper=1> prior_only;
 
 // offset
-  vector[n] offset; 
+  vector[n] input_offset; 
 
 // connectivity matrix: row-standardized for spatial lag of X 
   int<lower=0> dwx;
-  int wx_idx[dwx ? dwx : 1];
+  array[dwx ? dwx : 1] int wx_idx;
   int<lower=0> dw_nonzero;
   vector[dw_nonzero] W_w;
-  int W_v[dw_nonzero];
-  int W_u[dwx ? n + 1 : 1];
+  array[dw_nonzero] int W_v;
+  array[dwx ? n + 1 : 1] int W_u;
 
 // covariates and observational error information
     // lower, upper bounds for bounded data models
-  real bounds[2];
+  array[2] real bounds;
     // no. columns
   int<lower=0> dx_obs;
   int<lower=0> dx_me;
-  int<lower=0,upper=1> use_logit[dx_me];
+  array[dx_me] int<lower=0,upper=1> use_logit;
     // indices matching columns of observed and ME data matrices to columns of raw data matrix x (and parameter x_all)
-  int<lower=0> x_obs_idx[dx_obs ? dx_obs : 1]; 
-  int<lower=0> x_me_idx[dx_me ? dx_me : 1];  
+  array[dx_obs ? dx_obs : 1] int<lower=0> x_obs_idx; 
+  array[dx_me ? dx_me : 1] int<lower=0> x_me_idx;  
     // covariates observed with practical certainty 
   matrix[n, dx_obs ? dx_obs : 0] x_obs;
     // covariates observed with uncertainty, and standard errors
-  vector[n] x_me[dx_me];
-  vector<lower=0>[n] sigma_me[dx_me];
+  array[dx_me] vector[n] x_me;
+  array[dx_me] vector<lower=0>[n] sigma_me;
   //  priors for x_true
   vector[dx_me] prior_nux_true_alpha;
   vector[dx_me] prior_nux_true_beta;
@@ -60,17 +60,17 @@
   int nAx_w;
   int nC;
   vector[nAx_w] Ax_w;
-  int Ax_v[nAx_w];
-  int Ax_u[n + 1];
-  int Cidx[nC];
+  array[nAx_w] int Ax_v;
+  array[n + 1] int Ax_u;
+  array[nC] int Cidx;
   vector[n] Delta_inv;
   real log_det_Delta_inv;
   vector[n] lambda;
 
 // non-spatial partial pooling 
   int<lower=0,upper=1> has_re; // has varying intercept?
   int<lower=0> n_ids; // number of units
-  int<lower=0,upper=n_ids> id[n]; // identifier for the observational units associated with the varying intercepts
+  array[n] int<lower=0,upper=n_ids> id; // identifier for the observational units associated with the varying intercepts
 // priors
   vector[2] prior_alpha; // prior on the intercept
   int<lower=0> dbeta_prior;
@@ -84,15 +84,15 @@
   // ICAR 
   int<lower=0,upper=3> type; // 0=glm, 1=icar, 2=bym, 3=bym2
   int<lower=1> k; // no. of groups
-  int group_size[k]; // observational units per group
-  int group_idx[n]; // index of observations, ordered by group
+  array[k] int group_size; // observational units per group
+  array[n] int group_idx; // index of observations, ordered by group
   int<lower=0> m; // no of components requiring additional intercepts
   matrix[n, m] A; // dummy variables for any extra graph component intercepts  
   int<lower=1> n_edges; 
-  int<lower=1, upper=n> node1[n_edges];
-  int<lower=1, upper=n> node2[n_edges];
+  array[n_edges] int<lower=1, upper=n> node1;
+  array[n_edges] int<lower=1, upper=n> node2;
   vector[n_edges] weight;
-  int<lower=1, upper=k> comp_id[n]; 
+  array[n] int<lower=1, upper=k> comp_id; 
   vector[k] inv_sqrt_scale_factor; // can be a vector of ones, as a placeholder
 
   // ESF 
@@ -103,17 +103,17 @@
   real<lower=0> slab_df;
 
   // CAR
-  real car_rho_lims[2];
+  array[2] real car_rho_lims;
   int<lower=0,upper=1> car;
 
   // SAR
   int nImW_w;
   int nW;
   vector[nImW_w] ImW_w;
-  int ImW_v[nImW_w];
-  int ImW_u[n + 1];
-  int Widx[nW];
+  array[nImW_w] int ImW_v;
+  array[n + 1] int ImW_u;
+  array[nW] int Widx;
   vector[n] eigenvalues_w;
-  real sar_rho_lims[2];
+  array[2] real sar_rho_lims;
   int<lower=0,upper=1> sar;
 

inst/stan/parts/model.stan

@@ -1,5 +1,5 @@
 // parameter models
-  vector[n] x_true_transform[dx_me];
+  array[dx_me] vector[n] x_true_transform;
   x_true_transform = x_true;
   target += normal_lpdf(intercept | prior_alpha[1], prior_alpha[2]);
   if (dx_all) target += normal_lpdf(append_row(gamma, beta) | prior_beta_location, prior_beta_scale);

inst/stan/parts/params.stan

@@ -1,35 +1,35 @@
   //ICAR
   vector[m] alpha_phi;  
   vector[type ? n : 0] phi_tilde;
-  real<lower=0> spatial_scale[type ? 1 : 0];
+  array[type ? 1 : 0] real<lower=0> spatial_scale;
   vector[type > 1 ? n : 0] theta_tilde;
-  real<lower=0> theta_scale[type == 2];
-  real<lower=0,upper=1> rho[type == 3];
+  array[type == 2] real<lower=0> theta_scale;
+  array[type == 3] real<lower=0,upper=1> rho;
   // ESF
-  real<lower=0> aux1_global[dev ? 1 : 0];
-  real<lower=0> aux2_global[dev ? 1 : 0];
+  array[dev ? 1 : 0] real<lower=0> aux1_global;
+  array[dev ? 1 : 0] real<lower=0> aux2_global;
   vector<lower=0>[dev] aux1_local;
   vector<lower=0>[dev] aux2_local;
-  real<lower=0> caux[dev ? 1 : 0];
+  array[dev ? 1 : 0] real<lower=0> caux;
   vector[dev] z;
   // CAR/SAR
   vector[(car||sar) && !is_auto_gaussian ? n : 0] log_lambda;
-  real<lower=0> car_scale[car ? 1 : 0];
-  real<lower=car_rho_lims[1], upper=car_rho_lims[2]> car_rho[car ? 1 : 0];     
-  real<lower=0> sar_scale[sar ? 1 : 0];
-  real<lower=sar_rho_lims[1], upper=sar_rho_lims[2]> sar_rho[sar ? 1 : 0];     
+  array[car ? 1 : 0] real<lower=0> car_scale;
+  array[car ? 1 : 0] real<lower=car_rho_lims[1], upper=car_rho_lims[2]> car_rho;     
+  array[sar ? 1 : 0] real<lower=0> sar_scale;
+  array[sar ? 1 : 0] real<lower=sar_rho_lims[1], upper=sar_rho_lims[2]> sar_rho;     
   // GLM
 // parameters for the process model //
   real intercept;
   vector[dwx] gamma;
   vector[dx_all] beta;
-  real<lower=0> nu[is_student]; 
-  real<lower=0> sigma[has_sigma];
+  array[is_student] real<lower=0> nu; 
+  array[has_sigma] real<lower=0> sigma;
 // for partial pooling across groups/geographies
   vector[n_ids] alpha_re;
-  real<lower=0> alpha_tau[has_re];
+  array[has_re] real<lower=0> alpha_tau;
 // observational error models //
-  vector<lower=bounds[1],upper=bounds[2]>[n] x_true[dx_me];
+  array[dx_me] vector<lower=bounds[1],upper=bounds[2]>[n] x_true;
 //  vector<lower=bounds[1],upper=bounds[2]>[dx_me] mu_x_true;
     vector[dx_me] mu_x_true;
   vector<lower=0>[dx_me] sigma_x_true;

inst/stan/parts/priors.stan

@@ -42,9 +42,9 @@
 		      real sigma,
 		      real rho,
 		      vector ImW,
-		      int[] ImW_v,
-		      int[] ImW_u,
-		      int[] Widx,
+		      array[] int ImW_v,
+		      array[] int ImW_u,
+		      array[] int Widx,
 		      vector lambda,
 		      int n) {
   vector[n] z = y - mu;
@@ -80,7 +80,7 @@
 */
 real car_normal_lpdf(vector y, vector mu,
 		     real tau, real rho,
-		     vector ImC, int[] ImC_v, int[] ImC_u, int[] Cidx,
+		     vector ImC, array[] int ImC_v, array[] int ImC_u, array[] int Cidx,
 		     vector D_inv, real log_det_D_inv, vector lambda,
 		     int n) {
   vector[n] z = y - mu;
@@ -119,7 +119,7 @@ real car_normal_lpdf(vector y, vector mu,
  */
 real wcar_normal_lpdf(vector y, vector mu,
 		      real tau, real rho,
-		      vector A_w, int[] A_v, int[] A_u,
+		      vector A_w, array[] int A_v, array[] int A_u,
 		      vector D_inv, real log_det_D_inv,
 		      vector lambda,
 		      int n) {
@@ -143,8 +143,8 @@ real wcar_normal_lpdf(vector y, vector mu,
  */
 real auto_normal_lpdf(vector y, vector mu,
 		      real tau, real rho,
-		      vector Ax_w, int[] Ax_v, int[] Ax_u,
-		      int[] Cidx,
+		      vector Ax_w, array[] int Ax_v, array[] int Ax_u,
+		      array[] int Cidx,
 		      vector D_inv, real log_det_D_inv,
 		      vector lambda,
 		      int n, int WCAR) {
@@ -171,8 +171,8 @@ real auto_normal_lpdf(vector y, vector mu,
  * @return Log probability density of ICAR prior up to additive constant
  **/
 real icar_normal_lpdf(vector phi, real spatial_scale,
-              int[] node1, int[] node2, 
-              int k, int[] group_size, int[] group_idx,
+              array[] int node1, array[] int node2, 
+              int k, array[] int group_size, array[] int group_idx,
               int has_theta) {
   real lp;
   int pos=1;
@@ -218,7 +218,7 @@ vector make_phi(vector phi_tilde, real phi_scale,
 		real rho,
 		vector inv_sqrt_scale_factor,
 		int n, int k,
-		int[] group_size, int[] group_idx
+		array[] int group_size, array[] int group_idx
               ) {
   vector[n] phi;
   int pos=1;
@@ -244,7 +244,7 @@ vector make_phi(vector phi_tilde, real phi_scale,
  */
 vector convolve_bym(vector phi, vector theta,
               int n, int k,
-              int[] group_size, int[] group_idx
+              array[] int group_size, array[] int group_idx
               ) {
   vector[n] convolution;
   int pos=1;
@@ -279,7 +279,7 @@ vector convolve_bym(vector phi, vector theta,
 vector convolve_bym2(vector phi_tilde, vector theta_tilde,
           real spatial_scale,
 		      int n, int k,
-		      int[] group_size, int[] group_idx,
+		      array[] int group_size, array[] int group_idx,
 		      real rho, vector inv_sqrt_scale_factor
 		      ) {
   vector[n] convolution;

inst/stan/parts/trans_data.stan

@@ -17,7 +17,7 @@
   is_binomial = family == 4;
   is_auto_gaussian = family == 5 || family == 6;  
   has_sigma  =  family < 3;
-  has_offset = sum(offset) != 0;
+  has_offset = sum(input_offset) != 0;
   dx_all = dx_obs + dx_me;
   has_me = dx_me > 0;
 

inst/stan/parts/trans_params_declaration.stan

@@ -1,7 +1,7 @@
   // ESF
   vector[dev] beta_ev;
   vector[dev ? n : 0] esf;  
-  real error_scale[dev ? 1 : 0];  
+  array[dev ? 1 : 0] real error_scale;  
   // ICAR
   vector[type ? n : 0] phi;
   vector[type > 1 ? n : 0] theta;
@@ -13,6 +13,6 @@
   if (dx_obs) x_all[,x_obs_idx] = x_obs;
   if (dx_me) for (j in 1:dx_me) x_all[ ,x_me_idx[j]] = x_true[j];
   if (center_x) for (j in 1:dx_all) x_all[,j] = x_all[,j] - mean(x_all[,j]);
-  if (!car) fitted = offset + intercept;  
+  if (!car) fitted = input_offset + intercept;  
 
 

inst/stan/parts/trans_params_expression_auto-model.stan

@@ -10,8 +10,8 @@
       }
     if (dx_all) log_lambda_mu += x_all * beta;
     if (is_auto_gaussian) {
-      fitted = offset + log_lambda_mu;
+      fitted = input_offset + log_lambda_mu;
     } else {
-      fitted = offset + log_lambda;
+      fitted = input_offset + log_lambda;
     }
   }

inst/stan/parts/trans_params_expression_car.stan

@@ -10,8 +10,8 @@
       }
     if (dx_all) log_lambda_mu += x_all * beta;
     if (is_auto_gaussian) {
-      fitted = offset + log_lambda_mu;
+      fitted = input_offset + log_lambda_mu;
     } else {
-      fitted = offset + log_lambda;
+      fitted = input_offset + log_lambda;
     }
   }