Merge 1353bd5 into b003127

DESCRIPTION

@@ -109,7 +109,7 @@ Imports:
     R.utils (>= 2.0.0),
     Rcpp (>= 0.12.0), 
     rlang (>= 0.4.7),
-    rstan (>= 2.21.1),
+    rstan (>= 2.26.0),
     rstantools (>= 2.2.0),
     runner,
     scales,
@@ -135,8 +135,8 @@ LinkingTo:
     Rcpp (>= 0.12.0),
     RcppEigen (>= 0.3.3.3.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.21.1),
-    StanHeaders (>= 2.21.0-5)
+    rstan (>= 2.26.0),
+    StanHeaders (>= 2.26.0)
 Biarch: true
 Config/testthat/edition: 3
 Encoding: UTF-8

inst/stan/data/delays.stan

@@ -1,18 +1,18 @@
   int<lower = 0> delay_n;                     // number of delay distributions
   int<lower = 0> delay_n_p;                   // number of parametric delay distributions
   int<lower = 0> delay_n_np;                  // number of nonparametric delay distributions
-  real delay_mean_mean[delay_n_p];            // prior mean of mean delay distribution
-  real<lower = 0> delay_mean_sd[delay_n_p];   // prior sd of mean delay distribution
-  real<lower = 0> delay_sd_mean[delay_n_p];   // prior sd of sd of delay distribution
-  real<lower = 0> delay_sd_sd[delay_n_p];     // prior sd of sd of delay distribution
-  int<lower = 1> delay_max[delay_n_p];        // maximum delay distribution
-  int<lower = 0> delay_dist[delay_n_p];       // 0 = lognormal; 1 = gamma
+  array[delay_n_p] real delay_mean_mean;            // prior mean of mean delay distribution
+  array[delay_n_p] real<lower = 0> delay_mean_sd;   // prior sd of mean delay distribution
+  array[delay_n_p] real<lower = 0> delay_sd_mean;   // prior sd of sd of delay distribution
+  array[delay_n_p] real<lower = 0> delay_sd_sd;     // prior sd of sd of delay distribution
+  array[delay_n_p] int<lower = 1> delay_max;        // maximum delay distribution
+  array[delay_n_p] int<lower = 0> delay_dist;       // 0 = lognormal; 1 = gamma
   int<lower = 0> delay_np_pmf_max;            // number of nonparametric pmf elements
   vector<lower = 0, upper = 1>[delay_np_pmf_max] delay_np_pmf; // ragged array of fixed PMFs
-  int<lower = 1> delay_np_pmf_groups[delay_n_np + 1];              // links to ragged array
-  int<lower = 0> delay_weight[delay_n_p];
+  array[delay_n_np + 1] int<lower = 1> delay_np_pmf_groups;              // links to ragged array
+  array[delay_n_p] int<lower = 0> delay_weight;
 
   int<lower = 0> delay_types;                     // number of delay types
-  int<lower = 0> delay_types_p[delay_n];          // whether delay types are parametric
-  int<lower = 0> delay_types_id[delay_n];          // whether delay types are parametric
-  int<lower = 0> delay_types_groups[delay_types + 1]; // index of each delay (parametric or non)
+  array[delay_n] int<lower = 0> delay_types_p;          // whether delay types are parametric
+  array[delay_n] int<lower = 0> delay_types_id;          // whether delay types are parametric
+  array[delay_types + 1] int<lower = 0> delay_types_groups; // index of each delay (parametric or non)

inst/stan/data/generation_time.stan

@@ -1,8 +1,8 @@
-  real gt_mean_sd[1];                   // prior sd of mean generation time
-  real gt_mean_mean[1];                 // prior mean of mean generation time
-  real gt_sd_mean[1];                   // prior mean of sd of generation time
-  real gt_sd_sd[1];                     // prior sd of sd of generation time
-  int<lower = 1> gt_max[1];                        // maximum generation time
-  int gt_fixed[1];                      // 0 = variable gt; 1 = fixed gt
-  int gt_dist[1];                    // distribution (0 = lognormal, 1 = gamma)
+  array[1] real gt_mean_sd;                   // prior sd of mean generation time
+  array[1] real gt_mean_mean;                 // prior mean of mean generation time
+  array[1] real gt_sd_mean;                   // prior mean of sd of generation time
+  array[1] real gt_sd_sd;                     // prior sd of sd of generation time
+  array[1] int<lower = 1> gt_max;                        // maximum generation time
+  array[1] int gt_fixed;                      // 0 = variable gt; 1 = fixed gt
+  array[1] int gt_dist;                    // distribution (0 = lognormal, 1 = gamma)
   int gt_weight;                   

