global.R

library(tidyr)
library(skimr)
library(chRoma)
library(tibble)
library(dplyr)
library(data.table)
library(data.tree)
library(plotly)
library(classInt)
library(aws.s3)
library(digest)
library(shiny)
library(DT)
library(shinyTree)
library(shinyWidgets)
library(bs4Dash)
library(purrr)
library(shinyFeedback)

color <- read.csv("data/Microplastics_Color.csv")


#Data for embeddings generation via chRoma
items_vectorDB <- readRDS(file = "data/items_vectorDB.rda")
materials_vectorDB <- readRDS(file = "data/materials_vectorDB.rda")
Sys.setenv(OPENAI_API_KEY = readLines("data/openai.txt", warn = FALSE))
creds <- read.csv("data/s3_cred.csv")
Sys.setenv(
  "AWS_ACCESS_KEY_ID" = creds$Access.key.ID,
  "AWS_SECRET_ACCESS_KEY" = creds$Secret.access.key,
  "AWS_DEFAULT_REGION" = "us-east-2"
)

merge_terms <- function(file_paths, materials_vectorDB, items_vectorDB, alias, aliasi, use_cases, prime_unclassifiable) {
  clean_and_prepare_data <- function(dataframe, column, descriptor) {
    if (!column %in% colnames(dataframe)) {
      stop(paste("Column", column, "not found in dataframe"))
    }
    
    if (!"Alias" %in% colnames(descriptor)) {
      stop("Column 'Alias' not found in alias dataframe")
    }
    
    dataframe <- dataframe %>%
      mutate(!!sym(column) := as.character(!!sym(column))) %>%
      mutate(!!sym(column) := cleantext(!!sym(column)))
    
    key <- inner_join(dataframe %>% select(!!sym(column)), descriptor, by = setNames("Alias", column)) %>% distinct()
    
    left <- anti_join(dataframe %>% select(!!sym(column)), descriptor, by = setNames("Alias", column)) %>%
      distinct() %>% drop_na() %>% rename(text = !!sym(column)) %>% as.data.table()
    
    return(list(dataframe = dataframe, key = key, left = left))
  }
  
  process_unknowns <- function(left, vectorDB, descriptor, column, key) {
    if (nrow(left) > 0) {
      new_vDB <- add_collection(metadata = left)
      unknown_key <- query_collection(db = vectorDB, query_embeddings = new_vDB, top_n = 5, type = "dotproduct") %>%
        left_join(vectorDB$metadata, by = c("db_id" = "id")) %>%
        rename(Alias = text) %>%
        left_join(new_vDB$metadata, by = c("query_id" = "id")) %>%
        rename(!!sym(column) := text)
      
      embedding <- unknown_key %>%
        add_column(rank = rep(seq(5), nrow(unknown_key) / 5)) %>%
        select(Alias, !!sym(column), rank) %>%
        pivot_wider(names_from = rank, values_from = Alias) %>%
        rename(!!paste0(column, "_raw") := 1,
               !!paste0(column, "_match_1") := 2,
               !!paste0(column, "_match_2") := 3,
               !!paste0(column, "_match_3") := 4,
               !!paste0(column, "_match_4") := 5,
               !!paste0(column, "_match_5") := 6)
      
      unknown_key <- unknown_key %>% distinct(!!sym(column), .keep_all = TRUE) %>%
        select(Alias, !!sym(column))
      
      readable_col <- ifelse(column == "material", "Material", "Item")
      key_full <- unknown_key %>%
        inner_join(descriptor, by = c("Alias")) %>%
        select(!!sym(column), !!sym(readable_col), readable) %>%
        bind_rows(key)
      
      return(list(key = key_full, embedding = embedding))
    }
    return(NULL)
  }
  
  # Clean and prepare data
  material_data <- clean_and_prepare_data(file_paths, "material", alias)
  morphology_data <- clean_and_prepare_data(file_paths, "morphology", aliasi)
  
  # Process unknown materials
  
  material_results <- process_unknowns(material_data$left, materials_vectorDB, alias, "material", material_data$key)
  if (!is.null(material_results)) {
    material_data$key <- material_results$key 
    material_embedding <- material_results$embedding
  }
  
  # Process unknown morphologies
  morphology_results <- process_unknowns(morphology_data$left, items_vectorDB, aliasi, "morphology", morphology_data$key)
  if (!is.null(morphology_results)) {
    morphology_data$key <- morphology_results$key
    morphology_embedding <- morphology_results$embedding
  }
  
  material_data$dataframe$material <- cleantext(material_data$dataframe$material)
  material_data$dataframe$morphology <- cleantext(material_data$dataframe$morphology)
  
  # Replace old material with merged material
  dataframeclean <- material_data$dataframe %>%
    left_join(material_data$key, by = "material") %>%
    select(-Material) %>%
    rename(material_new = readable) %>%
    left_join(morphology_data$key, by = c("morphology")) %>%
    select(-Item) %>%
    rename(morphology_new = readable) %>%
    select("material", "morphology", "material_new", "morphology_new") %>%
    unique()
  
  if (nrow(material_data$left) > 0) {
    dataframeclean <- dataframeclean %>% left_join(material_embedding, by = c("material" = "material_raw"))
  }
  if (nrow(morphology_data$left) > 0) {
    dataframeclean <- dataframeclean %>% left_join(morphology_embedding, by = c("morphology" = "morphology_raw"))
  }
  
  materialclean <- cleantext(file_paths$material)
  morphologyclean <- cleantext(file_paths$morphology)
  
  dataframe <- file_paths %>%
    rename(material_raw = material, morphology_raw = morphology) %>%
    add_column(material = materialclean, morphology = morphologyclean) %>%
    left_join(dataframeclean, by = c("material", "morphology")) %>%
    select(-c(material, morphology)) %>%
    rename(material = material_new, morphology = morphology_new)
  
  dataframe <- dataframe %>%
    select(morphology_raw, everything()) %>%
    select(material_raw, everything()) %>%
    select(morphology, everything()) %>%
    select(material, everything())
  return(dataframe)
}


merge_data <- function(file_paths, materials_vectorDB, items_vectorDB, alias, aliasi, use_cases, prime_unclassifiable){
  .confidence_interval_width <- function(proportion, sample_size, population_size){
    1.96*sqrt((1/sample_size)*proportion * (1-proportion) * (population_size-sample_size)/(population_size-1))
  }
  sample_size = length(file_paths)
  population_size = 100000 
  
  dataframe <- file_paths %>%
    mutate(proportion = count/sum(count))%>%
    select(material, items, count, proportion) %>%
    mutate(material = as.character(material),
           items = as.character(items), 
           count = as.numeric(count),
           proportion = as.numeric(proportion)) %>%
    group_by(material, items) %>%
    summarise(proportion = sum(proportion),
              count = sum(count)) %>%
    ungroup() %>%
    mutate(proportion = proportion/sum(proportion))
  
  dataframeclean <- dataframe %>%
    mutate(material = cleantext(material), 
           items = cleantext(items)) 
  
  material_key <- inner_join(dataframeclean %>% select(material), 
                             alias, by = c("material" = "Alias")) %>%
    distinct()
  
  materials_left <- anti_join(dataframeclean %>% select(material), 
                              alias, by = c("material" = "Alias")) %>%
    distinct() %>%
    rename(text = material) %>%
    as.data.table()
  
  if(nrow(materials_left) > 0){
    new_material_vDB <- add_collection(metadata = materials_left)
    
    material_key <- query_collection(db = materials_vectorDB, query_embeddings = new_material_vDB, top_n = 1, type = "dotproduct") %>%
      left_join(materials_vectorDB$metadata, by = c("db_id" = "id")) %>%
      rename(Alias = text) %>%
      left_join(new_material_vDB$metadata, by = c("query_id" = "id")) %>%
      rename(material = text) %>%
      inner_join(alias, by = c("Alias")) %>%
      select(material, Material, readable)   %>%
      bind_rows(material_key)
  }
  
  
  #Run for items
  items_key <- inner_join(dataframeclean %>% select(items), 
                          aliasi, by = c("items" = "Alias")) %>%
    distinct()
  
  items_left <- anti_join(dataframeclean %>% select(items), 
                          aliasi, by = c("items" = "Alias")) %>%
    distinct() %>%
    rename(text = items) %>%
    as.data.table()
  
  if(nrow(items_left) > 0){
    new_items_vDB <- add_collection(metadata = items_left)
    
    items_key <- query_collection(db = items_vectorDB, query_embeddings = new_items_vDB, top_n = 1, type = "dotproduct") %>%
      left_join(items_vectorDB$metadata, by = c("db_id" = "id")) %>%
      rename(Alias = text) %>%
      left_join(new_items_vDB$metadata, by = c("query_id" = "id")) %>%
      rename(items = text) %>%
      inner_join(aliasi, by = "Alias") %>%
      select(items, Item, readable)   %>%
      bind_rows(items_key)
  }
  
  
  #Replace old material with merged material
  dataframeclean <- dataframeclean %>% 
    left_join(material_key, by = "material") %>%
    select(-c(material, readable)) %>%
    rename(material = Material)
  
  #Replace old item with merged item
  dataframeclean <- dataframeclean %>% 
    left_join(items_key, by = "items") %>%
    select(-c(items, readable)) %>%
    rename(items = Item)
  
  dataframeclean <- dataframeclean %>%
    left_join(use_cases, by = c("items"="Item"), keep = NULL)
  
  #Combine any new identical terms
  dataframeclean2 <- dataframeclean %>%
    mutate(count = as.numeric(count)) %>%
    group_by(material, items, Use) %>%
    summarise(count = sum(count), 
              proportion = sum(proportion)) %>%
    ungroup() %>% 
    mutate(proportion_width = .confidence_interval_width(proportion = proportion, sample_size = sample_size, population_size = population_size)) %>%
    mutate(min_proportion = proportion - proportion_width, 
           max_proportion = proportion + proportion_width)
  
  dataframeclean2 <- setkey(setDT(dataframeclean2), items) 
  dataframeclean2[aliasi, readable := i.readable]
  dataframeclean2 <- dataframeclean2 %>% select(-items) %>% rename(items = readable)
  
  dataframeclean2 <- setkey(setDT(dataframeclean2), material) 
  dataframeclean2[alias, readable := i.readable]
  dataframeclean2 <- dataframeclean2 %>% select(-material) %>% rename(material = readable)
  
  dataframeclean2 <- setkey(setDT(dataframeclean2), Use) 
  dataframeclean2[prime_unclassifiable, readable := i.readable]
  dataframeclean2 <- dataframeclean2 %>% select(-Use) %>% rename(use = readable) %>% select(-proportion_width)
  
  dataframeclean2$min_proportion <- as.numeric(dataframeclean2$min_proportion)
  dataframeclean2$max_proportion <- as.numeric(dataframeclean2$max_proportion)
  
  for(x in 1:nrow(dataframeclean2)){
    if(dataframeclean2$min_proportion[[x]] < 0){
      dataframeclean2$min_proportion[[x]] <- 0
    }
    if(dataframeclean2$max_proportion[[x]] > 1){
      dataframeclean2$max_proportion[[x]] <- 1
    }
  }
  
  dataframeclean2 <- dataframeclean2[,c("material","items","use","count","proportion","min_proportion","max_proportion")]
  
  return(dataframeclean2)
}

