From fc9e031dd508d41a941ea88d1273f46b39faaa29 Mon Sep 17 00:00:00 2001
From: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Date: Mon, 11 Mar 2024 11:33:04 +0000
Subject: [PATCH] use inv_square instead of 1/pow(..., 2) (#603)

---
 inst/stan/functions/generated_quantities.stan | 2 +-
 inst/stan/functions/observation_model.stan    | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/inst/stan/functions/generated_quantities.stan b/inst/stan/functions/generated_quantities.stan
index 003fbee52..39411a999 100644
--- a/inst/stan/functions/generated_quantities.stan
+++ b/inst/stan/functions/generated_quantities.stan
@@ -33,7 +33,7 @@ array[] real R_to_growth(vector R, real gt_mean, real gt_var) {
   int t = num_elements(R);
   array[t] real r;
   if (gt_var > 0) {
-    real k = gt_var / pow(gt_mean, 2);
+    real k = gt_var * inv_square(gt_mean);
     for (s in 1:t) {
       r[s] = (pow(R[s], k) - 1) / (k * gt_mean);
     }
diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan
index 97ca63967..ac54496c7 100644
--- a/inst/stan/functions/observation_model.stan
+++ b/inst/stan/functions/observation_model.stan
@@ -77,7 +77,7 @@ void report_lp(array[] int cases, array[] int cases_time, vector reports,
     obs_cases = cases;
   }
   if (model_type) {
-    real dispersion = 1 / pow(phi_sd > 0 ? rep_phi[model_type] : phi_mean, 2);
+    real dispersion = inv_square(phi_sd > 0 ? rep_phi[model_type] : phi_mean);
     if (phi_sd > 0) {
       rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,];
     }
@@ -108,7 +108,7 @@ vector report_log_lik(array[] int cases, vector reports,
       log_lik[i] = poisson_lpmf(cases[i] | reports[i]) * weight;
     }
   } else {
-    real dispersion = 1 / pow(rep_phi[model_type], 2);
+    real dispersion = inv_square(rep_phi[model_type]);
     for (i in 1:t) {
       log_lik[i] = neg_binomial_2_lpmf(cases[i] | reports[i], dispersion) * weight;
     }
@@ -121,7 +121,7 @@ array[] int report_rng(vector reports, array[] real rep_phi, int model_type) {
   array[t] int sampled_reports;
   real dispersion = 1e5;
   if (model_type) {
-    dispersion = 1 / pow(rep_phi[model_type], 2);
+    dispersion = inv_square(rep_phi[model_type]);
   }
 
   for (s in 1:t) {