inst/stan/data/observation_model.stan

@@ -1,4 +1,4 @@
-  int day_of_week[t - seeding_time]; // day of the week indicator (1 - 7)
+  array[t - seeding_time] int day_of_week; // day of the week indicator (1 - 7)
   int model_type;                    // type of model: 0 = poisson otherwise negative binomial
   real phi_mean;                     // Mean and sd of the normal prior for the
   real phi_sd;                       // reporting process

inst/stan/data/observations.stan

@@ -2,5 +2,5 @@
   int seeding_time;                                 // time period used for seeding and not observed
   int horizon;                                      // forecast horizon
   int future_time;                                  // time in future for Rt
-  int<lower = 0> cases[t - horizon - seeding_time]; // observed cases
+  array[t - horizon - seeding_time] int<lower = 0> cases; // observed cases
   vector<lower = 0>[t] shifted_cases;               // prior infections (for backcalculation)

inst/stan/data/rt.stan

@@ -4,7 +4,7 @@
   real <lower = 0> r_mean;           // prior mean of reproduction number
   real <lower = 0> r_sd;             // prior standard deviation of reproduction number
   int bp_n;                          // no of breakpoints (0 = no breakpoints)
-  int breakpoints[t - seeding_time]; // when do breakpoints occur
+  array[t - seeding_time] int breakpoints; // when do breakpoints occur
   int future_fixed;                  // is underlying future Rt assumed to be fixed
   int fixed_from;                    // Reference date for when Rt estimation should be fixed
   int pop;                           // Initial susceptible population

inst/stan/data/simulation_delays.stan

@@ -1,18 +1,18 @@
   int<lower = 0> delay_n;                  // number of delay distribution distributions
   int<lower = 0> delay_n_p;                // number of parametric delay distributions
   int<lower = 0> delay_n_np;                // number of nonparametric delay distributions
-  real delay_mean[n, delay_n_p]; // prior mean of mean delay distribution
-  real<lower = 0> delay_sd[n, delay_n_p];   // prior sd of sd of delay distribution
-  int<lower = 1> delay_max[delay_n_p];          // maximum delay distribution
-  int<lower = 0> delay_dist[delay_n_p];       // 0 = lognormal; 1 = gamma
+  array[n, delay_n_p] real delay_mean; // prior mean of mean delay distribution
+  array[n, delay_n_p] real<lower = 0> delay_sd;   // prior sd of sd of delay distribution
+  array[delay_n_p] int<lower = 1> delay_max;          // maximum delay distribution
+  array[delay_n_p] int<lower = 0> delay_dist;       // 0 = lognormal; 1 = gamma
   int<lower = 0> delay_np_pmf_max;          // number of nonparametric pmf elements
   vector<lower = 0, upper = 1>[delay_np_pmf_max] delay_np_pmf; // ragged array of fixed PMFs
-  int<lower = 1> delay_np_pmf_groups[delay_n_np + 1];      // links to ragged array
-  int delay_weight[delay_n_p];
+  array[delay_n_np + 1] int<lower = 1> delay_np_pmf_groups;      // links to ragged array
+  array[delay_n_p] int delay_weight;
 
   int<lower = 0> delay_types;                     // number of delay types
-  int<lower = 0> delay_types_p[delay_n];          // whether delay types are parametric
-  int<lower = 0> delay_types_id[delay_n];          // whether delay types are parametric
-  int<lower = 0> delay_types_groups[delay_types + 1]; // index of each delay (parametric or non)
+  array[delay_n] int<lower = 0> delay_types_p;          // whether delay types are parametric
+  array[delay_n] int<lower = 0> delay_types_id;          // whether delay types are parametric
+  array[delay_types + 1] int<lower = 0> delay_types_groups; // index of each delay (parametric or non)
 
   int<lower = 0> delay_id; // id of generation time

inst/stan/data/simulation_observation_model.stan

@@ -1,8 +1,8 @@
-  int day_of_week[t - seeding_time]; // day of the week indicator (1 - 7)
+  array[t - seeding_time] int day_of_week; // day of the week indicator (1 - 7)
   int week_effect;                   // should a day of the week effect be estimated
