#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)

#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.

quit(save = "no")