#Code LCA Group Project D 
library(ggplot2)
library(dplyr)
library(lubridate)
library(rPref)
library(MASS)
library(scales)  
library(mco)


#Global Parameters

start.Date <- ymd("2026-01-01")
lifetime <- 80

event_duration <- c(
  SR.n = 12, VI.n = 1, DR.n = 28,
  SR.s = 5, VI.s = 1, DR.s = 10,
  SR.p = 7, VI.p = 1, DR.p = 15,
  SR.t = 10, VI.t = 1, DR.t = 21
)


#Event generations

generate_events <- function(events_list, lifetime){
  df_sys <- data.frame(Freq = numeric(), Event = character(), System = character())
  for(sys in names(events_list)){
    evt <- events_list[[sys]]
    for(e in names(evt)){
      occ <- seq(0, lifetime, by = evt[e])
      occ <- occ[occ != 0]  # ignor DC (0)
      df <- data.frame(Freq = occ, Event = e, System = sys)
      df_sys <- rbind(df_sys, df)
    }
  }
  dc_end <- do.call(rbind, lapply(names(events_list), function(sys){
    data.frame(Freq = c(0, lifetime), Event = c("DC","END"), System = sys)
  }))
  df_sys <- rbind(df_sys, dc_end)
  
  df_sys <- df_sys %>%
    mutate(Duration = ifelse(Event %in% names(event_duration),
                             event_duration[Event], 0))
  
  return(df_sys)
}



compute_total_days <- function(events_list, lifetime){
  global.events <- generate_events(events_list, lifetime)
  
  total_days <- global.events %>%
    filter(!Event %in% c("DC", "END")) %>%
    filter(Freq < lifetime) %>%  # ignore interventions exactly at the end
    group_by(Freq) %>%
    summarise(Day = max(Duration)) %>%
    summarise(Total_days = sum(Day)) %>%
    pull(Total_days)
  
  return(total_days)
}


# Create the plot

plot_timeline <- function(global.events, title){
  system_levels <- data.frame(System = unique(global.events$System),
                              y = seq(length(unique(global.events$System)),1,-1))
  global.events <- merge(global.events, system_levels, by="System")
  
  # Color by suffixe
  all_events <- unique(global.events$Event)
  event_colors <- c("DC"="black", "END"="black")
  suffix_colors <- list(
    ".n" = c(0.50, 0.55, 0.62),
    ".p" = c(0.35, 0.22, 0.45),
    ".s" = c(0.99, 0.15, 0.11),
    ".t" = c(0.75, 0.80, 0.85)
  )
  
  for(suf in names(suffix_colors)){
    ev_suf <- all_events[grepl(suf, all_events)]
    n_ev <- length(ev_suf)
    hues <- suffix_colors[[suf]]
    if(n_ev > length(hues)) hues <- seq(hues[1], hues[length(hues)], length.out = n_ev)
    pal <- sapply(hues[1:n_ev], function(h) hsv(h, 0.8, 0.9))
    event_colors[ev_suf] <- pal
  }
  
  # Priority of END
  segments <- global.events %>%
    group_by(Freq, System) %>%
    mutate(has_end = "END" %in% Event) %>%
    filter(!(has_end & Event != "END")) %>%
    mutate(
      n_events = n(),
      row_id = row_number(),
      ymin = y - 0.3 + 0.6 * (row_id-1)/n_events,
      ymax = y - 0.3 + 0.6 * row_id/n_events
    ) %>%
    ungroup()
  
  all_years <- seq(0, lifetime, by = 1)
  every5_years <- seq(0, lifetime, by = 5)
  
  legend_order <- c(
    "DC", "END",
    names(events_base$Nuclear),
    names(events_base$Sidewalk),
    names(events_base$Parking),
    names(events_base$Tunnel)
  )
  legend_order <- legend_order[legend_order %in% names(event_colors)]
  
  
  ggplot(segments, aes(x = Freq)) +
    geom_vline(xintercept = all_years, color = "grey85", size = 0.3) +
    geom_vline(xintercept = every5_years, color = "grey50", size = 0.7) +
    geom_segment(aes(y = ymin, yend = ymax, xend = Freq, color = Event), size = 1.5) +
    scale_x_continuous(
      breaks = every5_years,
      labels = function(x) year(start.Date + years(x))
    ) +
    scale_y_continuous(breaks = system_levels$y, labels = system_levels$System) +
    scale_color_manual(values = event_colors, breaks = legend_order) +
    labs(x = "Years", y = "", title = title) +
    theme_minimal() +
    theme(
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank(),
      axis.title.y = element_blank(),
      plot.title = element_text(face = "bold", hjust = 0.5),
      legend.title = element_blank()
    )
}


