library(tidyverse)
library(gt) #For nice color coded capacity limits table.
	if(!exists("SAVE_RES_LOC")){SAVE_RES_LOC <- "./Results/Primary_Simulation_Results/Main_Results/"}
	dir.create(SAVE_RES_LOC, recursive = TRUE, showWarnings = FALSE)


###Process the simulations and save the main percentile results by year
	RES <- read_csv("Results/Simulations/Kemmerer_2024_Simulation.csv")
	RES[,"Year"] <- RES[,"Year"]
	RES<-	RES %>% filter(!is.na(Year))
	RES <-	RES %>% filter(!is.na(Population))
	HIST <-	readRDS("Data/Cleaned_Data/Population_Data/RDS/Kemmerer_Diamondville_Population_Data.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>=1940)
	######Population
	####Fan New
	GET_DATA <- function(RES,COL_NUM){ 

	YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE)
	FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),seq(0.01,0.99,by=0.01))}))  %>% as_tibble %>% mutate(Year=YEARS)
	FAN_DATA <- 	rbind(FAN_DATA[1,],FAN_DATA)
	START_VALUE <- (HIST %>% filter(Year==2024))[,COL_NUM+1] %>% as.numeric
	FAN_DATA <- 	FAN_DATA  %>% pivot_longer(colnames(FAN_DATA %>% select(-Year)),names_to="Percentile") %>% filter(Year>2024) %>% unique
	NUM_YEARS <- length(unique(FAN_DATA$Year) )
	FAN_DATA$Group <- rep(c(1:49,0,rev(1:49)),NUM_YEARS)
	FAN_DATA <- FAN_DATA %>% group_by(Year,Group)  %>% summarize(MIN=min(value),MAX=max(value))
	TEMP <- FAN_DATA %>% filter(Year==2025) %>% mutate(Year=2024) %>% ungroup
	TEMP[,3:4] <- START_VALUE 
	FAN_DATA <- rbind(TEMP,FAN_DATA %>% ungroup) %>% as_tibble
	return(FAN_DATA)
	}
	RES %>% pull(Sim_UUID) %>% unique %>% length()