# Function to calculate volume and mass
calculate_volume_and_mass <- function(df, morphology_shape, polymer_density) {
  meas_min <- 0.95
  meas_max <- 1.05
  # Join morphology and density data
  df <- df %>%
    left_join(morphology_shape, by = "morphology") %>%
    left_join(polymer_density, by = "material") %>%
    mutate(material = readable) %>%
    select(-readable)
  
  length_um <- df$length_um
  df <- df %>%
    add_column(L_min = length_um * meas_min,
               L_max = length_um * meas_max)
  
  # Calculate width and height if not provided
  if ("width_um" %in% colnames(df)) {
    width_um <- df$width_um
    df <- df %>%
      mutate(W_mean = ifelse(is.na(width_um), (W_min * length_um + W_max * length_um) / 2, width_um),
             W_min = ifelse(is.na(width_um), W_min * length_um, width_um * meas_min),
             W_max = ifelse(is.na(width_um), W_max * length_um, width_um * meas_max))
  } else {
    df <- df %>%
      mutate(W_mean = (W_min * length_um + W_max * length_um) / 2,
             W_min = W_min * length_um,
             W_max = W_max * length_um)
  }
  
  if ("height_um" %in% colnames(df)) {
    height_um <- df$height_um
    df <- df %>%
      mutate(H_mean = ifelse(is.na(height_um), (H_min * length_um + H_max * length_um) / 2, height_um),
             H_min = ifelse(is.na(height_um), H_min * length_um, height_um * meas_min),
             H_max = ifelse(is.na(height_um), H_max * length_um, height_um * meas_max))
  } else {
    df <- df %>%
      mutate(H_mean = (H_min * length_um + H_max * length_um) / 2,
             H_min = H_min * length_um,
             H_max = H_max * length_um)
  }
  
  # Calculate volume based on morphology
  df <- df %>%
    mutate(volume_min_um_3 = case_when(
      morphology == "fiber" ~ pi * ((W_min + H_min) / 4)^2 * L_min,
      morphology == "sphere" ~ 4/3 * pi * ((W_mean + H_mean + L_min) / 6)^2,
      TRUE ~ L_min * W_min * H_min),
      volume_mean_um_3 = case_when(
        morphology == "fiber" ~ pi * ((W_mean + H_mean) / 4)^2 * length_um,
        morphology == "sphere" ~ 4/3 * pi * ((W_mean + H_mean + length_um) / 6)^2,
        TRUE ~ length_um * W_mean * H_mean),
      volume_max_um_3 = case_when(
        morphology == "fiber" ~ pi * ((W_max + H_max) / 4)^2 * L_max,
        morphology == "sphere" ~ 4/3 * pi * ((W_max + H_max + L_max) / 6)^2,
        TRUE ~ L_max * W_max * H_max))
  
  # Calculate mass
  df <- df %>%
    mutate(min_mass_mg = density_min * volume_min_um_3,
           mean_mass_mg = density_mg_um_3 * volume_mean_um_3,
           max_mass_mg = density_max * volume_max_um_3)
  
  return(df)
}