# Example

events_base <- list(
  Nuclear = c(VI.n = 4, SR.n = 11, DR.n = 26),
  Sidewalk = c(VI.s = 3, SR.s = 10, DR.s = 20),
  Parking = c(VI.p = 6, SR.p = 12, DR.p = 20),
  Tunnel = c(VI.t = 5, SR.t = 10, DR.t = 25)
)

global.events_base <- generate_events(events_base, lifetime)
compute_total_days(events_base, lifetime) 
plot_timeline(global.events_base, "Combined Systems Planning Maintenance Strategies")


# Modifications


#test 1: 
events_test1 <- list(
Nuclear = c(VI.n = 4, SR.n = 11, DR.n = 26),
Sidewalk = c(VI.s = 4, SR.s = 10, DR.s = 20),
Parking = c(VI.p = 6, SR.p = 12, DR.p = 20),
Tunnel = c(VI.t = 5, SR.t = 10, DR.t = 25)
)


global.events_test1 <- generate_events(events_test1, lifetime)
compute_total_days(events_test1, lifetime)  
plot_timeline(global.events_test1, "Test Scenario n°1 – One-Year Longer Interval of the Most Recurrent Intervention (SR.s)")

events_test2 <- list(
Nuclear = c(VI.n = 4, SR.n = 11, DR.n = 26),
Sidewalk = c(VI.s = 2, SR.s = 10, DR.s = 20),
Parking = c(VI.p = 6, SR.p = 12, DR.p = 20),
Tunnel = c(VI.t = 5, SR.t = 10, DR.t = 25)
)

global.events_test2 <- generate_events(events_test2, lifetime)
compute_total_days(events_test2, lifetime)  
plot_timeline(global.events_test2, "Test Scenario n°2 – One-Year Shorter Interval of the Most Recurrent Intervention (SR.s)")

events_test3 <- list(   
Nuclear = c(VI.n = 3, SR.n = 11, DR.n = 26),
Sidewalk = c(VI.s = 3, SR.s = 10, DR.s = 20),
Parking = c(VI.p = 3, SR.p = 12, DR.p = 20),
Tunnel = c(VI.t = 3, SR.t = 10, DR.t = 25)
)


global.events_test3 <- generate_events(events_test3, lifetime)
compute_total_days(events_test3, lifetime)  
plot_timeline(global.events_test3, "Test Scenario n°3 – Harmonization of All Surface Repairs (SR) to a 3-Year Interval")


events_test4 <- list(
Nuclear = c(VI.n = 3, SR.n = 10, DR.n = 20),
Sidewalk = c(VI.s = 3, SR.s = 10, DR.s = 20),
Parking = c(VI.p = 3, SR.p = 10, DR.p = 20),
Tunnel = c(VI.t = 3, SR.t = 10, DR.t = 20)
)


global.events_test4 <- generate_events(events_test4, lifetime)
compute_total_days(events_test4, lifetime)  # <- total recalculé automatiquement
plot_timeline(global.events_test4, "Test Scenario n°4 – Minimum Interval for All Maintenance Interventions")