-  real<lower = 0> day_of_week_simplex[n, week_effect];
+  array[n, week_effect] real<lower = 0> day_of_week_simplex;
   int obs_scale;
-  real<lower = 0, upper = 1> frac_obs[n, obs_scale];
+  array[n, obs_scale] real<lower = 0, upper = 1> frac_obs;
   int model_type;
-  real<lower = 0> rep_phi[n, model_type];  // overdispersion of the reporting process
+  array[n, model_type] real<lower = 0> rep_phi;  // overdispersion of the reporting process
   int<lower = 0> trunc_id; // id of truncation

inst/stan/data/simulation_rt.stan

@@ -1,5 +1,5 @@
-  real initial_infections[seeding_time ? n : 0, 1]; // initial logged infections
-  real initial_growth[seeding_time > 1 ? n : 0, 1]; //initial growth
+  array[seeding_time ? n : 0, 1] real initial_infections; // initial logged infections
+  array[seeding_time > 1 ? n : 0, 1] real initial_growth; //initial growth
 
   matrix[n, t - seeding_time] R; // reproduction number
   int pop;                       // susceptible population

inst/stan/dist_fit.stan

@@ -3,15 +3,15 @@ data {
   int N;
   vector[N] low;
   vector[N] up;
-  real lam_mean[dist == 0];
-  real prior_mean[dist > 0];
-  real prior_sd[dist > 0];
-  real<lower = 0> par_sigma[dist == 2];
+  array[dist == 0] real lam_mean;
+  array[dist > 0] real prior_mean;
+  array[dist > 0] real prior_sd;
+  array[dist == 2] real<lower = 0> par_sigma;
 }
 
 transformed data {
-  real prior_alpha[dist == 2];
-  real prior_beta[dist == 2];
+  array[dist == 2] real prior_alpha;
+  array[dist == 2] real prior_beta;
 
   if (dist == 2) {
     prior_alpha[1] = (prior_mean[1] / prior_sd[1])^2;
@@ -20,16 +20,16 @@ transformed data {
 }
 
 parameters {
-  real<lower = 0> lambda[dist == 0];
-  real mu[dist == 1];
-  real<lower = 0> sigma[dist == 1];
-  real<lower = 0> alpha_raw[dist == 2];
-  real<lower = 0> beta_raw[dist == 2];
+  array[dist == 0] real<lower = 0> lambda;
+  array[dist == 1] real mu;
+  array[dist == 1] real<lower = 0> sigma;
+  array[dist == 2] real<lower = 0> alpha_raw;
+  array[dist == 2] real<lower = 0> beta_raw;
 }
 
 transformed parameters{
-  real<lower = 0> alpha[dist == 2];
-  real<lower = 0> beta[dist == 2];
+  array[dist == 2] real<lower = 0> alpha;
+  array[dist == 2] real<lower = 0> beta;
 
   if (dist == 2) {
     alpha[1] = prior_alpha[1] + par_sigma[1] * alpha_raw[1];
@@ -50,19 +50,19 @@ model {
 
   for(i in 1:N){
     if (dist == 0) {
-      target += log(
-        exponential_cdf(up[i] , lambda) -
-        exponential_cdf(low[i], lambda)
+      target += log_diff_exp(
+        exponential_lcdf(up[i] | lambda),
+        exponential_lcdf(low[i] | lambda)
       );
     } else if (dist == 1) {
-      target += log(
-        lognormal_cdf(up[i], mu, sigma) -
-        lognormal_cdf(low[i], mu, sigma)
+      target += log_diff_exp(
+        lognormal_lcdf(up[i] | mu, sigma),
+        lognormal_lcdf(low[i] | mu, sigma)
       );
     } else if (dist == 2) {
-      target += log(
-        gamma_cdf(up[i], alpha, beta) -
-        gamma_cdf(low[i], alpha, beta)
+      target += log_diff_exp(
+        gamma_lcdf(up[i] | alpha, beta),
+        gamma_lcdf(low[i] | alpha, beta)
       );
     }
   }

inst/stan/estimate_infections.stan

@@ -30,29 +30,29 @@ transformed data{
   real r_logmean = log(r_mean^2 / sqrt(r_sd^2 + r_mean^2));
   real r_logsd = sqrt(log(1 + (r_sd^2 / r_mean^2)));
 
-  int delay_type_max[delay_types] = get_delay_type_max(
+  array[delay_types] int delay_type_max = get_delay_type_max(
     delay_types, delay_types_p, delay_types_id,
     delay_types_groups, delay_max, delay_np_pmf_groups
   );
 }
 
 parameters{
   // gaussian process
-  real<lower = ls_min,upper=ls_max> rho[fixed ? 0 : 1];  // length scale of noise GP
-  real<lower = 0> alpha[fixed ? 0 : 1];    // scale of of noise GP
+  array[fixed ? 0 : 1] real<lower = ls_min,upper=ls_max> rho;  // length scale of noise GP
+  array[fixed ? 0 : 1] real<lower = 0> alpha;    // scale of of noise GP
   vector[fixed ? 0 : M] eta;               // unconstrained noise
   // Rt
   vector[estimate_r] log_R;                // baseline reproduction number estimate (log)
-  real initial_infections[estimate_r] ;    // seed infections
-  real initial_growth[estimate_r && seeding_time > 1 ? 1 : 0]; // seed growth rate
-  real<lower = 0> bp_sd[bp_n > 0 ? 1 : 0]; // standard deviation of breakpoint effect
-  real bp_effects[bp_n];                   // Rt breakpoint effects
+  array[estimate_r] real initial_infections ;    // seed infections
+  array[estimate_r && seeding_time > 1 ? 1 : 0] real initial_growth; // seed growth rate
+  array[bp_n > 0 ? 1 : 0] real<lower = 0> bp_sd; // standard deviation of breakpoint effect
+  array[bp_n] real bp_effects;                   // Rt breakpoint effects
   // observation model
-  real delay_mean[delay_n_p];         // mean of delays
-  real<lower = 0> delay_sd[delay_n_p];  // sd of delays
+  array[delay_n_p] real delay_mean;         // mean of delays
+  array[delay_n_p] real<lower = 0> delay_sd;  // sd of delays
   simplex[week_effect] day_of_week_simplex;// day of week reporting effect
-  real<lower = 0, upper = 1> frac_obs[obs_scale];     // fraction of cases that are ultimately observed
-  real<lower = 0> rep_phi[model_type];     // overdispersion of the reporting process
+  array[obs_scale] real<lower = 0, upper = 1> frac_obs;     // fraction of cases that are ultimately observed
+  array[model_type] real<lower = 0> rep_phi;     // overdispersion of the reporting process
 }
 
 transformed parameters {
@@ -154,9 +154,9 @@ model {
 }
 
 generated quantities {
-  int imputed_reports[ot_h];
+  array[ot_h] int imputed_reports;
   vector[estimate_r > 0 ? 0: ot_h] gen_R;
-  real r[ot_h];
+  array[ot_h] real r;
   real gt_mean;
   real gt_var;
   vector[return_likelihood ? ot : 0] log_lik;
@@ -167,9 +167,9 @@ generated quantities {
     r = R_to_growth(R, gt_mean, gt_var);
   } else {
     // sample generation time
-    real delay_mean_sample[delay_n_p] =
+    array[delay_n_p] real delay_mean_sample =
       normal_rng(delay_mean_mean, delay_mean_sd);
-    real delay_sd_sample[delay_n_p] =
+    array[delay_n_p] real delay_sd_sample =
       normal_rng(delay_sd_mean, delay_sd_sd);
     vector[delay_type_max[gt_id]] sampled_gt_rev_pmf = get_delay_rev_pmf(
       gt_id, delay_type_max[gt_id], delay_types_p, delay_types_id,

inst/stan/estimate_secondary.stan

@@ -8,7 +8,7 @@ functions {
 
 data {
   int t;                             // time of observations
-  int<lower = 0> obs[t];             // observed secondary data
+  array[t] int<lower = 0> obs;             // observed secondary data
   vector[t] primary;                 // observed primary data
   int burn_in;                       // time period to not use for fitting
 #include data/secondary.stan
@@ -17,19 +17,19 @@ data {
 }
 
 transformed data{
-  int delay_type_max[delay_types] = get_delay_type_max(
+  array[delay_types] int delay_type_max = get_delay_type_max(
     delay_types, delay_types_p, delay_types_id,
     delay_types_groups, delay_max, delay_np_pmf_groups
   );
 }
 
 parameters{
   // observation model
-  real delay_mean[delay_n_p];
-  real<lower = 0> delay_sd[delay_n_p];      // sd of delays
+  array[delay_n_p] real delay_mean;
+  array[delay_n_p] real<lower = 0> delay_sd;      // sd of delays
   simplex[week_effect] day_of_week_simplex;  // day of week reporting effect
-  real<lower = 0, upper = 1> frac_obs[obs_scale];   // fraction of cases that are ultimately observed
-  real<lower = 0> rep_phi[model_type];   // overdispersion of the reporting process
+  array[obs_scale] real<lower = 0, upper = 1> frac_obs;   // fraction of cases that are ultimately observed
+  array[model_type] real<lower = 0> rep_phi;   // overdispersion of the reporting process
 }
 
 transformed parameters {
@@ -89,7 +89,7 @@ model {
 }
 
 generated quantities {
-  int sim_secondary[t - burn_in];
+  array[t - burn_in] int sim_secondary;
   vector[return_likelihood > 1 ? t - burn_in : 0] log_lik;
   // simulate secondary reports
   sim_secondary = report_rng(secondary[(burn_in + 1):t], rep_phi, model_type);

inst/stan/estimate_truncation.stan

@@ -7,16 +7,16 @@ functions {
 data {
   int t;
   int obs_sets;
-  int obs[t, obs_sets];
-  int obs_dist[obs_sets];
+  array[t, obs_sets] int obs;
+  array[obs_sets] int obs_dist;
 #include data/delays.stan
 }
 transformed data{
   int trunc_id = 1;
-  int<lower = 1> end_t[obs_sets];
-  int<lower = 1> start_t[obs_sets];
+  array[obs_sets] int<lower = 1> end_t;
+  array[obs_sets] int<lower = 1> start_t;
 
-  int delay_type_max[delay_types];
+  array[delay_types] int delay_type_max;
   delay_type_max = get_delay_type_max(
     delay_types, delay_types_p, delay_types_id,
     delay_types_groups, delay_max, delay_np_pmf_groups
@@ -28,8 +28,8 @@ transformed data{
   }
 }
 parameters {
-  real delay_mean[delay_n_p];
-  real<lower = 0> delay_sd[delay_n_p];      // sd of delays
+  array[delay_n_p] real delay_mean;
+  array[delay_n_p] real<lower = 0> delay_sd;      // sd of delays
   real<lower=0> phi;
   real<lower=0> sigma;
 }

inst/stan/functions/delays.stan

@@ -1,8 +1,8 @@
-int[] get_delay_type_max(
-  int delay_types, int[] delay_types_p, int[] delay_types_id,
-  int[] delay_types_groups, int[] delay_max, int[] delay_np_pmf_groups
+array[] int get_delay_type_max(
+  int delay_types, array[] int delay_types_p, array[] int delay_types_id,
+  array[] int delay_types_groups, array[] int delay_max, array[] int delay_np_pmf_groups
 ) {
-  int ret[delay_types];
+  array[delay_types] int ret;
   for (i in 1:delay_types) {
     ret[i] = 1;
     for (j in delay_types_groups[i]:(delay_types_groups[i + 1] - 1)) {
@@ -18,10 +18,10 @@ int[] get_delay_type_max(
 }
 
 vector get_delay_rev_pmf(
-  int delay_id, int len, int[] delay_types_p, int[] delay_types_id,
-  int[] delay_types_groups, int[] delay_max,
-  vector delay_np_pmf, int[] delay_np_pmf_groups,
-  real[] delay_mean, real[] delay_sigma, int[] delay_dist,
+  int delay_id, int len, array[] int delay_types_p, array[] int delay_types_id,
+  array[] int delay_types_groups, array[] int delay_max,
+  vector delay_np_pmf, array[] int delay_np_pmf_groups,
+  array[] real delay_mean, array[] real delay_sigma, array[] int delay_dist,
   int left_truncate, int reverse_pmf, int cumulative
 ) {
   // loop over delays
@@ -75,9 +75,9 @@ vector get_delay_rev_pmf(
 }
 
 
-void delays_lp(real[] delay_mean, real[] delay_mean_mean, real[] delay_mean_sd,
-               real[] delay_sd, real[] delay_sd_mean, real[] delay_sd_sd,
-               int[] delay_dist, int[] weight) {
+void delays_lp(array[] real delay_mean, array[] real delay_mean_mean, array[] real delay_mean_sd,
+               array[] real delay_sd, array[] real delay_sd_mean, array[] real delay_sd_sd,
+               array[] int delay_dist, array[] int weight) {
     int mean_delays = num_elements(delay_mean);
     int sd_delays = num_elements(delay_sd);
     if (mean_delays) {