###### Paper: ARPALData: an R package for retrieving and analyzing air quality and weather data from ARPA Lombardia (Italy).
###### Authors: After the revision process
###### Journal: Environmental and Ecological Statistics
###### Date: 25/09/2023

##### The following code reproduces the Monte Carlo experiment to compute the computational burden
#####   using serial and parallel download strategies and discussed at Section 4.2 of the paper


##### Setup 
library(ARPALData)
library(tidyverse)

setwd("~/Ricerca/ARPALData_package")
'%notin%' <- Negate('%in%')

##### Define the list of monitoring stations for the experiment
# Currently active
# Activated before 2014
reg <- get_ARPA_Lombardia_AQ_registry()
reg_red <- reg %>%
  filter(is.na(DateStop),DateStart <= "2014-01-01") 

##### Define the sequence of dates for the experiment
Date_seq <- character(length = 40)
for (t in 1:length(Date_seq)) {
  Date_seq[t] <- as.character(as_datetime("2014-01-01") + months(t))
}

##### Simulations parameters
iterations <- 50
NumStats <- c(10,20,50,0)
NumMonths <- c(1,3,6,12,36)


##### Stations for loop
for (s in 1:length(NumStats)) {
  ##### Months for loop
  for (m in 1:length(NumMonths)) {
    num_stats_down <- NumStats[s]
    months_down <- NumMonths[m]
    res <- matrix(data = NA, nrow = iterations, ncol = 12)
    colnames(res) <- c("NumStats", "Months", "Iter", "Seed", "Date_begin", "Date_end",
                       "Begin_np", "End_np", "Time_np", "Begin_p", "End_p", "Time_p")
    ##### Iterations for loop
    for (iter in 1:iterations) {
      cat(paste0("NumStats: ",num_stats_down," & NumMonths: ",months_down, " - Iteration ",iter))
      seed <- iter + 1000
      set.seed(seed)
      ##### Random sampling station IDs
      if (num_stats_down != 0) {
        stats <- sample(x = reg_red$IDStation, size = num_stats_down, replace = F)
      } else {
        stats <- NULL
      }
      ##### Random sampling date sequence
      Date_begin <- sample(x = Date_seq, size = 1, replace = F)
      Date_end <- as.character(as_datetime(Date_begin) + months(months_down))
      stats
      Date_begin
      Date_end
      ##### Non parallel download
      b_download_np <- Sys.time()
      d <- get_ARPA_Lombardia_AQ_data(ID_station = stats,
                                      Date_begin = Date_begin,
                                      Date_end = Date_end,
                                      parallel = FALSE)
      e_download_np <- Sys.time()
      t_download_np <- e_download_np - b_download_np
      ##### Parallel download
      b_download_p <- Sys.time()
      d <- get_ARPA_Lombardia_AQ_data(ID_station = stats,
                                      Date_begin = Date_begin,
                                      Date_end = Date_end,
                                      parallel = TRUE)
      e_download_p <- Sys.time()
      t_download_p <- e_download_p - b_download_p
      
      res[iter,] <- c(num_stats_down, months_down, iter, seed, Date_begin, Date_end,
                      as.character(b_download_np), as.character(e_download_np), t_download_np,
                      as.character(b_download_p), as.character(e_download_p), t_download_p)
    }
    #### Transform to df
    res <- data.frame(res)
    res <- res %>%
      mutate(across(c("NumStats", "Months", "Iter", "Seed", "Time_np", "Time_p"), as.numeric),
             across(c("Date_begin", "Date_end","Begin_np", "End_np", "Begin_p", "End_p"), lubridate::as_datetime))
    save(res, file = paste0("SimsRes_NumStats",num_stats_down,"_NumMonths",months_down,".Rdata"))
    gc()
  }
}

##### Output analysis
NumMonths <- c(1,3,6,12,36)
NumStats <- c(0,10,20,50)
combinations <- expand.grid(NumMonths,NumStats)
ResSims <- vector(mode = "list", length = dim(combinations)[1])
for (i in 1:dim(combinations)[1]) {
  file_name <- paste0("~/R/Sims_timing/SimsRes_NumStats",combinations[i,2],"_NumMonths",combinations[i,1],".RData")
  if (file.exists(file_name)) {
    load(file_name)
    ResSims[[i]] <- res
  } else {
    ResSims[[i]] <- NULL
  }
}
ResSims <- bind_rows(ResSims)

##### Non-parallel: average computing time (seconds)
ResSims %>%
  mutate(Time_np = End_np - Begin_np,
         Time_p = End_p - Begin_p,
         NumStats = factor(NumStats,levels = c("10","20","50","0"), ordered = T),
         Months = factor(Months,levels = c("1","3","6","12","36"), ordered = T)) %>%
  group_by(NumStats,Months) %>%
  summarise(m_Time_np = mean(Time_np)) %>%
  pivot_wider(names_from = Months, values_from = m_Time_np)

##### Non-parallel: average computing time (seconds)
ResSims %>%
  mutate(Time_np = End_np - Begin_np,
         Time_p = End_p - Begin_p,
         NumStats = factor(NumStats,levels = c("10","20","50","0"), ordered = T),
         Months = factor(Months,levels = c("1","3","6","12","36"), ordered = T)) %>%
  group_by(NumStats,Months) %>%
  summarise(s_Time_np = sd(Time_np)) %>%
  pivot_wider(names_from = Months, values_from = s_Time_np)

##### Parallel: average computing time (seconds)
ResSims %>%
  mutate(Time_np = End_np - Begin_np,
         Time_p = End_p - Begin_p,
         NumStats = factor(NumStats,levels = c("10","20","50","0"), ordered = T),
         Months = factor(Months,levels = c("1","3","6","12","36"), ordered = T)) %>%
  group_by(NumStats,Months) %>%
  summarise(m_Time_p = mean(Time_p)) %>%
  pivot_wider(names_from = Months, values_from = m_Time_p)

##### Parallel: average computing time (seconds)
ResSims %>%
  mutate(Time_np = End_np - Begin_np,
         Time_p = End_p - Begin_p,
         NumStats = factor(NumStats,levels = c("10","20","50","0"), ordered = T),
         Months = factor(Months,levels = c("1","3","6","12","36"), ordered = T)) %>%
  group_by(NumStats,Months) %>%
  summarise(s_Time_p = sd(Time_p)) %>%
  pivot_wider(names_from = Months, values_from = s_Time_p)