design_explore_multi <- function(events_list_grid, lifetime, n_samples = 6000) {
  
  param_names <- names(events_list_grid)
  results <- data.frame()
  
  for(k in 1:n_samples) {
    

    sampled_vals <- lapply(events_list_grid, function(v) sample(v, 1))
    
    events_scenario <- list(
      Nuclear = c(
        VI.n = sampled_vals$Nuclear.VI.n,
        SR.n = sampled_vals$Nuclear.SR.n,
        DR.n = sampled_vals$Nuclear.DR.n
      ),
      Sidewalk = c(
        VI.s = sampled_vals$Sidewalk.VI.s,
        SR.s = sampled_vals$Sidewalk.SR.s,
        DR.s = sampled_vals$Sidewalk.DR.s
      ),
      Parking = c(
        VI.p = sampled_vals$Parking.VI.p,
        SR.p = sampled_vals$Parking.SR.p,
        DR.p = sampled_vals$Parking.DR.p
      ),
      Tunnel = c(
        VI.t = sampled_vals$Tunnel.VI.t,
        SR.t = sampled_vals$Tunnel.SR.t,
        DR.t = sampled_vals$Tunnel.DR.t
      )
    )
    
    total_days <- compute_total_days(events_scenario, lifetime)
    global.events <- generate_events(events_scenario, lifetime)
    min_dist <- min(diff(sort(unique(global.events$Freq))))
    
    results <- rbind(results,
                     cbind(as.data.frame(sampled_vals),
                           Total_days = total_days,
                           Min_interval = min_dist))
  }
  
  return(results)
}



build_events_grid <- function(events_base, variation_pct = 0.2, step = 1) {
  
  grid_list <- list()
  
  for(system in names(events_base)) {
    for(event in names(events_base[[system]])) {
      
      base_val <- events_base[[system]][event]
      
      min_val <- max(1, floor(base_val * (1 - variation_pct)))
      max_val <- ceiling(base_val * (1 + variation_pct))
      
      seq_vals <- seq(min_val, max_val, by = step)
      
      grid_name <- paste(system, event, sep = ".")
      grid_list[[grid_name]] <- seq_vals
    }
  }
  
  return(grid_list)
}

pourcentage = 0.2
events_list_grid <- build_events_grid(events_base, variation_pct = pourcentage, step = 1)

str(events_list_grid)
response.space <- design_explore_multi(events_list_grid, lifetime, n_samples = 6000)
p <- low(Total_days) * high(Min_interval)
sky <- psel(response.space, p)
pareto2 <- psel(response.space, p, top = nrow(response.space))





getwd()
LCI.materials <- read.csv("LCImaterials.csv")

event.scope.map <- data.frame(
  event = c("SR.n","VI.n","DR.n",
            "SR.s","VI.s","DR.s",
            "SR.p","VI.p","DR.p",
            "SR.t","VI.t","DR.t"),
  scope = c("RC","RC","RC",   # Nuclear Containment
            "RC","RC","RC",   # Sidewalk
            "RC","RC","RC",   # Parking
            "RC","RC","RC"),  # Tunnel
  stringsAsFactors = FALSE
)

base_intervals <- c(
  VI.n = 4, SR.n = 11, DR.n = 26,
  VI.s = 3, SR.s = 10, DR.s = 20,
  VI.p = 6, SR.p = 12, DR.p = 20,
  VI.t = 5, SR.t = 10, DR.t = 25
)


LCA_RC_system <- function(volume, thickness, ref_thickness, materials, dist.event, event_prefix, base_intervals, lifetime){
  
  materials.split <- split(materials, materials$scope)
  RC.mat <- materials.split$RC
  

  SR_name <- paste0("SR.", event_prefix)
  DR_name <- paste0("DR.", event_prefix)
  VI_name <- paste0("VI.", event_prefix)
  

  days_per_year <- 365
  
  interval_SR <- base_intervals[SR_name] * (ref_thickness / thickness) * days_per_year
  interval_DR <- base_intervals[DR_name] * (ref_thickness / thickness) * days_per_year
  interval_VI <- base_intervals[VI_name] * (ref_thickness / thickness) * days_per_year
  
  n_SR <- ceiling(lifetime * days_per_year / interval_SR)
  n_DR <- ceiling(lifetime * days_per_year / interval_DR)
  n_VI <- ceiling(lifetime * days_per_year / interval_VI)
  
  

  interventions_ACV <- 0.4*n_SR + 0.8*n_DR + 0.1*n_VI
  

  RC.LCA <- dplyr::mutate(RC.mat,
                          Quantity = quantities * volume / 1000,
                          LC.mult = Quantity * interventions_ACV,
                          Energy.LC = LC.mult * energy,
                          CO2.LC    = LC.mult * CO2,
                          NOX.LC    = LC.mult * NOX,
                          SO2.LC    = LC.mult * SO2)
  
  return(list(
    Energy = sum(RC.LCA$Energy.LC),
    CO2    = sum(RC.LCA$CO2.LC),
    NOX    = sum(RC.LCA$NOX.LC),
    SO2    = sum(RC.LCA$SO2.LC),
    n_SR = n_SR,
    n_DR = n_DR,
    n_VI = n_VI
  ))
}




