# '
# Creates the Air pollution summary statistics in Supplementary Table SI1.
# Uses panel data stored in "bld/data",
# variables names from "src/code/specifications.json"
# and saves the table as a tex file to "bld/tables".
# 
# '

rm(list=ls())
options(scipen = 999) # Turn off scientific notation
options(warn = 0) # 0 to turn on again, -1 to turn off

# Load libraries
library(tidyverse) # More intuitive data wrangling
library(xtable) # Create flexible tex files
library(rjson) # Read JSON files in R

# Load model specification
model <- fromJSON(file = paste0("src/code/specifications.json"))

# Load analysis data
analysis <- readRDS(paste0("bld/data/analysis_year.RDS"))

# Drop duplicates and keep only obs at 100km bandwidth
dat <- analysis %>%
  filter(icd == "J") %>%
  # distinct(mcode, date, .keep_all = TRUE) %>%
  arrange(mcode, year)
dat <- subset(dat, dat$dist2cutoff <= model$distance_to_border & dat$year > 2002)

# Write table header
final_table <- data.frame((matrix(nrow = 1, ncol = 6)))
final_table[1] <- c("Variable")
final_table[2] <- c("Outside")
final_table[3] <- c("Inside")
final_table[4] <- c("Outside")
final_table[5] <- c("Inside")
final_table[6] <- c("quality guidelines")

# Before period
d <- subset(dat, dat$year < 2007)
# Fill data into table conditional on model.
model_table = data.frame(matrix(nrow = 22, ncol = 6)) 
model_table[1] <- sapply(3:(2 * length(model$mechanisms)), function(i)
  if ((i %% 2) == 0) {
    paste("")
  } else {
    paste(sub("\\\\ ", "", model$mechanisms[[(1+floor(i/2))]], fixed=TRUE))
  }
)
model_table[2] <- sapply(3:(2 * length(model$mechanisms)), function(i)
  if ((i %% 2) == 0) {
    paste0("(", round(sd(d[[names(model$mechanisms[(floor(i/2))])]][d$treated == 0 & d$treated.fires == 0], na.rm = TRUE), 2), ")")
  } else {
    round(mean(d[[names(model$mechanisms[(1+floor(i/2))])]][d$treated == 0 & d$treated.fires == 0], na.rm = TRUE), 2)
  }
)
model_table[3] <- sapply(3:(2 * length(model$mechanisms)), function(i)
  if ((i %% 2) == 0) {
    paste0("(", round(sd(d[[names(model$mechanisms[(floor(i/2))])]][d$treated == 1 & d$treated.fires == 1], na.rm = TRUE), 2), ")")
  } else {
    round(mean(d[[names(model$mechanisms[(1+floor(i/2))])]][d$treated == 1 & d$treated.fires == 1], na.rm = TRUE), 2)
  }
)

# After period
d <- subset(dat, dat$year > 2006)
# Fill data into table conditional on model.
model_table[4] <- sapply(3:(2 * length(model$mechanisms)), function(i)
  if ((i %% 2) == 0) {
    paste0("(", round(sd(d[[names(model$mechanisms[(floor(i/2))])]][d$treated == 0 & d$treated.fires == 0], na.rm = TRUE), 2), ")")
  } else {
    round(mean(d[[names(model$mechanisms[(1+floor(i/2))])]][d$treated == 0 & d$treated.fires == 0], na.rm = TRUE), 2)
  }
)
model_table[5] <- sapply(3:(2 * length(model$mechanisms)), function(i)
  if ((i %% 2) == 0) {
    paste0("(", round(sd(d[[names(model$mechanisms[(floor(i/2))])]][d$treated == 1 & d$treated.fires == 1], na.rm = TRUE), 2), ")")
  } else {
    round(mean(d[[names(model$mechanisms[(1+floor(i/2))])]][d$treated == 1 & d$treated.fires == 1], na.rm = TRUE), 2)
  }
)

# Add information on WHO air quality guidelines column to table
model_table[6] <- 
  c(
    "", "",
    "\\multicolumn{1}{c}{5 (annual average)}",
    "\\multicolumn{1}{c}{15 (24h average)}", "", "", #pm2.5
    "\\multicolumn{1}{c}{4,000 (24h average)}", 
    "\\multicolumn{1}{c}{10,000 (8h average)}", "", "", #CO
    "\\multicolumn{1}{c}{10 (annual average)}", 
    "\\multicolumn{1}{c}{25 (24h average)}", "", "", #NO2
    "\\multicolumn{1}{c}{60 (6 months average)}", 
    "\\multicolumn{1}{c}{100 (8h average)}", "", "", #O3
    "\\multicolumn{1}{c}{40 (24h average)}", "", "", ""#, # SO2
    #"\\multicolumn{1}{c}{20 (24h mean)}",
    # "", "",
    # "", "",
    # "", "",
    # "", ""
)

# Add units description to table
model_table[,1] <- gsub('\\$\\}', ' \\\\,(\\\\mu gm^{-3})$}', model_table[,1])
model_table[,1] <- gsub('guideline', 'guideline (\\\\%)', model_table[,1])
# model_table[,1] <- gsub('expansion', 'expansion $(km^2)$', model_table[,1])
# model_table[,1] <- gsub('quantity', 'quantity (1000 t)', model_table[,1])

# Append both tables
final_table = rbind(final_table, model_table)

# Create latex table
tex_final_table = xtable(
  final_table,
  #caption = "Pre- and Post-ASM Mechanism Variable Averages by Treatment",
  #label = "table:who",
  #align = "llrrrrc"
)

# Add note and top row
comment          <- list()
comment$pos      <- as.list(c(0, nrow(final_table)))
comment$command  <- 
  as.vector(
    c(
      paste0(
        "\\begin{tabular}{lrrrrc} \n",
        "\\\\[-1.8ex]\\hline \n", "\\hline \\\\[-1.8ex] \n",
        "& \\multicolumn{2}{c}{Before} & \\multicolumn{2}{c}{After} & \\multicolumn{1}{c}{WHO \\citep{WHO.2021.Book} air} \\\\ \n"
      ),
      paste0(
        "\\hline \n", "\\hline \\\\[-1.8ex] \n",
        "\\end{tabular} \n",
        "\\begin{tablenotes}[para,flushleft] \n",
        "\\textit{Note:} Numbers represent averages of each group in each period. The outside group refers to municipalities whose capitals
        lie outside the Amazon biome but no further than 100 km away from the border (correspondingly for the inside group).
        The before period is from 2003 to 2006, and the after period from 2007 to 2017.
        Number of fires are the sum of fires per year per municipality. 
        Air pollution concentrations (in micrograms per cubic meter) are annual averages from six hour interval measurements. 
        Above guideline refers to the percentage of time intervals per year where the respective air pollution levels exceeded
        the WHO short-term exposure guidelines.
        PM\\textsubscript{2.5} is particulate matter with a diameter of less than 2.5 $\\mu m$; $CO$ is carbon monoxide;
        NO\\textsubscript{2} is nitrogen dioxide; O\\textsubscript{3} is ozone; and SO\\textsubscript{2} is sulfur dioxide.
        Standard deviations in parentheses. \\\\ \n",
        "\\end{tablenotes} \n"
      )
    ), mode="character"
)

# Save the latex table
print(
  tex_final_table,
  sanitize.text.function = function(x){x},
  file = paste0("bld/tables/who_guideline.tex"),
  include.rownames = FALSE,
  include.colnames = FALSE,
  caption.placement = "top",
  add.to.row = comment,
  hline.after = c(1),
  floating = FALSE,
  only.contents = TRUE
)