MAKE_GRAPH <- function(GRAPH_DATA,ALPHA=0.03,COLOR='cadetblue',LINE_WIDTH=0.75){
	PLOT <- ggplot(data=GRAPH_DATA)
	for(i in 1:49){
	C_DATA <- GRAPH_DATA%>% filter(Group==i)
	     PLOT <- 	PLOT +geom_ribbon(data=C_DATA,aes(x=Year,ymin=MIN,ymax=MAX),alpha=ALPHA,fill=COLOR)

}
	CI_90 <- rbind(GRAPH_DATA%>% filter(Group==20) %>% mutate(Interval='80%'),GRAPH_DATA%>% filter(Group==5) %>% mutate(Interval='95%'),GRAPH_DATA%>% filter(Group==0) %>% mutate(Interval='Median Prediction'))
	PLOT <- PLOT+geom_line(aes(x=Year,y=MIN,linetype=Interval,color=Interval),linewidth=LINE_WIDTH, data=CI_90)+geom_line(aes(x=Year,y=MAX,group=Interval,linetype=Interval,color=Interval),linewidth=LINE_WIDTH ,data=CI_90)+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash"))
	return(PLOT)
}
POP_DATA <- GET_DATA(RES,3)
POP_PLOT <- MAKE_GRAPH(POP_DATA)
POP_PLOT <- POP_PLOT+geom_line(data=HIST,aes(x=Year,y=Population),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2060, by = 10),2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 500))+ggtitle("Kemmerer & Diamondville, Population Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Population")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
png(paste0(SAVE_RES_LOC,"Population_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600)
POP_PLOT 
dev.off()

BIRTH_DATA <- GET_DATA(RES,4)
BIRTH_PLOT <- MAKE_GRAPH(BIRTH_DATA)
BIRTH_PLOT <- BIRTH_PLOT+geom_line(data=HIST,aes(x=Year,y=Births),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Birth Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Births")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
png(paste0(SAVE_RES_LOC,"Birth_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600)
	BIRTH_PLOT 
dev.off()

DEATH_DATA <- GET_DATA(RES,5) %>% filter(!is.na(MIN))
DEATH_PLOT <- MAKE_GRAPH(DEATH_DATA)
DEATH_PLOT <-	DEATH_PLOT+geom_line(data=HIST,aes(x=Year,y=Deaths),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(0, 35000, by = 10))+ggtitle("Kemmerer & Diamondville, Mortality Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Deaths")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
png(paste0(SAVE_RES_LOC,"Mortality_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600)
	DEATH_PLOT
dev.off()

MIGRATION_DATA <- GET_DATA(RES,6) %>% filter(!is.na(MIN))
MIGRATION_PLOT <- MAKE_GRAPH(MIGRATION_DATA)
MIGRATION_PLOT <- 	MIGRATION_PLOT+geom_line(data=HIST,aes(x=Year,y=Migration),color='black',linewidth=0.75)+ scale_x_continuous(breaks = c(seq(2010, 2060, by = 5),2065),limits=c(2009,2065))+ scale_y_continuous(breaks = seq(-1000, 1000, by = 50))+ggtitle("Kemmerer & Diamondville, Net Migration Forecast")+ expand_limits( y = 0)+labs(color = "Prediction Interval",linetype="Prediction Interval",y="Migration")+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())
png(paste0(SAVE_RES_LOC,"Migration_Fan_Chart.png"), width = 12, height = 8, units = "in", res = 600)
	MIGRATION_PLOT
dev.off()

#####Key year table
KEY <- RES %>% filter(Year %in% c(2029,2030,2035,2045,2055,2065))
 AVG_VALUES <- KEY %>% group_by(Year) %>% summarize(MED=median(Population),MEAN=mean(Population))
AVG_VALUES <- rbind(AVG_VALUES[,1:2]%>% rename(Value=MED) %>% mutate('Summary Stat.'="Median"),AVG_VALUES[,c(1,3)] %>% rename(Value=MEAN) %>% mutate('Summary Stat.'="Mean"))

HISTOGRAM <- ggplot(KEY, aes(x = Population,group=-Year,Color=Year,fill=Year)) + geom_histogram(alpha=0.3,bins=800)+geom_vline(data = AVG_VALUES, aes(xintercept = Value,group=`Summary Stat.`,color = `Summary Stat.`), size = 0.75)+scale_color_manual(values=c("red","black","black"))+ facet_grid(rows=vars(Year))+ scale_x_continuous(breaks = c(seq(0, 10000, by = 500)))+ theme_bw()+ theme(legend.position = "top",panel.grid.minor = element_blank())+ylab("Number of Simulation")+guides(fill= guide_legend(nrow = 1))
png(paste0(SAVE_RES_LOC,"Population_Histogram.png"), width = 8, height = 12, units = "in", res = 600)
	HISTOGRAM
dev.off()

#rm(KEY_YEARS)
POP_LEVELS <- seq(2000,6000,100)
YEARS <- c(2030,2035,2045,2055,2065)
for(i in YEARS ){
	KEY <- RES %>% filter(Year==i )  %>% pull(Population)
	ECDF <- ecdf(KEY)
	ECDF_VALUES <- ECDF(POP_LEVELS)
	if(!exists("KEY_YEARS")){KEY_YEARS<- ECDF_VALUES} else{KEY_YEARS<- cbind(KEY_YEARS,ECDF_VALUES)}

}
colnames(KEY_YEARS) <-  YEARS
rownames(KEY_YEARS) <- POP_LEVELS 
PLOT_GREEN <- "forestgreen"
PLOT_YELLOW <- "yellow"
PLOT_RED <- "red"
KEY_YEARS <- KEY_YEARS%>% as.data.frame
Capacity_Risk <- KEY_YEARS%>%  gt(rownames_to_stub = TRUE,caption="Year")  %>%  data_color(    fn = scales::col_numeric( palette = c(PLOT_RED, PLOT_YELLOW, PLOT_GREEN),  domain = c(0, 1) ) ) %>%  fmt_percent( decimals = 1,  drop_trailing_zeros = FALSE) %>% tab_stubhead(label =c("Capacity"))
gtsave( data = Capacity_Risk ,  filename = "./Results/Primary_Simulation_Results/Main_Results/Capacity_Table.html")
system("wkhtmltopdf --disable-smart-shrinking --no-stop-slow-scripts --enable-local-file-access --page-width 85mm --page-height 328mm ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table.html ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table.pdf")