From cacb64a6578703e1fe5051381deec990cd0378e2 Mon Sep 17 00:00:00 2001 From: Alex Gebben Work Date: Mon, 22 Dec 2025 18:08:55 -0700 Subject: [PATCH] Fixed missing permanent jobs added --- 1B_Run_Full_Simulation.r | 37 +++++++----- TESTING.r | 124 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 147 insertions(+), 14 deletions(-) create mode 100644 TESTING.r diff --git a/1B_Run_Full_Simulation.r b/1B_Run_Full_Simulation.r index cce55e0..45ccde3 100644 --- a/1B_Run_Full_Simulation.r +++ b/1B_Run_Full_Simulation.r @@ -95,9 +95,10 @@ return(list(DEMO,BIRTH_DATA,c(TOTAL_POP,TOTAL_BIRTHS,TOTAL_DEATHS,TOTAL_MIGRATIO MIGRATION_ARIMA_MODEL <- MIGRATION_ARIMA #DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL,OPERATOR_TOTAL,CONSTRUCTION_TOTAL,MIGRATION_MULTIPLIERS -CONSTRUCTION_MIGRATION <- CONSTRUCTION_LIN_MIGRATION -MIGRATION_MULTIPLIERS <- INDUCED_MIGRATION_MULTIPLIERS -OPERATOR_MIGRATION <- OPERATOR_LIN_MIGRATION +#CONSTRUCTION_MIGRATION <- CONSTRUCTION_LIN_MIGRATION +#MIGRATION_MULTIPLIERS <- INDUCED_MIGRATION_MULTIPLIERS +#OPERATOR_MIGRATION <- OPERATOR_LIN_MIGRATION +MIGRATION_ARIMA=MIGRATION_ARIMA_MODEL;OPERATOR_MIGRATION=OPERATOR_LIN_MIGRATION;CONSTRUCTION_MIGRATION=CONSTRUCTION_LIN_MIGRATION;MIGRATION_MULTIPLIERS=INDUCED_MIGRATION_MULTIPLIERS SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL,OPERATOR_MIGRATION,CONSTRUCTION_MIGRATION,MIGRATION_MULTIPLIERS ){ @@ -109,11 +110,17 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL CONSTRUCTION_MIGRATION <- LOCAL_WORK_ADJ(CONSTRUCTION_MIGRATION,0.41) #Assume between 41%-100% operators live in Kemmerer CONSTRUCTION_MIGRATION[7] <- CONSTRUCTION_MIGRATION[7] - sum(CONSTRUCTION_MIGRATION ) CONSTRUCTION_POPULATION_ADDED <- cumsum(CONSTRUCTION_MIGRATION) - ADDED_FROM_BASELINE <- runif(1) #The percentage of the possible growth in industries like restaurants to add compared to the Kemmerer IMPLAN model which understates possible structural growth + ###########TEST + GROWTH_LEVEL <- sample(c("VERY_LOW","LOW","HIGH","HIGH"),1) +if(GROWTH_LEVEL=="VERY_LOW"){ADDED_FROM_BASELINE <- 0;MIGRATION_ADJUST <- 0} else if(GROWTH_LEVEL=="LOW"){ADDED_FROM_BASELINE <- runif(1,0,0.5);MIGRATION_ADJUST <- 0}else{ADDED_FROM_BASELINE <- runif(1,0.5,1);MIGRATION_ADJUST <-(ADDED_FROM_BASELINE-0.5)/0.5} + MIGRATION_SHIFT <- 55*(1-MIGRATION_ADJUST) + + ############## + #ADDED_FROM_BASELINE <- runif(1) #The percentage of the possible growth in industries like restaurants to add compared to the Kemmerer IMPLAN model which understates possible structural growth PERMANENT_TERRAPOWER_MIGRATION <- INDUCED_SIMULATION(CONSTRUCTION_MIGRATION,OPERATOR_MIGRATION,MIGRATION_MULTIPLIERS,ADDED_FROM_BASELINE)+OPERATOR_MIGRATION - -###############NOTE NEED TO USE THIS AT END TO ADJUST THE RESULTS WHILE LEAVING THE DEMOGRAPHIC MATRIX - TERRA_POWER_EFFECT[1:7] <- TERRA_POWER_EFFECT[1:7]+CONSTRUCTION_MIGRATION + PERMANENT_TERRAPOWER_MIGRATION <- round(PERMANENT_TERRAPOWER_MIGRATION*2.314328) #Households per job in IMPLAN + ###############NOTE NEED TO USE THIS AT END TO ADJUST THE RESULTS WHILE LEAVING THE DEMOGRAPHIC MATRIX + TERRA_POWER_EFFECT[1:7] <- TERRA_POWER_EFFECT[1:7]+CONSTRUCTION_MIGRATION +PERMANENT_TERRAPOWER_MIGRATION #Dry Pine Migration DRY_PINE <- round(runif(1,0.15,0.45)*c(126,176,-302)*1.15) DRY_PINE[3] <- DRY_PINE[3]-sum(DRY_PINE) @@ -123,7 +130,6 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL if(runif(1)>0.67){ MW_CAPACITY <- runif(1,20,120) PERM_WORKERS <- MW_CAPACITY/8+30 #Seems to be a good ration 32.5 workers at 20 capacity and 45 at max this is a heurstic I came up with - OTHER_PROJECT_EFFECTS_PERM DATA_CENTER_START_YEAR <- which.max(runif(YEARS_AHEAD-6))+3 PERM_DATA_CENTER_OP <- PERM_WORKERS*(1+(2.515334+19.52286*ADDED_FROM_BASELINE)/40) #Add induced #Estimate construction workers on data center @@ -139,18 +145,21 @@ SINGLE_SIM <- function(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA_MODEL OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]<- OTHER_PROJECT_EFFECTS_TEMP_POP_ADDED[(DATA_CENTER_START_YEAR+3):(DATA_CENTER_START_YEAR+5)]+TEMP_DATA_MIGRATION OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+6] <- OTHER_PROJECT_EFFECTS_PERM[DATA_CENTER_START_YEAR+3]+PERM_DATA_CENTER_OP } - MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))+runif(1,-55,0)+TERRA_POWER_EFFECT+OTHER_PROJECT_EFFECTS_PERM) + MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))-MIGRATION_SHIFT+TERRA_POWER_EFFECT+OTHER_PROJECT_EFFECTS_PERM) + #MIGRATION_SIM_VALUES <- round(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA_MODEL))+-55+TERRA_POWER_EFFECT+OTHER_PROJECT_EFFECTS_PERM) + #The runif applies a downshift ranging from the historic decline rate all the way to the Lincoln rate applied in the model - FINAL_REPORT_VALUES <- matrix(NA,ncol=6,nrow=YEARS_AHEAD) - colnames(FINAL_REPORT_VALUES ) <- c("Sim_UUID","Year","Population","Births","Deaths","Net_Migration") + FINAL_REPORT_VALUES <- matrix(NA,ncol=7,nrow=YEARS_AHEAD) + colnames(FINAL_REPORT_VALUES ) <- c("Sim_UUID","Year","Population","Births","Deaths","Net_Migration","Growth_Rate") FINAL_REPORT_VALUES[,1] <- UUIDgenerate() + FINAL_REPORT_VALUES[,7] <- GROWTH_LEVEL for(i in 1:YEARS_AHEAD){ C_YEAR <- ST_YEAR+i-1 C_RES <-SINGLE_YEAR_SIM(DEMO,BIRTH_DATA,i,SIMULATE_MORTALITY_RATE_TRENDS(),MIGRATION_SIM_VALUES[i]) DEMO <- C_RES[[1]] BIRTH_DATA <- C_RES[[2]] - FINAL_REPORT_VALUES[i,-1] <- c(C_YEAR,C_RES[[3]]) + FINAL_REPORT_VALUES[i,c(-1,-7)] <- c(C_YEAR,C_RES[[3]]) } FINAL_REPORT_VALUES[1:7,3] <- as.numeric(FINAL_REPORT_VALUES[1:7,3])+CONSTRUCTION_POPULATION_ADDED FINAL_REPORT_VALUES[1:7,6] <- as.numeric(FINAL_REPORT_VALUES[1:7,6]) +CONSTRUCTION_MIGRATION @@ -180,11 +189,11 @@ CONSTRUCTION_LIN_MIGRATION <- CONSTRUCTION %>% pull("Construction_Emp_Migrated") MIGRATION_MATRIX <- do.call(cbind, mclapply(1:BATCH_SIZE,function(x)(as.vector(simulate(nsim=YEARS_AHEAD,MIGRATION_ARIMA) )),mc.cores = detectCores()-1)) # rownames(MIGRATION_MATRIX) <- ST_YEAR:(ST_YEAR+YEARS_AHEAD-1) # colnames(MIGRATION_MATRIX) <- 1:BATCH_SIZE -#SINGLE_SIM(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA,OPERATOR_LIN_MIGRATION,CONSTRUCTION_LIN_MIGRATION,INDUCED_MIGRATION_MULTIPLIERS) +SINGLE_SIM(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA,OPERATOR_LIN_MIGRATION,CONSTRUCTION_LIN_MIGRATION,INDUCED_MIGRATION_MULTIPLIERS) try(RES <- do.call(rbind,mclapply(1:BATCH_SIZE,function(x){SINGLE_SIM(DEMO,BIRTH_DATA,ST_YEAR,YEARS_AHEAD,MIGRATION_ARIMA,OPERATOR_LIN_MIGRATION,CONSTRUCTION_LIN_MIGRATION,INDUCED_MIGRATION_MULTIPLIERS)},mc.cores=NCORES))) if(exists("RES")){ RES <- as.data.frame(RES) - RES[,-1] <- as.numeric(as.matrix(RES[,-1])) + RES[,c(-1,-7)] <- as.numeric(as.matrix(RES[,c(-1,-7)])) if(NEW_RES_FILE & i==1){write_csv(RES,SIM_RES_FILE)}else {write_csv(RES,SIM_RES_FILE,col_names=FALSE,append=TRUE)} rm(RES) } diff --git a/TESTING.r b/TESTING.r new file mode 100644 index 0000000..35d694e --- /dev/null +++ b/TESTING.r @@ -0,0 +1,124 @@ +library(tidyverse) +library(gt) #For nice color coded capacity limits table. +source("Scripts/Load_Custom_Functions/Fan_Chart_Creation_Functions.r") #Functions created to make nice graphs + 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_With_Data_Center_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 +#POP_DATA <- GET_DATA(RES,3) +#POP_PLOT <- MAKE_GRAPH(POP_DATA) +########TEST +SIM_START=2025 + LAST_HISTORIC <- SIM_START-1 + YEARS <- min(RES$Year,na.rm=TRUE):max(RES$Year,na.rm=TRUE) + LEVELS <- seq(0.00,1,by=0.005) + GROUPS <- floor(length(LEVELS)/2) + COL_NUM <-3 + RES %>% tail +#C_DATA<- RES %>% filter(Growth_Rate=='HIGH') +#C_DATA <- RES %>% filter(Growth_Rate!='HIGH') + C_DATA <- RES + +nrow(RES) +FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((C_DATA %>% filter(Year==x))[,COL_NUM])),LEVELS)})) %>% as_tibble %>% mutate(Year=YEARS) +FAN_DATA <- FAN_DATA %>% pivot_longer(!Year,names_to="Percentile",values_to="Population") %>% group_by(Year) %>% mutate(DIFF=c(NA,diff(Population))) +GRAPH_DATA <- FAN_DATA %>% group_by(Year) %>% mutate(MIN_RANGE=lag(parse_number(Percentile)/100),MAX_RANGE=parse_number(Percentile)/100,MIN_POP=lag(Population),MAX_POP=Population) %>% filter(!is.na(DIFF)) +GRAPH_DATA <- GRAPH_DATA %>% mutate(MAX_POP=ifelse(MAX_RANGE!=1,MAX_POP-0.001,MAX_POP)) + CI_80 <- GRAPH_DATA%>% filter(Percentile %in%c('10.0%','90.0%')) %>% group_by(Year) %>% summarize(MIN_POP=min(MIN_POP),MAX_POP=max(MAX_POP),Interval='80%') + CI_95 <- GRAPH_DATA%>% filter(Percentile %in%c('2.5%','97.5%')) %>% group_by(Year) %>% summarize(MIN_POP=min(MIN_POP),MAX_POP=max(MAX_POP),Interval='95%') +MEDIAN <- GRAPH_DATA%>% filter(Percentile %in%c('50.0%'))%>% group_by(Year) %>% summarize(MIN_POP=min(MIN_POP),MAX_POP=max(MAX_POP),Interval='Median') +RANGE_DATA <- rbind(MEDIAN ,CI_80,CI_95) +ggplot(GRAPH_DATA %>% filter(MIN_RANGE>=0.015,MAX_RANGE<=0.985) )+geom_ribbon(aes(group=Percentile,x=Year,ymin=MIN_POP,ymax=MAX_POP,alpha=-DIFF),fill='cadetblue')+geom_line(data=GRAPH_DATA %>% filter(Percentile=='50%'),aes(x=Year,y=MIN_POP))+geom_line(data=RANGE_DATA,aes(color=Interval,linetype=Interval,x=Year,y=MIN_POP))+geom_line(data=RANGE_DATA,aes(color=Interval,linetype=Interval,x=Year,y=MAX_POP))+scale_color_manual(values=c("grey50","grey80","black"))+scale_linetype_manual(values = c("solid","solid","dotdash"))+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())+scale_alpha(range = c(0, 1),guide="none") +C_DATA + + geom_ + +RANGE_DATA + + +QUANT <- quantile(as.numeric(t((RES %>% filter(Year==2025))[,3])),LEVELS) %>% as.numeric +RANGE <- QUANT %>% diff +RANGE +RANGE <- RANGE[-1] +RANGE <- RANGE[-length(RANGE)] +RANGE <- 1-RANGE/max(RANGE) +RANGE +scale(RANGE) +FAN_DATA <- do.call(rbind,lapply(YEARS,function(x){quantile(as.numeric(t((RES %>% filter(Year==x))[,COL_NUM])),LEVELS)})) %>% as_tibble %>% mutate(Year=YEARS) + FAN_DATA <- rbind(FAN_DATA[1,],FAN_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_Main_Results.png"), width = 12, height = 8, units = "in", res = 600) +POP_PLOT +dev.off() + +#Population Summary table for the key years +KEY_YEARS <- c(2027,2030,2035,2045,2055,2065) +KEY_YEAR_SUMMARY <- POP_DATA %>% filter(Year %in% KEY_YEARS,Group==5) %>% select(-Group) %>% rename(CI_95_Lower=MIN,CI_95_Upper=MAX) %>% left_join(POP_DATA %>% filter(Year %in% KEY_YEARS,Group==0) %>% select(-Group) %>% select(Year,Median=MIN)) %>% select(Year,CI_95_Lower,Median,CI_95_Upper) +KEY_YEAR_SUMMARY <- round(KEY_YEAR_SUMMARY ) +KEY_YEAR_SUMMARY <- t(KEY_YEAR_SUMMARY ) +write.table(KEY_YEAR_SUMMARY,paste0(SAVE_RES_LOC,"Key_Year_Summary_Main_Results.csv"), sep=",",col.names=FALSE) + + +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_Main_Results.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_Main_Results.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_Main_Results.png"), width = 12, height = 8, units = "in", res = 600) + MIGRATION_PLOT +dev.off() + +#####Key year table + +YEARS <- c(2027,2030,2035,2045,2055,2065) +KEY_YEAR_DATA <- RES %>% filter(Year %in% YEARS) + AVG_VALUES <- KEY_YEAR_DATA %>% 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_YEAR_DATA, 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.`), linewidth= 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_Main_Results.png"), width = 8, height = 12, units = "in", res = 600) + HISTOGRAM +dev.off() +#rm(KEY_YEARS) +POP_LEVELS <- seq(2000,6000,100) +YEARS +if(exists("KEY_YEARS")){rm(KEY_YEARS)} +for(i in YEARS ){ + CURRENT_YEAR_DATA <- RES %>% filter(Year==i ) %>% pull(Population) + ECDF <- ecdf(CURRENT_YEAR_DATA) + 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_Main_Results.html") +system("wkhtmltopdf --disable-smart-shrinking --no-stop-slow-scripts --enable-local-file-access --page-width 99mm --page-height 328mm ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.html ./Results/Primary_Simulation_Results/Main_Results/Capacity_Table_Main_Results.pdf") +