#dataframe <- read.csv("River_Database_Match_July2024_clean.csv")
#Update polymer
# polymer_class_rename <- function(df) {
#   df <- df %>%
#     mutate(material = ifelse(polymer == 0 & polymer_class != 0, polymer_class, polymer))
#   return(df)
# }
# df <- dataframe
polymer_class_rename <- function(df) {
  df$polymer[df$polymer == ""] <- 0
  df$polymer_class[df$polymer_class == ""] <- 0
  df <- df %>%
    mutate(material = ifelse(polymer == 0 & polymer_class != 0, polymer_class, polymer))
  return(df)
}
#dataframe <- dataframe2
# Main function to convert particle count to mass
particle_count_mass <- function(dataframe, morphology_shape, polymer_density, trash_mass_clean, polymer_avg_decision, morph_weight, sample_weight){
  if("length_um" %in% colnames(dataframe) == TRUE){dataframe$length_um <- as.numeric(dataframe$length_um)
  dataframe$length_um[is.na(dataframe$length_um)] <- 0
  }else{
    dataframe <- dataframe %>% add_column(length_um = NA)
  }
  dataframe$morphology <- as.character(dataframe$morphology)
  dataframe$material <- as.character(dataframe$material)
  if("width_um" %in% colnames(dataframe) == TRUE){dataframe$width_um <- as.numeric(dataframe$width_um)}
  if("height_um" %in% colnames(dataframe) == TRUE){dataframe$height_um <- as.numeric(dataframe$height_um)}
  if("density" %in% colnames(dataframe) == TRUE){dataframe$density <- as.numeric(dataframe$density)}
  if("sample_ID" %in% colnames(dataframe) == TRUE){dataframe$sample_ID <- as.character(dataframe$sample_ID)}
  dataframe <- dataframe %>%
    mutate(material = ifelse(material == "polystyrene" & morphology == "foam", "EPS", material))
  
  dataframe <- left_join(dataframe, polymer_density, by = c("material" = "readable"), copy = F)
  dataframe <- select(dataframe, -material.y)
  
  dataframeclean <- dataframe %>% mutate(material = cleantext(material),
                                         morphology = cleantext(morphology)) 
  
  dataframeclean <- left_join(dataframeclean, morphology_shape, by = "morphology", copy = F)
  
  
  if("width_um" %in% colnames(dataframeclean)){
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "width_um"])) {
        dataframeclean[x, "W_min"] <- as.numeric(dataframeclean[x, "W_min"]) * as.numeric(dataframeclean[x, "length_um"])
        dataframeclean[x, "W_max"] <- as.numeric(dataframeclean[x, "W_max"]) * as.numeric(dataframeclean[x, "length_um"])
      }else{
        dataframeclean[x, "W_min"] <- as.numeric(dataframeclean[x, "width_um"]) * meas_min
        dataframeclean[x, "W_mean"] <- as.numeric(dataframeclean[x, "width_um"])
        dataframeclean[x, "W_max"] <- as.numeric(dataframeclean[x, "width_um"]) * meas_max
      }
    }
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "width_um"])) {
        dataframeclean[x, "W_mean"] <- (as.numeric(dataframeclean[x, "W_min"]) + as.numeric(dataframeclean[x, "W_max"]))/2
      }
    }
  }else{
    for(x in 1:nrow(dataframeclean)){
      dataframeclean[x, "W_min"] <- as.numeric(dataframeclean[x, "W_min"]) * as.numeric(dataframeclean[x, "length_um"])
      dataframeclean[x, "W_max"] <- as.numeric(dataframeclean[x, "W_max"]) * as.numeric(dataframeclean[x, "length_um"])
    }
    for(x in 1:nrow(dataframeclean)){
      dataframeclean[x, "W_mean"] <- (as.numeric(dataframeclean[x, "W_min"]) + as.numeric(dataframeclean[x, "W_max"]))/2
    }
  }
  
  if("height_um" %in% colnames(dataframeclean)){
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "height_um"])) {
        dataframeclean[x, "H_min"] <- as.numeric(dataframeclean[x, "H_min"]) * as.numeric(dataframeclean[x, "length_um"])
        dataframeclean[x, "H_max"] <- as.numeric(dataframeclean[x, "H_max"]) * as.numeric(dataframeclean[x, "length_um"])
      }else{
        dataframeclean[x, "H_min"] <- as.numeric(dataframeclean[x, "height_um"]) * meas_min
        dataframeclean[x, "H_mean"] <- as.numeric(dataframeclean[x, "height_um"])
        dataframeclean[x, "H_max"] <- as.numeric(dataframeclean[x, "height_um"]) * meas_max
      }
    }
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "height_um"])) {
        dataframeclean[x, "H_mean"] <- (as.numeric(dataframeclean[x, "H_min"]) + as.numeric(dataframeclean[x,"H_max"]))/2
      }
    }
  }else{
    for(x in 1:nrow(dataframeclean)){
      dataframeclean[x, "H_min"] <- as.numeric(dataframeclean[x, "H_min"]) * as.numeric(dataframeclean[x, "length_um"])
      dataframeclean[x, "H_max"] <- as.numeric(dataframeclean[x, "H_max"]) * as.numeric(dataframeclean[x, "length_um"])
    }
    for(x in 1:nrow(dataframeclean)){
      dataframeclean[x, "H_mean"] <- (as.numeric(dataframeclean[x, "H_min"]) + as.numeric(dataframeclean[x,"H_max"]))/2
    }
    
    for(x in 1:nrow(dataframeclean)){
      print(x)
      if(!is.na(dataframeclean[x, "H_mean"]) && dataframeclean[x, "H_mean"] > dataframeclean[x, "W_mean"]){
        dataframeclean[x, "H_max"] <- as.numeric(dataframeclean[x, "W_mean"])
        dataframeclean[x, "H_mean"] <- (as.numeric(dataframeclean[x, "H_min"]) + as.numeric(dataframeclean[x,"H_max"]))/2
      }
      
    }
    
  }
  
  dataframeclean <- dataframeclean %>%
    mutate(L_min = as.numeric(L_min) * as.numeric(length_um),
           L_mean = as.numeric(length_um),
           L_max = as.numeric(L_max) * as.numeric(length_um))
  
  # Calculate volume based on morphology
  dataframeclean <- dataframeclean %>%
    mutate(volume_min_um_3 = case_when(
      morphology == "fiber" ~ pi * ((W_min + H_min) / 4)^2 * L_min,
      morphology == "sphere" ~ 4/3 * pi * ((W_min + H_min + L_min) / 6)^3,
      TRUE ~ L_min * W_min * H_min),
      volume_mean_um_3 = case_when(
        morphology == "fiber" ~ pi * ((W_mean + H_mean) / 4)^2 * L_mean,
        morphology == "sphere" ~ 4/3 * pi * ((W_mean + H_mean + L_mean) / 6)^3,
        TRUE ~ L_mean * W_mean * H_mean),
      volume_max_um_3 = case_when(
        morphology == "fiber" ~ pi * ((W_max + H_max) / 4)^2 * L_max,
        morphology == "sphere" ~ 4/3 * pi * ((W_max + H_max + L_max) / 6)^3,
        TRUE ~ L_max * W_max * H_max))
  
  
  
  if("density" %in% colnames(dataframeclean)){
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "density"]) && is.na(dataframeclean[x, "density_max"])) {
        dataframeclean[x, "min_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x, "volume_min_um_3"])
        dataframeclean[x, "mean_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x,"volume_mean_um_3"])
        dataframeclean[x, "max_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x, "volume_max_um_3"])
      }else if(is.na(dataframeclean[x, "density"]) && !is.na(dataframeclean[x, "density_max"])){
        dataframeclean[x, "min_mass_mg"] <- as.numeric(dataframeclean[x, "density_min"]) * as.numeric(dataframeclean[x, "volume_min_um_3"])
        dataframeclean[x, "mean_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x,"volume_mean_um_3"])
        dataframeclean[x, "max_mass_mg"] <- as.numeric(dataframeclean[x, "density_max"]) * as.numeric(dataframeclean[x, "volume_max_um_3"])
      }else{
        dataframeclean[x, "min_mass_mg"] <- as.numeric(dataframeclean[x, "density"]) * as.numeric(dataframeclean[x, "volume_min_um_3"])
        dataframeclean[x, "mean_mass_mg"] <- as.numeric(dataframeclean[x, "density"]) * as.numeric(dataframeclean[x,"volume_mean_um_3"])
        dataframeclean[x, "max_mass_mg"] <- as.numeric(dataframeclean[x, "density"]) * as.numeric(dataframeclean[x, "volume_max_um_3"])
      }
    }
  }else{
    for(x in 1:nrow(dataframeclean)){
      if(is.na(dataframeclean[x, "density_max"]) && is.na(dataframeclean[x, "density_min"])){
        dataframeclean[x, "min_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x, "volume_min_um_3"])
        dataframeclean[x, "mean_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x,"volume_mean_um_3"])
        dataframeclean[x, "max_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x, "volume_max_um_3"])
      }else{
        dataframeclean[x, "min_mass_mg"] <- as.numeric(dataframeclean[x, "density_min"]) * as.numeric(dataframeclean[x, "volume_min_um_3"])
        dataframeclean[x, "mean_mass_mg"] <- as.numeric(dataframeclean[x, "density_mg_um_3"]) * as.numeric(dataframeclean[x,"volume_mean_um_3"])
        dataframeclean[x, "max_mass_mg"] <- as.numeric(dataframeclean[x, "density_max"]) * as.numeric(dataframeclean[x, "volume_max_um_3"])
      }
    }
  }
  
  dataframeclean_particles <- data.frame(dataframeclean)
  
  dataframeclean_particles <- dataframeclean_particles %>%
    left_join(trash_mass_clean, by = c("morphology" = "items",
                                       "material" = "material"))
  dataframeclean_particles$weight_estimate_g <- as.numeric(dataframeclean_particles$weight_estimate_g)
  for(x in 1:nrow(dataframeclean_particles)){
    if(! is.na(dataframeclean_particles[x, "weight_estimate_g"])){
      dataframeclean_particles[x, "mean_mass_mg"] <- (dataframeclean_particles[x, "weight_estimate_g"]*1000)
    }
  }
  
  if(polymer_avg_decision == T){
    dataframeclean_particles$density_mg_um_3 <- as.numeric(dataframeclean_particles$density_mg_um_3)
    dataframeclean_particles$volume_mean_um_3 <- as.numeric(dataframeclean_particles$volume_mean_um_3)
    dataframeclean_particles$mean_mass_mg <- as.numeric(dataframeclean_particles$mean_mass_mg)
    dataframeclean_particles$density_max <- as.numeric(dataframeclean_particles$density_max)
    dataframeclean_particles$density_min <- as.numeric(dataframeclean_particles$density_min)
    dataframeclean_particles_avg <- dataframeclean_particles
    for(x in 1:nrow(dataframeclean_particles)){
      if(is.na(dataframeclean_particles$mean_mass_mg[[x]]) && ! is.na(dataframeclean_particles$volume_mean_um_3[[x]])){
        if(morph_weight == T && sample_weight == T && "sample_ID" %in% colnames(dataframeclean_particles) == TRUE){
          dataframeclean_particles_avg <- dataframeclean_particles_avg %>% filter(sample_ID == dataframeclean_particles$sample_ID[[x]])
          density <- mean(dataframeclean_particles_avg$density_mg_um_3, na.rm = T)
          dataframeclean_particles$density_max[[x]] <- mean(dataframeclean_particles_avg$density_max, na.rm = T)
          dataframeclean_particles$density_min[[x]] <- mean(dataframeclean_particles_avg$density_min, na.rm = T)
          dataframeclean_particles$density_mg_um_3[[x]] <- density
          dataframeclean_particles$mean_mass_mg[[x]] <- (dataframeclean_particles$volume_mean_um_3[[x]])*density
          upper_diff <- (dataframeclean_particles_avg$density_max) - (dataframeclean_particles_avg$density_mg_um_3)
          lower_diff <- (dataframeclean_particles_avg$density_mg_um_3) - (dataframeclean_particles_avg$density_min)
        }else if(morph_weight == T || morph_weight == T && sample_weight == T && !("sample_ID" %in% colnames(dataframeclean_particles))){
          dataframeclean_particles_avg <- dataframeclean_particles %>% filter(morphology == dataframeclean_particles$morphology[[x]])
          density <- mean(dataframeclean_particles_avg$density_mg_um_3, na.rm = T)
          dataframeclean_particles$density_max[[x]] <- mean(dataframeclean_particles_avg$density_max, na.rm = T)
          dataframeclean_particles$density_min[[x]] <- mean(dataframeclean_particles_avg$density_min, na.rm = T)
          dataframeclean_particles$density_mg_um_3[[x]] <- density
          dataframeclean_particles$mean_mass_mg[[x]] <- (dataframeclean_particles$volume_mean_um_3[[x]])*density
          upper_diff <- (dataframeclean_particles_avg$density_max) - (dataframeclean_particles_avg$density_mg_um_3)
          lower_diff <- (dataframeclean_particles_avg$density_mg_um_3) - (dataframeclean_particles_avg$density_min)
        }else if(morph_weight == T && "sample_ID" %in% colnames(dataframeclean_particles) == TRUE){
          dataframeclean_particles_avg <- dataframeclean_particles %>% filter(morphology == dataframeclean_particles$morphology[[x]])
          dataframeclean_particles_avg <- dataframeclean_particles_avg %>% filter(sample_ID == dataframeclean_particles$sample_ID[[x]])
          density <- mean(dataframeclean_particles_avg$density_mg_um_3, na.rm = T)
          dataframeclean_particles$density_max[[x]] <- mean(dataframeclean_particles_avg$density_max, na.rm = T)
          dataframeclean_particles$density_min[[x]] <- mean(dataframeclean_particles_avg$density_min, na.rm = T)
          dataframeclean_particles$density_mg_um_3[[x]] <- density
          dataframeclean_particles$mean_mass_mg[[x]] <- (dataframeclean_particles$volume_mean_um_3[[x]])*density
          upper_diff <- (dataframeclean_particles_avg$density_max) - (dataframeclean_particles_avg$density_mg_um_3)
          lower_diff <- (dataframeclean_particles_avg$density_mg_um_3) - (dataframeclean_particles_avg$density_min)
        }else{
          density <- mean(dataframeclean_particles$density_mg_um_3, na.rm = T)
          dataframeclean_particles$mean_mass_mg[[x]] <- (dataframeclean_particles$volume_mean_um_3[[x]])*density
          upper_diff <- (dataframeclean_particles$density_max) - (dataframeclean_particles$density_mg_um_3)
          lower_diff <- (dataframeclean_particles$density_mg_um_3) - (dataframeclean_particles$density_min)
        }
        upper_conf <- mean(upper_diff, na.rm = T) 
        lower_conf <- mean(lower_diff, na.rm = T)
        dataframeclean_particles$max_mass_mg[[x]] <- dataframeclean_particles$mean_mass_mg[[x]] + upper_conf
        dataframeclean_particles$min_mass_mg[[x]] <- dataframeclean_particles$mean_mass_mg[[x]] + lower_conf
      }
    }
  }
  
  dataframeclean_particles <- dataframeclean_particles %>% 
    select(-c(weight_estimate_g))
  
  if("material_raw" %in% colnames(dataframe) && "morphology_raw" %in% colnames(dataframe)){
    dataframeclean_particles <- dataframeclean_particles %>% 
      select(-c(material, morphology))
    
    dataframe <- dataframe %>% select(material, morphology, material_raw, morphology_raw) %>%
      distinct()
    
    dataframeclean_particles <- dataframeclean_particles %>% 
      left_join(dataframe, by = c("material_raw", "morphology_raw"))
  }
  
  
  return(dataframeclean_particles)
}


concentration_count_mass <- function(dataframe, morphology_shape, polymer_density, corrected_DF, trash_mass_clean){
  dataframe$concentration_particle_vol <- as.numeric(dataframe$concentration_particle_vol)
  dataframe$sample_ID <- as.character(dataframe$sample_ID)
  if("morphology" %in% colnames(dataframe) == TRUE){dataframe$morphology <- as.character(dataframe$morphology)}
  if("material" %in% colnames(dataframe) == TRUE){dataframe$material <- as.character(dataframe$material)}
  if("density" %in% colnames(dataframe) == TRUE){dataframe$density <- as.numeric(dataframe$density)}
  if("error_SD" %in% colnames(dataframe) == TRUE){dataframe$error_SD <- as.numeric(dataframe$error_SD)}
  if("error_upper" %in% colnames(dataframe) == TRUE){dataframe$error_upper <- as.numeric(dataframe$error_upper)}
  if("error_lower" %in% colnames(dataframe) == TRUE){dataframe$error_lower <- as.numeric(dataframe$error_lower)}
  dataframeclean <- dataframe
  
  corrected_DF <- corrected_DF %>% select(sample_ID, alpha, alpha_upper, alpha_lower)
  corrected_DF$alpha <- as.numeric(corrected_DF$alpha)
  corrected_DF$alpha_upper <- as.numeric(corrected_DF$alpha_upper)
  corrected_DF$alpha_lower <- as.numeric(corrected_DF$alpha_lower)
  dataframeclean <- dataframeclean %>% left_join(corrected_DF, by = "sample_ID")
  
  dataframeclean$size_min <- as.numeric(dataframeclean$size_min)
  dataframeclean$size_max <- as.numeric(dataframeclean$size_max)
  
  unique_IDs <- dataframeclean %>% select(sample_ID) %>% unique()
  concentration_mass <- data.frame()
  concentration_mass_trash <- data.frame()
  for(x in 1:nrow(unique_IDs)){
    #x = 1
    print(x)
    ID <- unique_IDs[x, "sample_ID"]
    data_ID <- dataframeclean %>% filter(sample_ID == ID)
    data_ID <- discard(data_ID, ~all(is.na(.)))
    
    if("size_min" %in% colnames(data_ID) == TRUE && "size_max" %in% colnames(data_ID) == TRUE){
      concentration_ID <- data_ID %>% select(sample_ID, concentration_particle_vol, size_min, size_max, alpha) %>% filter(!is.na(size_min)) %>% filter(!is.na(size_max)) %>% unique()
    }else{
      concentration_ID <- data_ID %>% select(sample_ID, concentration_particle_vol) %>% unique()
    }
    
    if("morphology" %in% colnames(data_ID) == TRUE && "morphology_percent" %in% colnames(data_ID) == TRUE){
      data_ID$morphology_percent <- as.numeric(data_ID$morphology_percent)
      morph_ID <- data_ID %>% select(morphology, morphology_percent) %>% mutate(morphology_percent = morphology_percent/100) %>% unique() %>% filter_all(all_vars(!is.na(.)))
    }else{
      morphology <- c("fiber", "fragment", "sphere", "film", "foam")
      morphology_percent <- c(0.485/0.95, 0.31/0.95, 0.065/0.95, 0.055/0.95, 0.035/0.95)
      morph_ID <- data.frame(morphology, morphology_percent)
    }
    
    if("material" %in% colnames(data_ID) == TRUE && "material_percent" %in% colnames(data_ID) == TRUE){
      data_ID$material_percent <- as.numeric(data_ID$material_percent)
      material_ID <- data_ID %>% select(material, material_percent) %>% mutate(material_percent = material_percent/100) %>% unique() %>% filter_all(all_vars(!is.na(.)))
    }else{
      material <- c("PE", "PP", "PET", "PA", "PS", "PVC", "PVA")
      material_percent <- c(0.25/0.845, 0.145/0.845, 0.165/0.845, 0.12/0.845, 0.085/0.845, 0.02/0.845, 0.06/0.845)
      material_ID <- data.frame(material, material_percent)
    }
    
    # Function to generate power law distributed random numbers
    generate_power_law <- function(n, xmin, xmax, alpha) {
      u <- runif(n) # Generate uniform random numbers
      x <- ((xmax^(1 - alpha) - xmin^(1 - alpha)) * u + xmin^(1 - alpha))^(1 / (1 - alpha))
      return(x)
    }
    
    concentration_ID <- concentration_ID %>% add_column(min_concentration_um3_vol = NA, mean_concentration_um3_vol = NA, max_concentration_um3_vol = NA, min_concentration_mg_vol = NA, mean_concentration_mg_vol = NA, max_concentration_mg_vol = NA)
    #y = 1
    for(y in 1:nrow(concentration_ID)){
      if("size_min" %in% colnames(data_ID) == TRUE && "size_max" %in% colnames(data_ID) == TRUE){
        # Parameters
        total <- as.numeric(concentration_ID[y, "concentration_particle_vol"])
        ifelse(total <10 || total > 1000, n <- 1000, n <- total)
        xmin <- concentration_ID[y, "size_min"]
        xmax <- concentration_ID[y, "size_max"]
        alpha <- concentration_ID[y, "alpha"]
        
        # Generate the data
        set.seed(123)
        particle_lengths <- generate_power_law(n, xmin, xmax, alpha)
        
        # Calculate the number of particles for each morphology and material
        morph_ID <- morph_ID %>% mutate(count = round(n * morphology_percent))
        material_ID <- material_ID %>% mutate(count = round(n * material_percent))
        
        # Create a vector of morphologies and materials according to the calculated counts
        morphologies <- rep(morph_ID$morphology, morph_ID$count)
        if(length(morphologies) < n){morphologies <- append(morphologies, sample(morphologies, n-length(morphologies)))}
        if(length(morphologies) > n){length(morphologies) <- n}
        materials <- rep(material_ID$material, material_ID$count)
        if(length(materials) < n){materials <- append(materials, sample(materials, n-length(materials)))}
        if(length(materials) > n){length(materials) <- n}
        # Randomly shuffle the morphologies and materials to mix them up
        set.seed(123) # For reproducibility
        morphologies <- sample(morphologies)
        set.seed(123) # For reproducibility
        materials <- sample(materials)
        
      }else{
        #For trash data, don't perform the inverse power dist, don't mix up material and morphology types
        particle_lengths = NA
        # Calculate the number of particles for each morphology and material
        n <- as.numeric(concentration_ID[y, "concentration_particle_vol"])
        morph_ID <- morph_ID %>% mutate(count = round(n * morphology_percent))
        material_ID <- material_ID %>% mutate(count = round(n * material_percent))
        
        # Create a vector of morphologies and materials according to the calculated counts
        morphologies <- rep(morph_ID$morphology, morph_ID$count)
        materials <- rep(material_ID$material, material_ID$count)
      }
      
      # Combine particle lengths with morphologies
      particle_data <- data.frame(
        length_um = particle_lengths,
        morphology = morphologies,
        material = materials
      )
      
      
      if("density" %in% colnames(data_ID) == T){
        density_ID <- data_ID %>% select(material, density) %>% unique() %>% filter_all(all_vars(!is.na(.)))
        particle_data <- particle_data %>% left_join(density_ID, by = "material")
      }
      
      particle_mass <- particle_count_mass(dataframe = particle_data, morphology_shape = morphology_shape, polymer_density = polymer_density, trash_mass_clean = trash_mass_clean, 
                          polymer_avg_decision = F, morph_weight = F, sample_weight = F)
      particle_mass <- particle_mass %>% filter(!is.na(mean_mass_mg))
      
      concentration_ID[y, "min_concentration_um3_vol"] <- sum(particle_mass$volume_min_um_3)
      concentration_ID[y, "mean_concentration_um3_vol"] <- sum(particle_mass$volume_mean_um_3)
      concentration_ID[y, "max_concentration_um3_vol"] <- sum(particle_mass$volume_max_um_3)
      concentration_ID[y, "min_concentration_mg_vol"] <- sum(particle_mass$min_mass_mg)
      concentration_ID[y, "mean_concentration_mg_vol"] <- sum(particle_mass$mean_mass_mg)
      concentration_ID[y, "max_concentration_mg_vol"] <- sum(particle_mass$max_mass_mg)
      
      if(total <1 || total > 1000){concentration_ID <- concentration_ID %>% mutate(min_concentration_um3_vol = (min_concentration_um3_vol/n)*total,
                                                      mean_concentration_um3_vol = (mean_concentration_um3_vol/n)*total,
                                                      max_concentration_um3_vol = (max_concentration_um3_vol/n)*total,
                                                      min_concentration_mg_vol = (min_concentration_mg_vol/n)*total,
                                                      mean_concentration_mg_vol = (mean_concentration_mg_vol/n)*total,
                                                      max_concentration_mg_vol = (max_concentration_mg_vol/n)*total)}
      
    }
    
    if("size_min" %in% colnames(concentration_ID) == FALSE && "size_max" %in% colnames(concentration_ID) == FALSE){
      concentration_mass_trash <- rbind(concentration_mass_trash, concentration_ID) %>% discard(~all(is.na(.)))
    }else{
      concentration_mass <- rbind(concentration_mass, concentration_ID)
      
    }
  }
  
  dataframe <- dataframe %>% left_join(concentration_mass, by = c("sample_ID", "concentration_particle_vol", "size_min", "size_max")) %>% select(-alpha)
  if(nrow(concentration_mass_trash) > 0){
    dataframe <- dataframe %>% 
      left_join(concentration_mass_trash, by = c("sample_ID", "concentration_particle_vol")) %>% 
      mutate(mean_concentration_mg_vol = coalesce(mean_concentration_mg_vol.y, mean_concentration_mg_vol.x)) %>%
      select(-c(mean_concentration_mg_vol.y, mean_concentration_mg_vol.x))
  }
  
  return(dataframe)
}

correctionFactor_conc <- function(dataframe, alpha_vals, metric, corrected_min, corrected_max){
  dataframe$concentration_particle_vol <- as.numeric(dataframe$concentration_particle_vol)
  dataframe$size_min <- as.numeric(dataframe$size_min)
  dataframe$size_max <- as.numeric(dataframe$size_max)
  dataframe$sample_ID <- as.character(dataframe$sample_ID)
  if("study_media" %in% colnames(dataframe) == TRUE){dataframe$study_media <- as.character(dataframe$study_media)}
  if("known_alpha" %in% colnames(dataframe) == TRUE){dataframe$known_alpha <- as.numeric(dataframe$known_alpha)}
  if("sample_ID"  %in% colnames(dataframe) == TRUE){dataframe$sample_ID <- as.character(dataframe$sample_ID)}
  if("error_SD" %in% colnames(dataframe) == TRUE){dataframe$error_SD <- as.numeric(dataframe$error_SD)}
  if("error_upper" %in% colnames(dataframe) == TRUE){dataframe$error_upper <- as.numeric(dataframe$error_upper)}
  if("error_lower" %in% colnames(dataframe) == TRUE){dataframe$error_lower <- as.numeric(dataframe$error_lower)}
  dataframeclean <- dataframe
  if("study_media" %in% colnames(dataframe) == TRUE){
    if(metric == "length (um)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "length", "length_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "length",
             "sd" = "length_sd")}
    if(metric == "width (um)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "width", "width_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "length",
             "sd" = "width_sd")}
    if(metric == "mass (ug)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "mass", "mass_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "mass",
             "sd" = "mass_sd")}
    if(metric == "volume (um3)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "volume", "volume_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "volume",
             "sd" = "volume_sd")}
    if(metric == "surface area (um2)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "surface_area", "surface_area_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "surface_area",
             "sd" = "surface_area_sd")}
    if(metric == "specific surface area (g/m2)"){dataframeclean <- merge(x = dataframeclean, y = alpha_vals[ , c("study_media", "specific_surface_area", "specific_surface_area_sd")], by = "study_media", all.x=TRUE)
    dataframeclean <- dataframeclean %>%
      rename("alpha" = "specific_surface_area",
             "sd" = "specific_surface_area_sd")}
    alpha_lower <- (dataframeclean$alpha) - (dataframeclean$sd)
    alpha_upper <- (dataframeclean$alpha) + (dataframeclean$sd)
    dataframeclean <- dataframeclean %>% 
      add_column(alpha_lower = alpha_lower,
                 alpha_upper = alpha_upper) %>%
      select(-study_media, -sd)
  }else{
    dataframeclean <- dataframeclean %>% 
      add_column(alpha = NA,
                 alpha_lower = NA,
                 alpha_upper = NA
      )
  }
  
  if("sample_ID" %in% colnames(dataframeclean) == TRUE){
    unique_bins <- dataframeclean %>% count(sample_ID, size_min, size_max)
    bin_numbers <- unique_bins %>% count(sample_ID)
    alpha_bins <- bin_numbers %>% filter(n >= 5)
    if(nrow(alpha_bins) >= 1){
      unique_alpha <- data.frame()
      for(x in 1:nrow(alpha_bins)){
        x = 1
        sample_name <- alpha_bins$sample_ID[[x]]
        sample_bins <- dataframeclean %>% filter(sample_ID == sample_name)
        midpoint <- (as.numeric(sample_bins$size_min) + as.numeric(sample_bins$size_max))/2
        sample_bins <- sample_bins %>%
          add_column(midpoint = midpoint,
                     alpha_calc_lower = NA,
                     alpha_calc_upper = NA) %>%
          select(sample_ID, midpoint, concentration_particle_vol, alpha, alpha_calc_lower, alpha_calc_upper)
        sample_bins$log_size <- log10(sample_bins$midpoint)
        sample_bins$concentration_particle_vol <- as.numeric(sample_bins$concentration_particle_vol)
        sample_bins$log_abundance <- log10(sample_bins$concentration_particle_vol)
        r1model <- lm(log_abundance ~ log_size, data = sample_bins)
        alpha <- -(as.numeric(coef(r1model)[2]))
        sample_bins$alpha <- alpha
        r1model_sum <- summary.lm(r1model)
        error <- as.numeric(coef(r1model_sum)[4])
        sample_bins$alpha_calc_lower <- alpha - error
        sample_bins$alpha_calc_upper <- alpha + error
        sample_bins <- sample_bins %>% select(sample_ID, alpha, alpha_calc_lower, alpha_calc_upper)
        sample_bins <- unique(sample_bins)
        unique_alpha <- rbind(unique_alpha, sample_bins)
        unique_alpha <- unique_alpha %>% rename(alpha_calc = alpha)
      }
      dataframeclean <- dataframeclean %>% left_join(unique_alpha, by = "sample_ID")
      for(x in 1:nrow(dataframeclean)){
        if(! is.na(dataframeclean[x, "alpha_calc"])){
          dataframeclean[x, "alpha"] <- dataframeclean[x, "alpha_calc"]
          dataframeclean[x, "alpha_lower"] <- dataframeclean[x, "alpha_calc_lower"]
          dataframeclean[x, "alpha_upper"] <- dataframeclean[x, "alpha_calc_upper"]
        }
      }
      dataframeclean <- dataframeclean %>% select(-c(alpha_calc, alpha_calc_lower, alpha_calc_upper))
    }
  }
  
  if("known_alpha" %in% colnames(dataframeclean) == TRUE){
    for(x in 1:nrow(dataframeclean)){
      if(! is.na(dataframeclean[x, "known_alpha"])){
        known_alpha_val <- as.numeric(dataframeclean[x, "known_alpha"])
        dataframeclean[x, "alpha"] <- (dataframeclean[x, "known_alpha"])
        if(metric == "length (um)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * length_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * length_coef)}
        if(metric == "width (um)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * width_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * width_coef)}
        if(metric == "mass (ug)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * mass_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * mass_coef)}
        if(metric == "volume (um3)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * volume_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * volume_coef)}
        if(metric == "surface area (um2)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * surface_area_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * surface_area_coef)}
        if(metric == "specific surface area (g/m2)"){dataframeclean[x, "alpha_lower"] <- known_alpha_val - (known_alpha_val * specific_surface_area_coef)
        dataframeclean[x, "alpha_upper"] <- known_alpha_val + (known_alpha_val * specific_surface_area_coef)}
      }
    }
    dataframeclean <- dataframeclean %>% select(-known_alpha)
  }
  
  for(x in 1:nrow(dataframeclean)){
    if(is.na(dataframeclean[x, "alpha"])){
      dataframeclean[x, "alpha"] <- 1.6
      dataframeclean[x, "alpha_lower"] <- 1.1 #1.6 - 0.5
      dataframeclean[x, "alpha_upper"] <- 2.1 #1.6 + 0.5
    }
  }
  
  dataframeclean <- dataframeclean %>%
    add_column(concentration_lower = NA,
               concentration_upper = NA)
  
  if("sample_ID" %in% colnames(dataframeclean) == TRUE){
    unique_bins <- dataframeclean %>% count(sample_ID, size_min, size_max)
    bin_numbers <- unique_bins %>% count(sample_ID)
    bins_add <- bin_numbers %>% filter(n >= 2)
    
    dataframeclean_ <- dataframeclean %>% distinct(sample_ID, .keep_all = T)
    
    if(nrow(bins_add) >= 1){
      sample_add <- data.frame()
      for(x in 1:nrow(bins_add)){
        #x = 1
        sample_name <- bins_add$sample_ID[[x]]
        sample_bins <- dataframeclean %>% filter(sample_ID == sample_name)
        sample_bins$size_min <- as.numeric(sample_bins$size_min)
        sample_bins$size_max <- as.numeric(sample_bins$size_max)
        sample_bins$concentration_particle_vol <- as.numeric(sample_bins$concentration_particle_vol)
        row <- which(grepl(sample_name, dataframeclean_$sample_ID))
        dataframeclean_$size_min[[row]] <- min(sample_bins$size_min)
        dataframeclean_$size_max[[row]] <- max(sample_bins$size_max)
        dataframeclean_$concentration_particle_vol[[row]] <- sum(sample_bins$concentration_particle_vol)
        if("error_SD" %in% colnames(dataframeclean) == TRUE && !(is.na(sample_bins$error_SD[[1]]))){
            sample_bins$error_SD <- as.numeric(sample_bins$error_SD)
            dataframeclean_$concentration_particle_vol <- as.numeric(dataframeclean_$concentration_particle_vol)
            dataframeclean_$concentration_upper[[row]] <- dataframeclean_$concentration_particle_vol[[row]] + (max(sample_bins$error_SD))
            dataframeclean_$concentration_lower[[row]] <- dataframeclean_$concentration_particle_vol[[row]] - (max(sample_bins$error_SD))
        }
        
        if("error_lower" %in% colnames(dataframeclean) == TRUE && "error_upper" %in% colnames(dataframeclean) == TRUE && !(is.na(sample_bins$error_upper[[1]])) && !(is.na(sample_bins$error_lower[[1]]))){
          sample_bins$error_upper <- as.numeric(sample_bins$error_upper)
          sample_bins$error_lower <- as.numeric(sample_bins$error_lower)
          dataframeclean_$concentration_upper[[row]] <- sum(sample_bins$error_upper)
          dataframeclean_$concentration_lower[[row]] <- sum(sample_bins$error_lower)
        }
      }
    }
    
    dataframeclean <- dataframeclean_
    ID <- dataframe %>% select(sample_ID)
    ID_clean <- mutate_all(ID, cleantext) 
    unique_ID <- unique(ID)
    unique_ID_clean <- unique(ID_clean)
    ID_raw_clean <- data.table(unique_ID = unique_ID, unique_ID_clean = unique_ID_clean) %>%
      rename(unique_ID = unique_ID.sample_ID,
             unique_ID_clean = unique_ID_clean.sample_ID)
    
    unique_ID_clean <- as.list(unique_ID_clean$sample_ID)
    dataframeclean <- dataframeclean %>% arrange(factor(sample_ID, levels = unique_ID_clean))
    dataframeclean$alpha <- as.numeric(dataframeclean$alpha)
  }
  
  for(x in 1:nrow(dataframeclean)){
    if(is.na(dataframeclean$concentration_upper[[x]])){
      dataframeclean$concentration_upper[[x]] <- dataframeclean$concentration_particle_vol[[x]]
    }
    if(is.na(dataframeclean$concentration_lower[[x]])){
      dataframeclean$concentration_lower[[x]] <- dataframeclean$concentration_particle_vol[[x]]
    }
  }
  
  dataframeclean <- dataframeclean %>%
    add_column(correction_factor = NA,
               correction_factor_lower = NA,
               correction_factor_upper = NA,
               corrected_concentration = NA,
               corrected_concentration_lower = NA,
               corrected_concentration_upper = NA)
  
  #Extrapolated parameters
  x1D_set = as.numeric(corrected_min) #lower limit default extrapolated range is 1 um
  x2D_set = as.numeric(corrected_max) #upper limit default extrapolated range is 5 mm
  
  for(x in 1:nrow(dataframeclean)) {
    #mean alpha
    x1M_set = as.numeric(dataframeclean$size_min[[x]])
    x2M_set = as.numeric(dataframeclean$size_max[[x]])
    alpha = as.numeric(dataframeclean$alpha[[x]])
    if(alpha == 1){
      CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
      CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
      CF <- (as.numeric(CFL) + as.numeric(CFU))/2
    }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                      x2M = x2M_set, #upper measured length
                      x1D = x1D_set, #default lower size range
                      x2D = x2D_set,  #default upper size range
                      a = alpha #alpha
                      )
    }
    
    dataframeclean$correction_factor[[x]] <- as.numeric(CF)
    dataframeclean$corrected_concentration[[x]] <- as.numeric(dataframeclean$correction_factor[[x]]) * as.numeric(dataframeclean$concentration_particle_vol[[x]])
    
    #min alpha
    alpha = as.numeric(dataframeclean$alpha_lower[[x]])
    if(alpha == 1){
      CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
      CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
      CF <- (as.numeric(CFL) + as.numeric(CFU))/2
    }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                      x2M = x2M_set, #upper measured length
                      x1D = x1D_set, #default lower size range
                      x2D = x2D_set,  #default upper size range
                      a = alpha #alpha
    )
    }
    dataframeclean$correction_factor_lower[[x]] <- as.numeric(CF)
    dataframeclean$corrected_concentration_lower[[x]] <- as.numeric(dataframeclean$correction_factor_lower[[x]]) * as.numeric(dataframeclean$concentration_lower[[x]])
    
    #max alpha
    alpha = as.numeric(dataframeclean$alpha_upper[[x]])
    if(alpha == 1){
      CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
      CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
      CF <- (as.numeric(CFL) + as.numeric(CFU))/2
    }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                      x2M = x2M_set, #upper measured length
                      x1D = x1D_set, #default lower size range
                      x2D = x2D_set,  #default upper size range
                      a = alpha #alpha
    )
    }
    dataframeclean$correction_factor_upper[[x]] <- as.numeric(CF)
    dataframeclean$corrected_concentration_upper[[x]] <- as.numeric(dataframeclean$correction_factor_upper[[x]]) * as.numeric(dataframeclean$concentration_upper[[x]])
  }
  dataframeclean <- dataframeclean %>% 
    select(-c(concentration_lower, concentration_upper)) 
  return(dataframeclean)
}

