Skip to content

Commit

Permalink
Update deprecated syntax for future rstan compatibility (#453)
Browse files Browse the repository at this point in the history
* Update array syntax

* CDF to LCDF
  • Loading branch information
andrjohns authored Sep 18, 2023
1 parent 5f77fb4 commit fac8567
Show file tree
Hide file tree
Showing 22 changed files with 135 additions and 135 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down
22 changes: 11 additions & 11 deletions inst/stan/data/delays.stan
Original file line number Diff line number Diff line change
@@ -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)
14 changes: 7 additions & 7 deletions inst/stan/data/generation_time.stan
Original file line number Diff line number Diff line change
@@ -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;
2 changes: 1 addition & 1 deletion inst/stan/data/observation_model.stan
Original file line number Diff line number Diff line change
@@ -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
Expand Down
2 changes: 1 addition & 1 deletion inst/stan/data/observations.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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)
2 changes: 1 addition & 1 deletion inst/stan/data/rt.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
18 changes: 9 additions & 9 deletions inst/stan/data/simulation_delays.stan
Original file line number Diff line number Diff line change
@@ -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
8 changes: 4 additions & 4 deletions inst/stan/data/simulation_observation_model.stan
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions inst/stan/data/simulation_rt.stan
Original file line number Diff line number Diff line change
@@ -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
Expand Down
44 changes: 22 additions & 22 deletions inst/stan/dist_fit.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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];
Expand All @@ -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)
);
}
}
Expand Down
30 changes: 15 additions & 15 deletions inst/stan/estimate_infections.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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;
Expand All @@ -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,
Expand Down
14 changes: 7 additions & 7 deletions inst/stan/estimate_secondary.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 {
Expand Down Expand Up @@ -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);
Expand Down
14 changes: 7 additions & 7 deletions inst/stan/estimate_truncation.stan
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
}
Expand Down
22 changes: 11 additions & 11 deletions inst/stan/functions/delays.stan
Original file line number Diff line number Diff line change
@@ -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)) {
Expand All @@ -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
Expand Down Expand Up @@ -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) {
Expand Down
Loading

0 comments on commit fac8567

Please sign in to comment.