volume_parking  <- function(t)  3 * 17.6 * 48 * t
volume_sidewalk <- function(t)  1.8 * 250 * t
volume_tunnel   <- function(t)  700 * pi * ((5 + t)^2 - 5^2)
volume_nuclear  <- function(t)  pi * 37 * 59 * t

  

  



fitness_full <- function(x, materials, lifetime, base_intervals){
  

  events_scenario <- list(
    Nuclear  = c(VI.n = x[1], SR.n = x[2], DR.n = x[3]),
    Sidewalk = c(VI.s = x[4], SR.s = x[5], DR.s = x[6]),
    Parking  = c(VI.p = x[7], SR.p = x[8], DR.p = x[9]),
    Tunnel   = c(VI.t = x[10], SR.t = x[11], DR.t = x[12])
  )
  
  total_days <- compute_total_days(events_scenario, lifetime)
  global.events <- generate_events(events_scenario, lifetime)
  min_dist <- min(diff(sort(unique(global.events$Freq))))*365
  
  

  parking.ev  <- subset(global.events, System == "Parking")
  tunnel.ev   <- subset(global.events, System == "Tunnel")
  sidewalk.ev <- subset(global.events, System == "Sidewalk")
  nuclear.ev  <- subset(global.events, System == "Nuclear")
  

  t_parking  <- x[13]
  t_tunnel   <- x[14]
  t_sidewalk <- x[15]
  t_nuclear  <- x[16]
  

  v_parking  <- volume_parking(t_parking)
  v_tunnel   <- volume_tunnel(t_tunnel)
  v_sidewalk <- volume_sidewalk(t_sidewalk)
  v_nuclear  <- volume_nuclear(t_nuclear)
  


  LCA_parking  <- LCA_RC_system(v_parking,  t_parking,  0.25, materials, parking.ev,  "p", base_intervals, lifetime)
  LCA_tunnel   <- LCA_RC_system(v_tunnel,   t_tunnel,   0.30, materials, tunnel.ev,   "t", base_intervals, lifetime)
  LCA_sidewalk <- LCA_RC_system(v_sidewalk, t_sidewalk, 0.12, materials, sidewalk.ev, "s", base_intervals, lifetime)
  LCA_nuclear  <- LCA_RC_system(v_nuclear,  t_nuclear,  1.00, materials, nuclear.ev,  "n", base_intervals, lifetime)
  

  Energy <- LCA_parking$Energy + LCA_tunnel$Energy + LCA_sidewalk$Energy + LCA_nuclear$Energy
  CO2    <- LCA_parking$CO2    + LCA_tunnel$CO2    + LCA_sidewalk$CO2    + LCA_nuclear$CO2
  NOX    <- LCA_parking$NOX    + LCA_tunnel$NOX    + LCA_sidewalk$NOX    + LCA_nuclear$NOX
  SO2    <- LCA_parking$SO2    + LCA_tunnel$SO2    + LCA_sidewalk$SO2    + LCA_nuclear$SO2
  

  Energy.costs <- 0.128
  CO2.unitcost <- 26
  NOX.unitCost <- 42
  SO2.unitCosts <- 85
  
  cost <- (Energy * Energy.costs + CO2*CO2.unitcost + NOX*NOX.unitCost + SO2*SO2.unitCosts)/1e9
  
  
  return(c(total_days, -min_dist, Energy, CO2, NOX, SO2, cost))
}