correctionFactor_particle <- function(dataframe, corrected_min, corrected_max, binning_type, bin_number){
  
  dataframe$length_um <- as.numeric(dataframe$length_um)
  dataframe$sample_ID <- as.character(dataframe$sample_ID)
  if("sample_volume" %in% colnames(dataframe) == TRUE){
    dataframe$sample_volume <- as.numeric(dataframe$sample_volume)}
  
  unique_sample_ID <- unique(dataframe$sample_ID)
  unique_sample_ID <- as.data.frame(unique_sample_ID) %>%
    rename(sample_ID = unique_sample_ID)
  unique_alpha <- data.frame()
  
  for(x in 1:nrow(unique_sample_ID)){
    subset <- filter(dataframe, sample_ID == unique_sample_ID$sample_ID[[x]])
    subset <- subset %>% filter(!is.na(length_um))
    sample_size <- as.numeric(nrow(subset))
    if(sample_size != 0){
      int <- classify_intervals(subset$length_um, n = bin_number, style = binning_type)
      midpoints <- get_midpoints(x = int)
      midpoints <- as.data.frame(midpoints)
      freq <- midpoints %>% count(midpoints) %>%
        rename(length_um = midpoints,
               abundance = n)
      freq$log_size <- log10(freq$length_um)
      freq$log_abundance <- log10(freq$abundance)
      r1model <- lm(log_abundance ~ log_size, data = freq)
      alpha <- -(as.numeric(coef(r1model)[2]))
      r1model_sum <- summary.lm(r1model)
      error <- as.numeric(coef(r1model_sum)[4])
      alpha_lower <- alpha - error
      alpha_upper <- alpha + error
      subset_ <- subset %>% add_column(alpha = alpha,
                                       alpha_lower = alpha_lower,
                                       alpha_upper = alpha_upper) %>% 
        select(sample_ID, alpha, alpha_lower, alpha_upper)
      subset_ <- unique(subset_)
      unique_alpha <- rbind(unique_alpha, subset_)
    }
    
  }
  unique_alpha <- unique_alpha %>% unique()
  dataframe <- left_join(dataframe, unique_alpha, by = "sample_ID")
  
  if("sample_volume" %in% colnames(dataframe) == TRUE){
    dataframe_ <- data.table()
    
    for(x in 1:nrow(unique_sample_ID)){
      subset <- filter(dataframe, sample_ID == unique_sample_ID$sample_ID[[x]])
      particle_number <- subset %>% count(sample_volume) %>%
        rename(particle_number = n)
      concentration <- (as.numeric(particle_number$particle_number))/(as.numeric(particle_number$sample_volume))
      subset <- subset %>% add_column(concentration = concentration,
                                      correction_factor = NA,
                                      correction_factor_lower = NA,
                                      correction_factor_upper = NA,
                                      corrected_concentration = NA,
                                      corrected_concentration_lower = NA,
                                      corrected_concentration_upper = NA)
      #Extrapolated parameters
      x1D_set = as.numeric(corrected_min) #lower limit default extrapolated range is 1 um
      x2D_set = as.numeric(corrected_max) #upper limit default extrapolated range is 5 mm
      
      #mean alpha
      x1M_set = as.numeric(min(subset$length_um))
      x2M_set = as.numeric(max(subset$length_um))
      alpha = as.numeric(subset$alpha[1])
      if(alpha == 1){
        CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
        CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
        CF <- (as.numeric(CFL) + as.numeric(CFU))/2
      }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                        x2M = x2M_set, #upper measured length
                        x1D = x1D_set, #default lower size range
                        x2D = x2D_set,  #default upper size range
                        a = alpha #alpha
      )
      }
      
      subset$correction_factor <- as.numeric(CF)
      subset$corrected_concentration <- as.numeric(subset$correction_factor) * as.numeric(subset$concentration)
      
      #min alpha
      alpha = as.numeric(subset$alpha_lower[1])
      if(alpha == 1){
        CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
        CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
        CF <- (as.numeric(CFL) + as.numeric(CFU))/2
      }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                        x2M = x2M_set, #upper measured length
                        x1D = x1D_set, #default lower size range
                        x2D = x2D_set,  #default upper size range
                        a = alpha #alpha
      )
      }
      subset$correction_factor_lower <- as.numeric(CF)
      subset$corrected_concentration_lower <- as.numeric(subset$correction_factor_lower) * as.numeric(subset$concentration)
      
      #max alpha
      alpha = as.numeric(subset$alpha_upper[1])
      if(alpha == 1){
        CFL <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 0.999)
        CFU <- CFfnx(x1M = x1M_set, x2M = x2M_set, x1D = x1D_set, x2D = x2D_set, a = 1.001)
        CF <- (as.numeric(CFL) + as.numeric(CFU))/2
      }else{CF <- CFfnx(x1M = x1M_set,#lower measured length
                        x2M = x2M_set, #upper measured length
                        x1D = x1D_set, #default lower size range
                        x2D = x2D_set,  #default upper size range
                        a = alpha #alpha
      )
      }
      subset$correction_factor_upper <- as.numeric(CF)
      subset$corrected_concentration_upper <- as.numeric(subset$correction_factor_upper) * as.numeric(subset$concentration)
      
      
      
      dataframe_ <- rbind(dataframe_, subset)
      
      
      
    }
    dataframe <- dataframe_
  }
  
  return(dataframe)
}

