#############################################################################
# This script contains the function used to optimize the effects of weather,
# given a weather time series of temperature and relative humidity normalised 
# measurements in order to achieve a certain variability
#############################################################################
compute_beta_var <- function(param, df, objective_var) {
  
  e_Te <- param
  e_RH <- param
  
  # Seasonal transmission rate
  x <- e_Te * df$Te_norm + e_RH * df$RH_norm
  
  # Numerical stability
  Beta_seas <- exp(pmin(50, x))
  
  mean_beta <- mean(Beta_seas, na.rm = TRUE)
  if (!is.finite(mean_beta) || mean_beta <= 0) return(1e12)
  
  # Variability metric
  p_low  <- quantile(Beta_seas, 0.05, na.rm = TRUE)
  p_high <- quantile(Beta_seas, 0.95, na.rm = TRUE)
  
  beta_var <- 100 * ((p_high - p_low) / 2) / mean_beta
  
  # Smooth calibration objective
  (beta_var - objective_var)^2
}

compute_TWO_beta_vars <- function(param, df, objective_var) {

  e_current <- param[1]
  e_lag     <- param[2]

  df <- df[-1,]
  x <- e_current * df$Te_norm +
       e_current * df$RH_norm +
       e_lag     * df$Te_norm_Lag +
       e_lag     * df$RH_norm_Lag

  Beta_seas <- exp(pmin(50, x))

  mean_beta <- mean(Beta_seas, na.rm = TRUE)

  if (!is.finite(mean_beta) || mean_beta <= 0){
    return(1e12)}

  p_low  <- quantile(Beta_seas, 0.05, na.rm = TRUE)
  p_high <- quantile(Beta_seas, 0.95, na.rm = TRUE)

  beta_var <- 100 * ((p_high - p_low) / 2) / mean_beta

  (beta_var - objective_var)^2
}

TMP_compute_beta_var <- function(param, df, objective_var) {
  
  e_Te <- param[1]
  e_RH <- param[2]
  
  # Seasonal transmission rate
  x <- e_Te * df$Te_norm + e_RH * df$RH_norm
  
  # Numerical stability
  Beta_seas <- exp(pmin(50, x))
  
  mean_beta <- mean(Beta_seas, na.rm = TRUE)
  if (!is.finite(mean_beta) || mean_beta <= 0) return(1e12)
  
  # Variability metric
  p_low  <- quantile(Beta_seas, 0.05, na.rm = TRUE)
  p_high <- quantile(Beta_seas, 0.95, na.rm = TRUE)
  
  beta_var <- 100 * ((p_high - p_low) / 2) / mean_beta
  
  # Smooth calibration objective
  (beta_var - objective_var)^2
}