library(tidyverse)
	NUM_YEARS_PROJECTED <-  50 #How many years into the future should each Monte Carlo run project to. For example 25 years if starting from 2025 and ending in 2050.
YEARS <- 2023:(2023+NUM_YEARS_PROJECTED)
#Setup save results
	RES_DIR <- "./Results"
	RAW_SIM_FILE <- paste0(RES_DIR,"/Raw_Simulations.csv")
	PERCENTILE_DATA <- paste0(RES_DIR,"/Percentile_Clean_Results.csv")

###Process the simulations and save the main percentile results by year
	RES <- read_csv(RAW_SIM_FILE)
	GRAPH_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(RES %>% filter(Year==x) %>% pull(Population),c(0.025,0.05,0.1,0.25,0.4,0.5,0.6,0.75,0.9,0.95,0.975))}))  %>% as_tibble
	YEARS <- 2023:(2023+NUM_YEARS_PROJECTED)
	GRAPH_DATA$Year <- YEARS 
	FAN_DATA <- GRAPH_DATA
	GRAPH_DATA <- GRAPH_DATA %>% pivot_longer(cols=!Year,names_to=c("Percentile"),values_to="Population")
	write_csv(GRAPH_DATA,PERCENTILE_DATA)
	#Add historic
	MEDIAN_PRED <- 	GRAPH_DATA %>% filter(Percentile=='50%')
	GRAPH_DATA <- 	GRAPH_DATA %>% filter(Percentile!='50%')

	HIST <-	readRDS("Data/Cleaned_Data/Wyoming_County_Population.Rds") %>% filter(County=='Lincoln') %>% mutate(Percentile="Actual Population") %>% filter(Year>1930)
	ALPHA=0.2
	COLOR <- 'black'
GRAPH_DATA$Percentile <-	factor(GRAPH_DATA$Percentile,levels=rev(c('2.5%','5%','10%','25%','40%','60%','75%','90%','95%','97.5%')))
GRAPH <-	ggplot(data=GRAPH_DATA)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`2.5%`,ymax=`97.5%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`5%`,ymax=`95%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`10%`,ymax=`90%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`25%`,ymax=`75%`),alpha=ALPHA,fill=COLOR)+geom_ribbon(data=FAN_DATA,aes(x=Year,ymin=`40%`,ymax=`60%`),alpha=ALPHA,fill=COLOR)+geom_line(aes(x=Year,y=Population,group=Percentile,color=Percentile))+geom_line(data=HIST,aes(x=Year,y=Population),color='black',size=0.75)+geom_line(data=MEDIAN_PRED,aes(x=Year,y=Population),color='black',linetype=4,size=0.75)+ scale_x_continuous(breaks = c(seq(1940, 2070, by = 10)))+ scale_y_continuous(breaks = seq(0, 35000, by = 5000))+theme_bw()+ggtitle("Lincoln County, Wyoming Population Forecast")
GRAPH 
length(RES$SIM_ID %>% unique)
ggsave(Lincoln_Forecast.png",GRAPH)