#create function to extract midpoints from binning outputs
get_midpoints <- function(x){
  lower <- as.numeric(gsub(",.*","",gsub("\\(|\\[|\\)|\\]","", x)))
  upper <- as.numeric(gsub(".*,","",gsub("\\(|\\[|\\)|\\]","", x)))
  return(round(lower+(upper-lower)/2, 2))
}

#Make df for alpha values
study_media <- c("marinesurface","freshwatersurface","marinesediment","freshwatersediment","effluent", "biota", "drinkingwater")
length <- c(2.07, 2.64, 2.57, 3.25, 2.54, 2.59, 1.64)
length_sd <- c(0.03, 0.01, 0.20, 0.19, 0.01, 0.04, 0.55)
width <- c(1.96, 2.70, 2.51, 2.87, 2.55, 2.52, NA)
width_sd <- c(0.026, 0.008, 0.091, 0.167, 0.015, 0.07, NA)
mass <- c(1.32, 1.65, 1.50, 1.56, 1.40, 1.41, NA)
mass_sd <- c(0.009, 0.071, 0.026, 0.077, 0.002, 0.060, NA)
volume <- c(1.48, 1.68, 1.50, 1.53, 1.45, 1.40, NA)
volume_sd <- c(0.063, 0.081, 0.023, 0.013, 0.004, 0.004, NA)
surface_area <- c(1.50, 2.00, 1.75, 1.89, 1.73, 1.69, NA)
surface_area_sd <- c(0.009, 0.065, 0.050, 0.055, 0.005, 0.023, NA)
specific_surface_area <- c(1.98, 2.71, 2.54, 2.82, 2.58, 2.46, NA)
specific_surface_area_sd <- c(0.297, 0.009, 0.082, 0.096, 0.013, 0.092, NA)

