#clear environment
rm(list = ls())


library(knitr) #used to print the results
library(tidyverse) #used for everything


# import data -------------------------------------------------------------

load("Rdata/isotopes.Rdata")

iso.df <- iso.df %>%
  select(Name,Species,Facility,Type,Matrix,isotope.raw,isotope.expected) %>%
  arrange(Species,Facility,Name)

write.csv(iso.df,"export/raw_data.csv",row.names=FALSE)

qc.df <- data.frame(matrix(nrow=2,ncol=0))
row <- 0

for (myspecies in c("N","C")){
  QC <- iso.df %>%
    filter(Type=="Quality Control") %>%
    filter(Species==myspecies) %>%
    pull(isotope.expected)
  row <- row+1
  qc.df[row,"Species"] <- myspecies
  qc.df[row,"min"] <- min(QC)
  qc.df[row,"max"] <- max(QC)
  qc.df[row,"mean"] <- mean(QC)
}

# normalization -----------------------------------------------------------

isorange.df <- data.frame()

norm.list <- iso.df %>%
  filter(Type=="Normalization") %>%
  pull(Name) %>%
  unique() %>%
  as.character()

onepoint.list <- combn(norm.list,1)
twopoint.list <- combn(norm.list,2)
threepoint.list <- combn(norm.list,3)
fourpoint.list <- combn(norm.list,4)

point.list <- list(onepoint.list,twopoint.list,threepoint.list,fourpoint.list)

name.list <- c("one.point","two.point","three.point","four.point")

facility.list <- iso.df %>%
  pull(Facility) %>%
  unique()

results.list <- list()

name_counter <- 0

for (list in point.list){
  name_counter <- name_counter+1
  print(name_counter)
  for (facility in facility.list){
    for (myspecies in c("N","C")){
      for (col in 1:ncol(list)){
        
        normalization <- list[,col]
        
        temp.df <- iso.df %>%
          filter(Name %in% normalization) %>%
          filter(Species==myspecies)%>%
          filter(Facility==facility)
        
        y <- temp.df$isotope.expected
        x <- temp.df$isotope.raw
        
        lm <- lm(y~x)
        slope <- coef(lm)[2]
        intercept <- coef(lm)[1]
        
        if (slope<0.5){
          slope <- 1
          intercept <- mean(y)-mean(x)
        }
        
        qc_temp <- qc.df %>%
          filter(Species==myspecies)
        
        if (min(y)>qc_temp["min"] | max(y)<qc_temp["max"]){
          extrapolation <- "Extrapolated"
        } else {
          extrapolation <- "Bounded"
        }
        
        method <- name.list[name_counter]
        
        results.df <- iso.df %>%
          select(-Matrix)
        
        matrix <- unique(temp.df$Matrix)
        
        if (length(matrix)>1){
          matrix <- "Mixed"
        }
        
        for (myrow in 1:nrow(results.df)){
          if (results.df[myrow,"Species"]==myspecies & results.df[myrow,"Facility"]==facility){
            results.df[myrow,"Method"] <- method
            results.df[myrow,"Normalization"] <- col
            results.df[myrow,"Extrapolation"] <- extrapolation
            results.df[myrow,"Matrix"] <- matrix
            results.df[myrow, "iso.range"] <- max(y)-min(y)
            results.df[myrow,"isotope.norm"] <- results.df[myrow, "isotope.raw"] * slope + intercept
            results.df[myrow,'isotope.deviation'] <- results.df[myrow,"isotope.norm"] - results.df[myrow,"isotope.expected"]
            results.df[myrow,'qc.difference'] <- results.df[myrow,"isotope.expected"] - qc_temp[4]
            results.df[myrow,"all.names"] <- paste(normalization,collapse = ' ')
          }
        }
        
        results.list <- append(list(results.df),results.list)
        
      }
    }
  }
}

normalizations.df <- bind_rows(results.list) %>%
  drop_na()

save(normalizations.df,file="Rdata/normalizations_02102023.Rdata")

# summarize the results ---------------------------------------------------

summary.df <- normalizations.df %>%
  filter(Type=="Quality Control") %>%
  select(Species,Facility,Method,isotope.deviation) %>%
  group_by(Species,Facility,Method) %>%
  summarize_all(mean) %>%
  mutate_if(is.numeric,round,digits=3)

epa_summary.df <- summary.df %>% 
  filter(Facility==facility.list[1]) %>%
  select(-Facility)

csi_summary.df <- summary.df %>% 
  filter(Facility==facility.list[2]) %>%
  select(-Facility)

summary.df <- full_join(epa_summary.df,csi_summary.df)

kable(summary.df)

summary.df <- normalizations.df %>%
  filter(Type=="Quality Control") %>%
  select(Species,Facility,Method,Matrix,Extrapolation,isotope.deviation,all.names) %>%
  group_by(Species,Method,Facility,Matrix,Extrapolation,all.names) %>%
  summarize_all(mean) %>%
  mutate_if(is.numeric,round,digits=3) %>%
  mutate(Method=factor(Method,levels=c("one.point","two.point","three.point","four.point"))) %>%
  arrange(Species,Facility,Method,isotope.deviation)

colnames(summary.df) <- c("Species","Method","Facility","Matrix","Extrapolation","Standards","QC Deviation (permil)")

kable(summary.df)

write.csv(summary.df,"export/all_normalizations.csv",row.names=FALSE)
