library(tidyverse)
##Value from table C-2 of Holtec report, cost categories
	#"Initial Annual Estimated Costs and 2019 Constant Dollar Values for the Various Activities for the Proposed CISF and the No-Action Alternative"

#########
#CURRENT_YEAR <- 30 
#END_YEAR <- 40
#PHASE_SIZE <- 5000
#ST_CAP <- 8680
#TRANSPORT_COST_RATIO <- 155462880/5000 # or 269883561/8680 since transportation cost scale linearly
#VOLUME <- 10^5/2+480
#TRANS_CONSTRAINT <- (ST_CAP+5000*19)/19
#PHASES_ADDED <- max(ceiling((VOLUME-ST_CAP)/PHASE_SIZE),0) #Make sure added phases are not negative

OPTIMAL_COST <- function(PHASES_ADDED,VOLUME,CURRENT_YEAR,DISCOUNT=0.05,STARTING_VOLUME=8680,PHASE_SIZE=5000,END_YEAR=40,PHASE_CONST_COST=103399272,TRANSPORT_COST_RATIO=155462880/5000,TRANS_CONSTRAINT =(8680+5000*19)/19,DECOM_COST_PER_TON =24822656/5000,INFLATION_ADJUST=1.2874){
	ADDED_VOLUME=max(VOLUME-STARTING_VOLUME,0)
if(ADDED_VOLUME/PHASE_SIZE>PHASES_ADDED){stop("Not enough capacity for the requested volume")}
	###Construction Cost to cover supplied SNF volume
	CONSTRUCT_COST <- PHASES_ADDED*PHASE_CONST_COST
	#Cost to transport SNF from reactors to the CIFS in the current year and to repository at end of project period
	SHIPPING_TIME<- VOLUME/TRANS_CONSTRAINT 
	SHIPPING_YEARS <- ceiling(SHIPPING_TIME )
	if(SHIPPING_YEARS>END_YEAR-CURRENT_YEAR){stop("Not enough time to ship the requested SNF volume")}
	AT_CAPACITY_VOLUME <- TRANS_CONSTRAINT*floor(SHIPPING_TIME )
	UNDER_CAPACITY_VOLUME <- VOLUME-AT_CAPACITY_VOLUME
	SHIPPING_SCHEDULE_OUT <- rep(AT_CAPACITY_VOLUME,SHIPPING_YEARS)
	if(UNDER_CAPACITY_VOLUME!=0){SHIPPING_SCHEDULE_OUT[1] <- UNDER_CAPACITY_VOLUME}
	SHIPPING_SCHEDULE_OUT <- TRANSPORT_COST_RATIO*SHIPPING_SCHEDULE_OUT
	SHIPPING_SCHEDULE_IN <- rev(SHIPPING_SCHEDULE_OUT)
	SHIPPING_DISCOUNT <- (1/((1+DISCOUNT)^(1:SHIPPING_YEARS)))
	SHIPPING_IN_TOTAL_COST <- sum(SHIPPING_SCHEDULE_IN*SHIPPING_DISCOUNT)
	OUT_DISCOUNT <- 1/((1+DISCOUNT)^(END_YEAR-CURRENT_YEAR-SHIPPING_YEARS))
	SHIPPING_OUT_TOTAL_COST <- OUT_DISCOUNT*sum(SHIPPING_SCHEDULE_OUT*SHIPPING_DISCOUNT)
	SHIPPING_IN_TOTAL_COST +SHIPPING_OUT_TOTAL_COST +CONSTRUCT_COST 
	##Decommsioning cost: See Holtec Report, the NRC applies a adjustment factor for the larger capacity which we apply also to the small capacity this is in Section C-2 (just after the table and above section C-3)
	YEARS_UNTIL_DECOM <- END_YEAR-CURRENT_YEAR+1
	#Assume decom is paid by phase not average volume
	DECOM_TOTAL_COST <-(DECOM_COST_PER_TON*PHASES_ADDED*PHASE_SIZE )/((1+DISCOUNT)^YEARS_UNTIL_DECOM)
	TOTAL_COST <- CONSTRUCT_COST+SHIPPING_IN_TOTAL_COST+SHIPPING_OUT_TOTAL_COST+ DECOM_TOTAL_COST 
	TOTAL_COST <- TOTAL_COST*INFLATION_ADJUST
	return(TOTAL_COST)
}
OPTIMAL_COST(20,10^4,20)