alpha_vals <- data.frame(study_media = study_media,
                         length = length,
                         length_sd = length_sd,
                         width = width,
                         width_sd = width_sd,
                         mass=mass,
                         mass_sd = mass_sd,
                         volume = volume,
                         volume_sd = volume_sd,
                         surface_area = surface_area,
                         surface_area_sd = surface_area_sd,
                         specific_surface_area = specific_surface_area,
                         specific_surface_area_sd = specific_surface_area_sd
)

#Find average coefficient of variability to create upper and lower alpha vals from given values
length_coef <- as.numeric(mean(length_sd/length)) #0.0289907
width_coef <- as.numeric(mean(width_sd/width)) #0.02405526
mass_coef <- as.numeric(mean(mass_sd/mass)) #0.02675376
volume_coef <- as.numeric(mean(volume_sd/volume)) #0.02003795
surface_area_coef <- as.numeric(mean(surface_area_sd/surface_area)) #0.0187786
specific_surface_area_coef <- as.numeric(mean(specific_surface_area_sd/specific_surface_area)) #0.0436807


use_cases <- read.csv("data/Item_Use_Case.csv")
prime_unclassifiable <- read.csv("data/PrimeUnclassifiable.csv")
PrimeUnclassifiable <- read.csv("data/PrimeUnclassifiable.csv")

primeItems <- read.csv("data/PrimeItems.csv")
primeMaterials <- read.csv("data/PrimeMaterials.csv")

#Build cleaning functions
cleantext <- function(x) {
  x <- tolower(gsub("[[:space:]]", "", x))
  ifelse(x == "", NA, x)
}

removeslash <- function(x){
  gsub("/", " OR ", x)
}

cleanmaterials <- function(x) {
  x <- x[!names(x) %in% c("X1", "X2", "X3", "X4", "X5")]
  if(is.list(x)) lapply(x, cleanmaterials)
}

cleanitems <- function(x) {
  x <- x[!names(x) %in% c("X1", "X2", "X3", "X4", "X5", "X6")]
  if(is.list(x)) lapply(x, cleanitems)
}

#Build bootstrapping functions

BootMean <- function(data) {
  B <- 10000
  mean <- numeric(B)
  n = length(data)
  
  set.seed(34347)
  for (i in 1:B) {
    boot <- sample(1:n, size=n, replace = TRUE)
    mean[i] <- mean(data[boot], na.rm = T)
  }
  return(quantile(mean, c(0.025, 0.5, 0.975), na.rm = T))
}

confidence_interval_width <- function(data){
  proportion = 0.95
  sample_size = length(data)
  population_size = 10000
  1.96*abs(sqrt((1/sample_size)*proportion * (1-proportion) * (population_size-sample_size)/(population_size-1)))
}