t_ref <- c(
  0.25,  # parking  (x[13])
  0.30,  # tunnel   (x[14])
  0.12,  # sidewalk (x[15])
  1.00   # nuclear  (x[16])
)

variation <- 0.20

lower_thickness <- t_ref * (1 - variation)
upper_thickness <- t_ref * (1 + variation)

base_x <- c(
  base_intervals["VI.n"], base_intervals["SR.n"], base_intervals["DR.n"],
  base_intervals["VI.s"], base_intervals["SR.s"], base_intervals["DR.s"],
  base_intervals["VI.p"], base_intervals["SR.p"], base_intervals["DR.p"],
  base_intervals["VI.t"], base_intervals["SR.t"], base_intervals["DR.t"]
)

lower_intervals <- floor(base_x * (1 - variation))
upper_intervals <- ceiling(base_x * (1 + variation))
lower_intervals <- pmax(lower_intervals, 1)

lower <- c(lower_intervals, lower_thickness)
upper <- c(upper_intervals, upper_thickness)


set.seed(12345)

r_final <- nsga2(
  function(x) fitness_full(x, LCI.materials, lifetime, base_intervals),
  idim = 16,
  odim = 7,
  lower.bounds = lower,
  upper.bounds = upper,
  popsize = 400,        
  generations = 50      
)


adjusted_intervals <- function(par_row, base_intervals) {
  x <- par_row
  interval_SR_n <- x[2] * (1 / x[16])  
  interval_DR_n <- x[3] * (1 / x[16]) 
  interval_VI_n <- x[1] * (1 / x[16]) 
  
  interval_SR_s <- x[5] * (0.12 / x[15]) 
  interval_DR_s <- x[6] * (0.12 / x[15]) 
  interval_VI_s <- x[4] * (0.12 / x[15]) 
  
  interval_SR_p <- x[8] * (0.25 / x[13]) 
  interval_DR_p <- x[9] * (0.25 / x[13]) 
  interval_VI_p <- x[7] * (0.25 / x[13]) 
  
  interval_SR_t <- x[11] * (0.30 / x[14]) 
  interval_DR_t <- x[12] * (0.30 / x[14])
  interval_VI_t <- x[10] * (0.30 / x[14]) 
  
  return(c(
    VI.n = interval_VI_n, SR.n = interval_SR_n, DR.n = interval_DR_n,
    VI.s = interval_VI_s, SR.s = interval_SR_s, DR.s = interval_DR_s,
    VI.p = interval_VI_p, SR.p = interval_SR_p, DR.p = interval_DR_p,
    VI.t = interval_VI_t, SR.t = interval_SR_t, DR.t = interval_DR_t,
    t.parking = x[13], t.tunnel = x[14], t.sidewalk = x[15], t.nuclear = x[16]
  ))
}


df_par_real <- t(apply(r_final$par, 1, adjusted_intervals))
df_par_real <- as.data.frame(df_par_real)


df_val <- as.data.frame(r_final$value)
colnames(df_val) <- c(
  "Total_days",
  "Min_distance",
  "Energy","CO2","NOX","SO2","Cost"
)
df_val$Min_distance <- -df_val$Min_distance  
df_val$Total_days <- df_val$Total_days - 183   


all.results <- cbind(df_par_real, df_val)
all.results$pareto <- r_final$pareto.optimal

cols_plot <- which(sapply(all.results, is.numeric) & colnames(all.results) != "pareto")


pdf("parcoord_all_parameters_and_objectives_real_intervals.pdf",
    width = 34, height = 8)

parcoord(
  all.results[, cols_plot],
  var.label = TRUE,
  col = ifelse(all.results$pareto, "indianred", "skyblue2"),
  lty = ifelse(all.results$pareto, 1, 3),
  lwd = ifelse(all.results$pareto, 3, 1)
)

dev.off()

pareto3 <- as.data.frame(all.results[all.results$pareto, ])
pareto3 <- pareto3[, c("Total_days", "Cost")]