CHECK_FEASIBLE_SHIPPING <- function(VOLUME,ST_YEAR){
RESULT <- try(OPTIMAL_COST(10^5,VOLUME,ST_YEAR),silent=TRUE)
return(class(RESULT)!="try-error")
}
FIND_FEASIBLE_LIMIT <- function(STARTING_TIME){
	ST_BOUND <-0
	END_BOUND <-140*10^3
	MID <- ceiling((END_BOUND+ST_BOUND)/2)
	while(END_BOUND-ST_BOUND>10){
		CHECK <-CHECK_FEASIBLE_SHIPPING(MID,STARTING_TIME)
		if(CHECK){ST_BOUND <- MID} else{END_BOUND <- MID}
		MID <- ceiling((END_BOUND+ST_BOUND)/2)
	}
	MAX_CAPACITY <- ST_BOUND+max(which(sapply(ST_BOUND:END_BOUND,CHECK_FEASIBLE_SHIPPING,ST_YEAR=STARTING_TIME)))
	return(MAX_CAPACITY)
}
SHIPPING_CAPACITY_LIMITS <- cbind(1:40,sapply(1:40,FIND_FEASIBLE_LIMIT)) %>% as_tibble %>% rename(Year=V1,Max_Capacity=V2) 
SHIPPING_CAPACITY_LIMITS  %>% print(n=100)
REACTOR_VALUES <- readRDS("Data/Cleaned_Data/Reactor_Values.Rds")
C_VALUES <- REACTOR_VALUES  %>% filter(Year==2026+40,Discount==0.05)
C_VALUES
TBL <- CURRENT
CIFS_SIZE <- ST_CAP
MAX_REV <- function(TBL,CIFS_SIZE){
#		TBL <- TBL %>% filter(Year==YEAR,Discount==DISCOUNT,Revenue/Total_Tons>SHIPPING_COST)
		REV <- TBL %>% pull(Revenue)
		VOL <- TBL %>% pull(Total_Tons)
		RES <-	lp(direction = "max", objective.in = REV, const.mat = matrix(VOL, nrow = 1,byrow=TRUE),const.dir = "<=", const.rhs = CIFS_SIZE, all.bin = TRUE)
		 return(RES)
	}
ST_CAP <- 8680
PHASE_SIZE <- 5000

RES <-	
	names(RES)
PROFIT_EST <- function(ADDED_PHASES,ST_YEAR,YEARS_AHEAD,DATA=REACTOR_VALUES,DISCOUNT_RATE=0.05,ST_CAP=8680,PHASE_SIZE=5000){
	CURRENT <-  DATA%>% filter(Year==ST_YEAR+YEARS_AHEAD,Discount==DISCOUNT_RATE)
	RES <- MAX_REV(CURRENT,ST_CAP+PHASE_SIZE*ADDED_PHASES)
	REVENUE <- RES$objval
	TONS_STORED <- 	sum(CURRENT$Total_Tons*RES$solution)
	PROFIT <- rbind(REVENUE,OPTIMAL_COST(ADDED_PHASES,TONS_STORED,YEARS_AHEAD))
	return(PROFIT)
}
MAX_REV(CURRENT,ST_CAP+PHASE_SIZE*1)

OPTIMAL_COST(0,TONS_STORED,YEARS_AHEAD)
t(sapply(0:3,PROFIT_EST,ST_YEAR=2026,YEARS_AHEAD=15)/10^6) %>% as_tibble  %>% rename(Rev=V1,Cost=V2) %>% mutate(Profit=Rev-Cost)

RES$solution
names(RES)
length(RES$objective)

CURRENT[RES[[9]],]
CURRENT[-RES[[9]],]

REACTOR_VALUES %>% filter(Year=[RES[[9]],]