AggregateTrees <- function(DF, Alias, Hierarchy){
  
  DF <- mutate_all(DF, cleantext) %>%
    mutate(Count = as.numeric(Count))
  
  colnames(DF) <- c("Alias", "Count")
  
  DF$Alias <- as.character(DF$Alias)
  
  Hierarchy <- mutate_all(Hierarchy, cleantext)
  
  colnames(Hierarchy) <- c("from", "Key")
  
  Alias <- mutate_all(Alias, cleantext) 
  
  DF <- DF %>% group_by(Alias) %>%
    summarise(across(Count, sum))
  
  
  DF_v2 <- DF %>%
    left_join(Alias) %>%
    dplyr::select(Key, Count) %>%
    mutate(Key = ifelse(is.na(Key), "other", Key)) %>%
    right_join(Hierarchy) %>%
    select(from, Key, Count) %>%
    add_row(from = "trash", Key = "missing", Count = sum(DF$Count, na.rm = T) - sum(.$Count, na.rm = T))%>%
    mutate(Count = Count/sum(Count, na.rm = T))
  
  
  DF_network <- FromDataFrameNetwork(DF_v2, check = "check")
  
  DF_network$Do(function(x) x$totalsum <- ifelse(is.null(x$Count), 0, x$Count) + sum(Get(x$children, "totalsum")), traversal = "post-order")
  
  Treedf <- ToDataFrameNetwork(DF_network, "totalsum") 
  
  Treedf %>%
    add_row(from = "trash", to = "", totalsum = sum(Treedf %>%
                                                      filter(from == "trash") %>% 
                                                      pull(totalsum)))
  
  
}

grouped_uncertainty <- function(DF_group, Group_Alias, Group_Hierarchy, type){
  
  df_join = data.frame(from = character(), 
                       to = character(), 
                       count = numeric())
  DF_group$Count = as.numeric(DF_group$Count)
  
  groups <- DF_group
  
  for(row in 1:nrow(groups)){
    df_subset <- DF_group %>%
      select(Class, Count)
    
    df_join <- AggregateTrees(DF = df_subset, Alias = Group_Alias, Hierarchy = Group_Hierarchy) %>%
      mutate(from = ifelse(from == "trash", type, from)) %>%
      bind_rows(df_join)
  }
  
  df_join_boot <- df_join %>%
    group_by(from, to) %>%
    summarise(mean_prop = mean(totalsum, na.rm = T), 
              min_prop = mean(totalsum, na.rm = T) - confidence_interval_width(totalsum), 
              max_prop = mean(totalsum, na.rm = T) + confidence_interval_width(totalsum))
  
}

sunburstplot <-function(df_join_boot){
  
  values <- paste(df_join_boot$to, 
                  "<br>", 
                  round(df_join_boot$mean_prop, 2) * 100, 
                  " (", 
                  round(df_join_boot$min_prop, 2) * 100, 
                  "-", 
                  round(df_join_boot$max_prop, 2) * 100, 
                  ")%", 
                  sep = "")
  
  values[df_join_boot$mean_prop < 0.07] <- NA
  
  plot_ly() %>%
    add_trace(
      labels = df_join_boot$to,
      parents = df_join_boot$from,
      type = "sunburst",
      maxdepth = 6,
      domain = list(column = 1), 
      branchvalues = 'total',
      texttemplate = values,
      values = df_join_boot$mean_prop) %>%
    layout(colorway = c("#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072", "#80B1D3", "#FDB462", "#B3DE69", "#FCCDE5", "#D9D9D9", "#BC80BD"))
    
}


###create function to derive correction factor (CF) from Koelmans et al (equation 2)
CFfnx = function(a, #default alpha from Koelmans et al (2020)
                 x2D, #set default values to convert ranges to (1-5,000 um) #5mm is upper default 
                 x1D, #1 um is lower default size
                 x2M, x1M){
  CF = (x2D^(1-a)-x1D^(1-a))/(x2M^(1-a)-x1M^(1-a)) 
  return(CF)}

#create function to locate mean value across bounded inverse power function (mean value theorum for integrals)
mean_val = function(a, #alpha, either solved for or default 1.6
                 x_b, #upper lim provided
                 x_a #lower lim provided
                 ){
  if(a == 1){
    c = ((1/(x_b - x_a))*(log(x_b) - log(x_a)))^(-1)
  }else{
    f_b = ((x_b)^(-a+1))/(-a+1)
    f_a = ((x_a)^(-a+1))/(-a+1)
    f_c = (1/(x_b - x_a))*(f_b - f_a)
    c = (abs(f_c))^(-(1/a)) 
  }
  return(c)}

#Files for tool
alias <- read.csv("data/PrimeMaterials.csv")
hierarchy <- read.csv("data/MaterialsHierarchyLower.csv")
aliasi <- read.csv("data/PrimeItems.csv")
hierarchyi <- read.csv("data/ITEMSHierarchyLower.csv")
microcolor <- read.csv("data/Microplastics_Color.csv")
trash_mass <- read.csv("data/trash_mass.csv")
aliasclean <- mutate_all(alias, cleantext)
aliascleani <- mutate_all(aliasi, cleantext)
hierarchyclean <- mutate_all(hierarchy, cleantext)
hierarchycleani <- mutate_all(hierarchyi, cleantext)
microcolorclean <- mutate_all(microcolor, cleantext)
trash_mass_clean <- mutate_all(trash_mass, cleantext)

#Creating materials hierarchy
Materials <- hierarchy
Materials[is.na(Materials)] <- ""
Materials <- mutate_all(Materials, removeslash)
Materials$pathString <- paste("Trash", Materials$X1, Materials$X2, Materials$X3, Materials$X4, Materials$X5, sep = "/")
Materials <- as.Node(Materials[,])
Materials <- as.list(Materials)
Materials <- Materials[-1]
Materials <- cleanmaterials(Materials)

Materials_hierarchy <- hierarchy
Materials_hierarchy[is.na(Materials_hierarchy)] <- ""
Materials_hierarchy <- mutate_all(Materials_hierarchy, removeslash) %>%
  mutate(key = "trash") %>%
  relocate(key) %>%
  unite(pathString, sep = "/")
Materials_hierarchy <- as.Node(Materials_hierarchy, pathDelimiter = "/")
Materials_hierarchy <- as.list(Materials_hierarchy)
Materials_hierarchy <- Materials_hierarchy[-1]

#Creating items hierarchy
Items <- hierarchyi
Items[is.na(Items)] <- ""
Items <- mutate_all(Items, removeslash)
Items$pathString <- paste("Trash", Items$X1, Items$X2, Items$X3, Items$X4, Items$X5, Items$X6, sep = "/")
Items <- as.Node(Items)
Items <- as.list(Items)
Items <- Items[-1]
Items <- cleanitems(Items)

Items_hierarchy <- hierarchyi
Items_hierarchy[is.na(Items_hierarchy)] <- ""
Items_hierarchy <- mutate_all(Items_hierarchy, removeslash) %>%
  mutate(key = "trash") %>%
  relocate(key) %>%
  unite(pathString, sep = "/")
Items_hierarchy <- as.Node(Items_hierarchy, pathDelimiter = "/")
Items_hierarchy <- as.list(Items_hierarchy)
Items_hierarchy <- Items_hierarchy[-1]

#Files for display
Materials_Alias <- read.csv("data/PrimeMaterials.csv")
Materials_Hierarchy <- read.csv("data/MaterialsHierarchyLower.csv")
Items_Alias <- read.csv("data/PrimeItems.csv")
Items_Hierarchy <- read.csv("data/ITEMSHierarchyLower.csv")
Material_Item_Relation <- read.csv("data/MaterialItemRelationship.csv")
Brand_Manufacturer_Relation <- read.csv("data/BrandManufacturer.csv")
Brand_Item_Relation <- read.csv("data/BrandItem.csv")
NOAA <- read.csv("data/NOAA.csv")
Micro_Color_Display <-read.csv("data/Microplastics_Color.csv")
MicroOnly <- read.csv("data/PremadeSurveys/Most_Specific_Microplastics.csv")
AllMore <- read.csv("data/PremadeSurveys/Most_Specific_All.csv")
AllLess <- read.csv("data/PremadeSurveys/Least_Specific_All.csv")
polymer_db <- read.csv("data/median_polymer_density.csv")
#polymer_db <- read.csv("data/10_90_Percentile_Polymer_Density.csv")

#make item and material pathstrings for merging tool
pathstrings_items <- data.frame(matrix(ncol=2, dimnames = list("", c("items", "pathString"))))
for(y in 1:ncol(hierarchycleani)){
  for(x in 1:nrow(hierarchycleani)){
    pathstrings_items[nrow(pathstrings_items) + 1, 1] <- hierarchycleani[x,y]
    if(!is.na(hierarchycleani[x,y]) && y == 1){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], sep = "/")}
    if(!is.na(hierarchycleani[x,y]) && y == 2){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], hierarchycleani[x,2], sep = "/")}
    if(!is.na(hierarchycleani[x,y]) && y == 3){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], hierarchycleani[x,2], hierarchycleani[x,3], sep = "/")}
    if(!is.na(hierarchycleani[x,y]) && y == 4){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], hierarchycleani[x,2], hierarchycleani[x,3], hierarchycleani[x,4], sep = "/")}
    if(!is.na(hierarchycleani[x,y]) && y == 5){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], hierarchycleani[x,2], hierarchycleani[x,3], hierarchycleani[x,4], hierarchycleani[x,5], sep = "/")}
    if(!is.na(hierarchycleani[x,y]) && y == 6){pathstrings_items[nrow(pathstrings_items), 2] <- paste("Trash", hierarchycleani[x,1], hierarchycleani[x,2], hierarchycleani[x,3], hierarchycleani[x,4], hierarchycleani[x,5], hierarchycleani[x,6], sep = "/")}
  }
  pathstrings_items <- pathstrings_items %>% distinct() %>% drop_na()
}

pathstrings_materials <- data.frame(matrix(ncol=2, dimnames = list("", c("material", "pathString"))))
for(y in 1:ncol(hierarchyclean)){
  for(x in 1:nrow(hierarchyclean)){
    pathstrings_materials[nrow(pathstrings_materials) + 1, 1] <- hierarchyclean[x,y]
    if(!is.na(hierarchyclean[x,y]) && y == 1){pathstrings_materials[nrow(pathstrings_materials), 2] <- paste("Trash", hierarchyclean[x,1], sep = "/")}
    if(!is.na(hierarchyclean[x,y]) && y == 2){pathstrings_materials[nrow(pathstrings_materials), 2] <- paste("Trash", hierarchyclean[x,1], hierarchyclean[x,2], sep = "/")}
    if(!is.na(hierarchyclean[x,y]) && y == 3){pathstrings_materials[nrow(pathstrings_materials), 2] <- paste("Trash", hierarchyclean[x,1], hierarchyclean[x,2], hierarchyclean[x,3], sep = "/")}
    if(!is.na(hierarchyclean[x,y]) && y == 4){pathstrings_materials[nrow(pathstrings_materials), 2] <- paste("Trash", hierarchyclean[x,1], hierarchyclean[x,2], hierarchyclean[x,3], hierarchyclean[x,4], sep = "/")}
    if(!is.na(hierarchyclean[x,y]) && y == 5){pathstrings_materials[nrow(pathstrings_materials), 2] <- paste("Trash", hierarchyclean[x,1], hierarchyclean[x,2], hierarchyclean[x,3], hierarchyclean[x,4], hierarchyclean[x,5], sep = "/")}
  }
  pathstrings_materials <- pathstrings_materials %>% distinct() %>% drop_na()
}