# Tri pour geom_line
pareto3 <- pareto3[order(pareto3$Total_days), ]

pdf("pareto_front_TotalDays_vs_Cost.pdf", width = 10, height = 7)
ggplot(all.results, aes(x = Total_days, y = Cost)) +
  geom_point(shape = 21, color = "skyblue2", fill = "skyblue2", alpha = 0.5) +
  geom_point(data = pareto3, aes(x = Total_days, y = Cost), size = 3, color = "red") +
  geom_line(data = pareto3, aes(x = Total_days, y = Cost), color = "blue", size = 1) +
  labs(title = "Pareto Front NSGA-II", x = "Total Days", y = "Cost [GCHF]") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
dev.off()


save(r_final, file = "NSGA2_results.RData")
saveRDS(all.results, file = "all_results.rds")
load("NSGA2_results.RData")
all.results <- readRDS("all_results.rds")



red_parcoord <- all.results[all.results$pareto, ]
print(red_parcoord[, cols_plot])
write.csv(
  red_parcoord[, cols_plot],
  file = "parcoord_red_solutions.csv",
  row.names = FALSE
)



print(pareto3)
red_parcoord_full <- all.results[all.results$pareto, ]
str(red_parcoord_full)
summary(red_parcoord_full)
write.csv(
  red_parcoord_full,
  file = "pareto_solutions_full_from_parcoord.csv",
  row.names = FALSE
)




cols_plot <- c("VI.n","SR.n","DR.n","VI.s","SR.s","DR.s","VI.p","SR.p","DR.p","VI.t","SR.t","DR.t",
               "Thickness_nuclear","Thickness_tunnel","Thickness_parking","Thickness_sidewalk",
               "Total_days","Min_interval","Energy","CO2","NOX","SO2","Cost")


pdf("parcoord_all_parameters_and_objectives_real_intervals.pdf", width = 34, height = 8)
parcoord(
  all.results[, cols_plot],
  var.label = TRUE,
  col = ifelse(all.results$pareto, "indianred", "skyblue2"),
  lty = ifelse(all.results$pareto, 1, 3),
  lwd = ifelse(all.results$pareto, 3, 1)
)
dev.off()


png("parcoord_all_parameters_and_objectives_real_intervals.png", width = 3400, height = 1000, res = 150)
parcoord(
  all.results[, cols_plot],
  var.label = TRUE,
  col = ifelse(all.results$pareto, "indianred", "skyblue2"),
  lty = ifelse(all.results$pareto, 1, 3),
  lwd = ifelse(all.results$pareto, 3, 1)
)
dev.off()



pareto3 <- as.data.frame(all.results[all.results$pareto, ])
pareto3 <- pareto3[, c("Total_days", "Cost")]


pareto3 <- pareto3[order(pareto3$Total_days), ]

# PDF
pdf("pareto_front_TotalDays_vs_Cost.pdf", width = 10, height = 7)
ggplot(all.results, aes(x = Total_days, y = Cost)) +
  geom_point(shape = 21, color = "skyblue2", fill = "skyblue2", alpha = 0.5) +
  geom_point(data = pareto3, aes(x = Total_days, y = Cost), size = 3, color = "red") +
  geom_line(data = pareto3, aes(x = Total_days, y = Cost), color = "blue", size = 1) +
  labs(title = "Pareto Front NSGA-II", x = "Total Days", y = "Cost [GCHF]") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
dev.off()

# PNG
png("pareto_front_TotalDays_vs_Cost.png", width = 1000, height = 700, res = 150)
ggplot(all.results, aes(x = Total_days, y = Cost)) +
  geom_point(shape = 21, color = "skyblue2", fill = "skyblue2", alpha = 0.5) +
  geom_point(data = pareto3, aes(x = Total_days, y = Cost), size = 3, color = "red") +
  geom_line(data = pareto3, aes(x = Total_days, y = Cost), color = "blue", size = 1) +
  labs(title = "Pareto Front NSGA-II", x = "Total Days", y = "Cost [GCHF]") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
dev.off()