#This is the analysis portion of the results vectors that get produced by the statelevel_modelMSI.R file. 

#Handle packages -- not all needed for this file alone.

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)

#Some key inputs
potentialBs = c(10,20,30,40,50,60,70,80,90,100,120,140,160,180,200,300,400,500,600,700)
boats.n.ij = readRDS(file="boats_adjSWF")

#We loop over the budget values we had...
for(budget in 1:length(potentialBs)) {
  
  #We need to generate the proper file name to append to the loading call below based on the current budget. 
  appendix = potentialBs[budget]
  
  #We first load in that result's vector of lakes chosen.
  sol1 = readRDS(paste0("statewide_model_results/lakes_selected_stateB", appendix))
  
  #We then just go into the boats.n.ij object and first slice out all the rows indicated by sol1--those were the lakes chosen.
  sol.slice1 = boats.n.ij[sol1, ]
  boats.sliced1 = boats.n.ij[-sol1, ]
  
  #We then also slice out the columns by those same values--this reflects the boats interfacing with the chosen lakes but going in the other direction.
  sol.slice2 = boats.sliced1[, sol1]
  boats.sliced2 = boats.sliced1[, -sol1]
  
  #Now, we can just sum up all the boats intercepted across slices 1 and 2.
  boats.stopped = sum(as.vector(sol.slice1), as.vector(sol.slice2))
  
  #Store the result (growing a vector, which is fine here.)
  if(budget == 1) {
    lakes_selected_state = boats.stopped
  } else {
    lakes_selected_state = c(lakes_selected_state, boats.stopped)
  }
  
}

#Save this to file in accordance with the outputs of the other modeling files.
saveRDS(lakes_selected_state, "lakes_selected_state")