#Files for bootstrapping routine and sunburst plots
ItemsHierarchy_sunburst <- data.frame(matrix(ncol=2, dimnames = list("", c("from", "to"))))
for(y in 1:ncol(hierarchycleani)){
  for(x in 1:nrow(hierarchycleani)){
    ItemsHierarchy_sunburst[nrow(ItemsHierarchy_sunburst) + 1, 2] <- hierarchycleani[x,y]
    if(!is.na(hierarchycleani[x,y]) && y == 1){ItemsHierarchy_sunburst[nrow(ItemsHierarchy_sunburst), 1] <- paste("Trash")}
    if(!is.na(hierarchycleani[x,y]) && y != 1){ItemsHierarchy_sunburst[nrow(ItemsHierarchy_sunburst), 1] <- paste(hierarchycleani[x, y -1])}
  }
  ItemsHierarchy_sunburst <- ItemsHierarchy_sunburst %>% distinct() %>% drop_na()
}

MaterialsHierarchy_sunburst <- data.frame(matrix(ncol=2, dimnames = list("", c("from", "to"))))
for(y in 1:ncol(hierarchyclean)){
  for(x in 1:nrow(hierarchyclean)){
    MaterialsHierarchy_sunburst[nrow(MaterialsHierarchy_sunburst) + 1, 2] <- hierarchyclean[x,y]
    if(!is.na(hierarchyclean[x,y]) && y == 1){MaterialsHierarchy_sunburst[nrow(MaterialsHierarchy_sunburst), 1] <- paste("Trash")}
    if(!is.na(hierarchyclean[x,y]) && y != 1){MaterialsHierarchy_sunburst[nrow(MaterialsHierarchy_sunburst), 1] <- paste(hierarchyclean[x, y -1])}
  }
  MaterialsHierarchy_sunburst <- MaterialsHierarchy_sunburst %>% distinct() %>% drop_na()
}

ItemsAlias_sunburst <- read.csv("data/PrimeItems.csv")%>%
  rename(Key = Item) %>%
  select(-readable)
MaterialsAlias_sunburst <- read.csv("data/PrimeMaterials.csv") %>%
  rename(Key = Material) %>%
  select(-readable)

#Test data
particle_testData <- read.csv("tests/rescaling_particle.csv", check.names = FALSE)
sample_testData <- read.csv("tests/count_mass_concentration.csv", check.names = FALSE)

#data tables for count to mass conversion
#Make polymer-density dataframe
#Output correct survey sheet
polymer_db_ <- data.frame(polymer_db)
polymer_db_$material <- cleantext(polymer_db_$material)
polymer_db_$density <- as.numeric(polymer_db_$density)
polymer_db_$lower_conf <- as.numeric(polymer_db_$lower_conf)
polymer_db_$upper_conf <- as.numeric(polymer_db_$upper_conf)
density_mg_um_3 <- polymer_db_$density * 1e-9
density_min <- polymer_db_$lower_conf * 1e-9
density_max <- polymer_db_$upper_conf * 1e-9
polymer_db_ <- polymer_db_ %>%
  mutate(density_mg_um_3 = density_mg_um_3,
         density_max = density_max,
         density_min = density_min)
polymer_db_ <- polymer_db_ %>% 
  left_join(Materials_Alias, by = c("material" = "Alias")) %>%
  select(-c("Material"))
polymer_density <- polymer_db_ %>%
  select(material, density_mg_um_3, density_max, density_min, readable)

#Average densities for parent polymer classes
polymer_parents <- polymer_db_ %>% left_join(MaterialsHierarchy_sunburst, by = c("material" = "to")) %>%
  rename("parent" = "from") %>%
  drop_na(parent)
polymer_parents <- polymer_parents %>%
  group_by(parent) %>%
  summarize(density_mg_um_3 = mean(density_mg_um_3, na.rm=TRUE), density_max = mean(density_max, na.rm=TRUE), density_min = mean(density_min, na.rm=TRUE)) %>%
  rename(Material = parent)
polymer_parents <- polymer_parents %>%
  left_join(primeMaterials, by = "Material") %>%
  select(-Alias) %>%
  rename(material = Material)
  
#Join density dbs
polymer_density <- rbind(polymer_density, polymer_parents)
polymer_density <- polymer_density %>% distinct(material, .keep_all = TRUE)

#Add +/- 5% error for measured particle dimensions
morphology <- c("fragment","sphere","fiber","film","foam")
meas_min <- 0.95
meas_max <- 1.05
L_min <- c(0.95, 0.95, 0.95, 0.95, 0.95)
L_max <- c(1.05, 1.05, 1.05, 1.05, 1.05)

#Min and max values given in Kooi Koelmans
W_min <- c(0.1,0.60,0.001,0.1,0.1)
W_max <- c(1,1,0.5,1,1)
H_min <- c(0.01,0.36,0.001,0.001,0.01)
H_max <- c(1,1,0.5,0.1,1)

morphology_base <- data.frame(morphology=morphology,
                               L_min=L_min,
                               L_max=L_max,
                               W_min=W_min,
                               W_max=W_max,
                               H_min=H_min,
                               H_max=H_max
)

#convert morphologies in TT to morphologies with defined dimensions
morphology_unique <- c("fiber", "nurdle", "foam", "sphere", "line", "bead", "pellet", "sheet", "film", "fragment", "rubberyfragment", "fiberbundle")
morphology <- c("fiber", "sphere", "foam", "sphere",  "fiber", "sphere", "sphere", "film", "film", "fragment", "fragment", "film")
morph_conversion <- data.frame(morphology = morphology,
                               morphology_unique = morphology_unique)
morphology_shape <- morph_conversion %>% left_join(morphology_base, by = "morphology") %>% select(-morphology) %>% rename(morphology = morphology_unique)

#polymer abundance from Burns et al 2018/Bond et al 2018 (ref in Kooi Koelmans 2019 for continuous polymer distribution)
polymer <- c("PE", "PP", "PET", "PA", "PS", "PVC", "PVA")
#Abundance provided does not equal 100, so taken as a proportion
abundance <- c(0.25/0.845, 0.145/0.845, 0.165/0.845, 0.12/0.845, 0.085/0.845, 0.02/0.845, 0.06/0.845)
polymer_abundance <- data.frame(polymer = polymer,
                                abundance = abundance)
polymer_abundance <- polymer_abundance %>% left_join(polymer_density, by = c("polymer" = "readable")) %>% select(-material)
polymer_abundance <- polymer_abundance %>% mutate(density_mg_um_3 = density_mg_um_3 * abundance,
                                                  density_max = density_max * abundance,
                                                  density_min = density_min * abundance)

polymer_abundance <- polymer_abundance %>% add_row(polymer = "average", 
                                                   density_mg_um_3 = sum(polymer_abundance$density_mg_um_3), 
                                                   density_max = sum(polymer_abundance$density_max), 
                                                   density_min = sum(polymer_abundance$density_min))
# c(1.148166e-09, 1.523589e-09, 1.037747e-09)



#morphology abundance from Burns et al 2018/Bond et al 2018 (ref in Kooi Koelmans 2019 for continuous Corey Shape Factor distribution)
morphology <- c("fiber", "fragment", "sphere", "film", "foam")
#Abundance provided does not equal 100, so taken as a proportion
abundance <- c(0.485/0.95, 0.31/0.95, 0.065/0.95, 0.055/0.95, 0.035/0.95)
morphology_abundance <- data.frame(morphology = morphology,
                                   abundance = abundance)
morphology_abundance <- morphology_abundance %>% left_join(morphology_shape, by = "morphology")
morphology_abundance <- morphology_abundance %>% mutate(L_min = L_min * abundance, 
                                                        L_max = L_max * abundance,
                                                        W_min = W_min * abundance,
                                                        W_max = W_max * abundance,
                                                        H_min = H_min * abundance,
                                                        H_max = H_max * abundance)
morphology_abundance <- morphology_abundance %>% add_row(morphology = "average", 
                                                         L_min = sum(morphology_abundance$L_min), 
                                                         L_max = sum(morphology_abundance$L_max), 
                                                         W_min = sum(morphology_abundance$W_min),
                                                         W_max = sum(morphology_abundance$W_max),
                                                         H_min = sum(morphology_abundance$H_min),
                                                         H_max = sum(morphology_abundance$H_max))
#c(0.950, 1.05000000, 0.0836684211, 0.74473684, 2.883158e-02, 0.692631579)


#UI Elements ----
titles <- c(
  "Materials Alias Table",
  "Materials Hierarchy Table",
  "Items Alias Table",
  "Items Hierarchy Table",
  "Material-Item Relational Table",
  "Material-Density Relational Table",
  "Manufacturer Brand Relational Table",
  "Item-Brand Relational Table",
  "Misaligned Categories Table",
  "Color Alias"
)

captions <- c(
  "This table describes the aliases that can be used to describe material types and links them to a key term. Each row represents a unique material and each column is an alias for that material.",
  "This table describes how the unique material types relate to one another in a hierarchical structure (ex: foam and rubber are a subset of plastic).",
  "This table describes the aliases that can be used to describe item types and links them to a key term. Each row represents a unique item and each column is an alias for that item.",
  "This table describes how the unique items relate to one another in a hierarchical structure (ex: forks, knives, and spoons all fall under utensils).",
  "This table relates the items, materials, and survey sheets used to make the other relational tables.",
  "This table relates materials to all known material densities. Sources also displayed.",
  "This table relates brand types to their respective manufacturer.",
  "This table relates brands to items.",
  "This table displays all categories that did not fit our item-material framework.",
  "This table describes the aliases that can be used to describe particle colors. Each row represents a unique color and each column is an alias for that color."
)

view_code <- c(
  "table1",
  "table2",
  "table3",
  "table4",
  "table5",
  "table9",
  "table6",
  "table7",
  "table8",
  "table12"
)

files <- list(Materials_Alias, 
              Materials_Hierarchy,
              Items_Alias,
              Items_Hierarchy,
              Material_Item_Relation, 
              polymer_db, 
              Brand_Manufacturer_Relation,
              Brand_Item_Relation,
              PrimeUnclassifiable,
              color
)