#THIS FILE HOLDS CODE FOR CONSTRUCTING A STATE-LEVEL MODEL WHEREIN A GENERIC STATE ACTOR ALLOCATES  INSPECTION STATIONS SO AS TO INTERCEPT THE MAXIMUM NUMBER OF RISKY BOAT TRIPS, IRRESPECTIVE OF JURISDICTIONAL BOUNDARIES. 

#ON A POWERFUL COMPUTER, THIS TAKES XX TO FULLY COMPLETE. YOU HAVE BEEN WARNED :). 

#STEP 0: LOAD PACKAGES. NOTE--NOT ALL OF THESE MAY STILL BE NEEDED, BUT THIS HASN'T BEEN VERIFIED. 
list.of.packages = c("pacman", "dplyr", "ROI", "ompr", "ompr.roi", "slam", "Rsymphony",
                     "ROI.plugin.symphony", "tidyr", "lattice", "Matrix", "igraph")
new.packages = list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if (length(new.packages) > 0) {
  r=getOption("repos")
  r["CRAN"] = "http://cran.rstudio.com/"
  options(repos=r)
  install.packages(new.packages)
}
pacman::p_load(dplyr, ROI, ompr, ompr.roi, slam, Rsymphony, ROI.plugin.symphony, tidyr, lattice, Matrix, igraph)

args = commandArgs(trailingOnly = TRUE)
iteration = as.integer(args[1])

#STEP 1.1: LOAD IN OUR NETWORK ADJACENCY MATRIX, OUR FIRST KEY INPUT OBJECT. THIS INCLUDES LAKES AS BOTH ROWS AND COLUMNS (AS "FROM" AND AS "TO" LAKES) WITH CELLS HOLDING THE NUMBERS OF BOATS MOVING FROM-->TO YEARLY. BY THIS POINT, EDGES THAT CANNOT RESULT IN ANY NEW INVASIVES (BIJ == 0) HAVE ALREADY BEEN REMOVED.
boats.n.ij = readRDS(file="boats_adjSWF")
#STEP 1.2: GET THE FIRST DIMENSION OF THIS TABLE--THIS IS THE NUMBER OF UNIQUE LAKES IN OUR SYSTEM.
n.ij = dim(boats.n.ij)[1] 

#STEP 2: SET UP THE MODEL OBJECT FOR THIS MODEL. THIS INCLUDES HAVING A BINARY DECISION VARIABLE X FOR WHICH LAKES ARE BEING SELECTED AND A BINARY DECISION MATRIX VARIABLE U WHICH CORRESPONDS TO WHICH CELLS IN OUR BOATS.N.IJ MATRIX ARE BEING INTERCEPTED BY INSPECTIONS. THEN, WE SET OUR OBJECTIVE AS THE MAX OF U*BOATS.N.IJ WHEREVER BOATS.N.IJ IS NON-ZERO. LASTLY, WE CONSTRAIN SUCH THAT WE CANNOT BE INTERCEPTING BOATS ALONG AN EDGE (A CELL IN BOATS.N.IJ) WITHOUT INSPECTING AT LEAST THE LAKE REPRESENTED BY EITHER THE ROW OR COLUMN. 

#***WARNING--THIS IS A _BIG_ MODEL, SO THIS PART TAKES A GOOD LONG TIME!
# start_time <- Sys.time() #Get a start time.
# mod.obj = MIPModel() %>%
#   add_variable(x[i], type = "binary", i = 1:n.ij) %>%
#   add_variable(u[i, j], type = "binary", i = 1:n.ij, j = 1:n.ij, boats.n.ij[i,j] > 0) %>%
#   set_objective(sum_over(boats.n.ij[i,j] * u[i,j], i=1:n.ij, j=1:n.ij, boats.n.ij[i,j] >0), "max") %>%
#   add_constraint(u[i,j] <= x[i] + x[j], i = 1:n.ij, j = 1:n.ij, boats.n.ij[i,j] > 0)
# end_time <- Sys.time() #Get an end time.
# print(end_time - start_time)
# 
# saveRDS(mod.obj, file = "mod.obj")

mod.obj = readRDS(file = "mod.obj") #For shortcutting.

#STEP 3 -- SET UP A STORAGE OBJECT, LAKES_SELECTED_STATE, TO HOLD THE SOLUTIONS GENERATED. 
lakes_selected_state = list()

#STEP 4 -- ITERATE OVER A SET OF DIFFERENT STATEWIDE BUDGET VALUES, B, SOLVING THE MODEL FOR EACH B.
potentialBs = c(10,20,30,40,50,60,70,80,90,100,120,140,160,180,200,300,400,500,600,700)
  
B = potentialBs[iteration] #Set this B for this rep.

#Just change the one constraint.
mod.obj.this = mod.obj %>% 
  add_constraint(sum_over(x[i], i = 1:n.ij) <= B)

#STEP 5 -- SOLVE THIS ITERATION OF THE MODEL, USING A TIME LIMIT OF 4 HOURS AND A GAP LIMIT OF 1%, WHICHEVER IS HIT FIRST.
start_time <- Sys.time() #Get a start time.

mod.solved = mod.obj.this %>% 
  solve_model(with_ROI(solver = "symphony", verbosity = -1, time_limit=14400, gap_limit=1))
end_time <- Sys.time() #Get an end time.

print(end_time - start_time)

solution2 = get_solution(mod.solved, x[i])
solution3 = solution2[solution2$value == 1,]

chosen.xs = solution3$i
print(chosen.xs)
lakes_selected_state = chosen.xs

#Save result
saveRDS(lakes_selected_state, file = paste0("lakes_selected_stateB", B))

#}
# 
#Leave R.
quit(save = "no")