R/smart_class.R

Defines functions FUN

#### SmartProject#############################################
#' SmartProject Class
#'
#' The \code{SmartProject} class implements the main class of
#' \pkg{smartR} package.
#'
#' @docType class
#'
#' @usage NULL
#'
#' @keywords data
#' @return Object of \code{\link{R6Class}} with attributes and methods to fullfill
#' a complete analisys with the SMART approach.
#'
#' @format \code{\link{R6Class}} object.
#'
#' @field rawDataSurvey Stores the raw survey data after being populated by \code{loadSurveyLFD()} method.
#' @field yearInSurvey Stores the distinct years in the \code{rawDataSurvey} time-serie.
#' @field specieInSurvey Stores the distinct species in the \code{rawDataSurvey} time-serie.
#' @field surveyBySpecie Stores a list of \code{\link{SurveyBySpecie}} objects, one for each species in the time-series.

#' @field rawDataFishery Stores the raw fishery data as is in the provided csv file. The attribute is populated by \code{loadFisheryLFD()} method.
#' @field yearInFishery Stores the distinct years in the \code{rawDataFishery} time-serie.
#' @field specieInFishery Stores the distinct species in the \code{rawDataFishery} time-serie.
#' @field fisheryBySpecie Stores a list of \code{\link{FisheryBySpecie}} objects, one for each species in the time-series.
#'
#' @field gooLstCoho Stores a list of plots of species cohorts spatial distribution .
#'
#' @field sampMap Stores the \code{\link[=SampleMap]{environment}} object.
#' @field fleet Stores the \code{\link{FishFleet}} object.
#'
#' @field simProd Stores the simulated pattern of production.
#' @field simEffo Stores the simulated pattern of effort.
#' @field simBanFG Stores a vector of fishable/banned fishing grounds.
#' @field simSpatialCost Stores the simulated pattern of spatial costs.
#' @field simEffortCost Stores the simulated pattern of effort costs.
#' @field simProdCost Stores the simulated pattern of production costs.
#' @field simTotalCost Stores the simulated pattern of total costs.
#' @field simRevenue Stores the simulated pattern of revenue by species and fishing ground.
#' @field simTotalRevenue Stores the simulated pattern of total revenues.
#' @field simCostRevenue Stores the simulated pattern of costs and revenues.
#' @field simResPlot Stores the plots with the simulation' results.

#' @field outGmat Stores the evolution of gains during the simulation.
#' @field outOptimEffo Stores the resulting pattern of effort.
#' @field outWeiProp Stores the annual proportion of fish by cohort and fishing ground.
#' @field outWeiPropQ Stores the seasonal proportion of fish by cohort and fishing ground.
#'
#'
#' @section Methods:
#' \describe{
#'   \item{\code{setCostInput()}}{This method is used to setup the required data for costs computation}
#'   \item{\code{setInProduction()}}{This method is used to setup the required data for production costs computation}
#'   \item{\code{setDaysAtSea()}}{This method is used to compute the number of Days at Sea of each vessel}
#'   \item{\code{setEffortIndex()}}{This method is used to compute the value of the Effort Index}
#'   \item{\code{setProductionIndex()}}{This method is used to compute the value of the Production Index}
#'   \item{\code{getHarbFgDist()}}{This method is used to compute the weighted average distance of every fishing ground to each harbour}
#'   \item{\code{setFgWeigDist()}}{This method is used as helper function to get the weighted average distance between fishing ground and harbours}
#'   \item{\code{setRegHarbBox()}}{This method is used to compute the distance of each harbour to every fishing ground centroid}
#'   \item{\code{loadSurveyLFD(csv_path)}}{This method is used to load the raw survey LFD data from a  csv file}
#'   \item{\code{loadFisheryLFD(csv_path)}}{This method is used to load the raw fishery LFD data from a  csv file}
#'   \item{\code{setYearSurvey()}}{This method is used to store the distinct year in the survey time-series}
#'   \item{\code{setYearFishery()}}{TThis method is used to store the distinct year in the fishery time-series}
#'   \item{\code{loadMap(map_path)}}{This method is used to load the Environmental Grid and initialize the Environment object}
#'   \item{\code{createFleet()}}{This method is used to initialize the Fleet object}
#'   \item{\code{setSpecieSurvey()}}{This method is used to store the distinct species in the survey dataset}
#'   \item{\code{setSpecieFishery()}}{This method is used to store the distinct species in the fishery dataset}
#'   \item{\code{splitSpecieSurvey()}}{This method is used to split the survey dataset by species}
#'   \item{\code{splitSpecieFishery()}}{This method is used to split the fishery dataset by species}
#'   \item{\code{addSpecieSurvey(sing_spe)}}{This method is used to initialize a new \code{surveyBySpecie} object}
#'   \item{\code{addSpecieFishery(sing_spe)}}{This method is used to initialize a new \code{fisheryBySpecie} object}
#'   \item{\code{setSpreaFishery()}}{This method is used to prepare the fishery LFD data for MCMC analysis}
#'   \item{\code{setSpatFishery()}}{This method is used to setup the plot with the spatial distribution of the fishery dataset}
#'   \item{\code{setSpreaSurvey()}}{This method is used to prepare the survey LFD data for MCMC analysis}
#'   \item{\code{setSpatSurvey()}}{This method is used to setup the plot with the spatial distribution of the survey dataset}
#'   \item{\code{setDepthSurvey()}}{This method is used to assign the depth of each survey tow}
#'   \item{\code{setStratumSurvey()}}{This method is used to assign a depth stratum to each survey tow}
#'   \item{\code{setAbuAvgAll()}}{This method is used to compute the spiecies abundances at each survey stratum}
#'   \item{\code{setMeditsIndex()}}{This method is used to compute the MEDITS index}
#'   \item{\code{setStrataAbu()}}{This method is used to compute the weighted number of individuals of each size in every stratum}
#'   \item{\code{loadFleeEffoDbs(effort_path, met_nam, onBox = TRUE, perOnBox = 1)}}{This method is used to extract the vms data from one or more vmsbase DB}
#'   \item{\code{ggplotRawPoints(year)}}{This method is used to plot the raw vms points}
#'   \item{\code{ggplotFgWeigDists()}}{This method is used to plot the weighted average distance between harbours and fishing grounds}
#'   \item{\code{setAvailData()}}{This method is used to gather the required data for the spatial clustering}
#'   \item{\code{predictProduction(species)}}{This method is used to compute the estimated production}
#'   \item{\code{simProdAll(selRow = numeric(0))}}{This method is used to compute the simulated production}
#'   \item{\code{genSimEffo(method = "flat", selRow = numeric(0), areaBan = numeric(0))}}{This method is used to create a simulated pattern of effort}
#'   \item{\code{getSimSpatialCost()}}{This method is used to compute the simulated spatial costs}
#'   \item{\code{getSimEffortCost()}}{This method is used to compute the simulated effort costs}
#'   \item{\code{getSimProdCost()}}{This method is used to compute the simulated production costs}
#'   \item{\code{getSimTotalCost()}}{This method is used to collect all the simulated costs}
#'   \item{\code{getSimRevenue(selRow = numeric(0), timeScale = "Year")}}{This method is used to compute the simulated revenues}
#'   \item{\code{getLWstat()}}{This method is used to compute the length/weight statistics for each fishing ground}
#'   \item{\code{simulateFishery(thr0 = 100, effoBan = numeric(0), timeStep = "Year")}}{This method is used to simulate one year of fishing}
#'   \item{\code{setSimResults()}}{This method is used to store the results of a simulation run}
#'   \item{\code{ggplotFishingPoints(year)}}{This method is used to plot the fishing points}
#'   \item{\code{setCellPoin()}}{This method is used assign a cell to each vms point}
#'   \item{\code{setTrackHarb()}}{This method is used to assign the harbour to each fishing trip}
#'   \item{\code{setFishGround(numCut)}}{This method is used to setup the fishing ground configuration}
#'   \item{\code{addFg2Fishery()}}{This method is used to add the fishing ground information to each fishery data point}
#'   \item{\code{addFg2Survey()}}{This method is used to add the fishing ground information to each survey data point}
#'   \item{\code{setWeekEffoMatrCell()}}{This method is used to combine the raw effort points in weekly aggregated effort by cell}
#'   \item{\code{setWeekEffoMatrGround()}}{This method is used to combine the raw effort points in weekly aggregated effort by fishing ground}
#'   \item{\code{ggplotGridEffort(year)}}{This method is used to plot the gridded fishing effort}
#'   \item{\code{getNnlsModel(species, minobs, thr_r2)}}{This method is used to compute the coefficients of the NNLS model}
#'   \item{\code{cohoDisPlot(whoSpe, whoCoh, whiYea, interp)}}{This method is used to store the spatial distribution of the species by cohort}
#'   }
#' @examples
#' # Initialize SmartProject
#' yourSmartRstudy <- SmartProject$new()
#' 
#' # Initialize fleet object
#' yourSmartRstudy$createFleet()
#' 
#' 
#' ######################
#' ## Environment Data ##
#' ######################
#' 
#' # Locate the example environment asset' file
#' envAssetPath <- system.file("extdata/mapAsset.RDS", package = "smartR")
#' 
#' # Load environment asset' data
#' yourSmartRstudy$importEnv(readRDS(envAssetPath))
#' 
#' # Setup case study' map
#' yourSmartRstudy$sampMap$getGooMap()
#' yourSmartRstudy$sampMap$setGooGrid()
#' yourSmartRstudy$sampMap$setGooBbox()
#' yourSmartRstudy$sampMap$setGgDepth()
#' yourSmartRstudy$sampMap$setGgBioDF()
#' # View case study' grid
#' print(yourSmartRstudy$sampMap$gooGrid)
#' 
#' 
#' ################
#' ## Fleet Data ##
#' ################
#' 
#' # Locate the example fleet asset' file
#' effAssetPath <- system.file("extdata/effAsset.RDS", package = "smartR")
#' 
#' # Load fleet asset' data
#' yourSmartRstudy$fleet$rawEffort <- readRDS(effAssetPath)
#' 
#' # Setup fishing vessel ids
#' yourSmartRstudy$fleet$setEffortIds()
#' 
#' # View speed distribution to setup fishing point filter
#' yourSmartRstudy$fleet$plotSpeedDepth(
#' which_year = "2012",
#' speed_range = c(2, 8),
#' depth_range = c(-20, -600)
#' )
#' 
#' # Setup fishing points' filter
#' yourSmartRstudy$fleet$setFishPoinPara(
#' speed_range = c(2, 8),
#' depth_range = c(-20, -600)
#' )
#' 
#' # Compute fishing points
#' yourSmartRstudy$fleet$setFishPoin()
#' 
#' # Assign cell id to each fishing point
#' yourSmartRstudy$setCellPoin()
#' 
#' # Add week and month number to each point
#' yourSmartRstudy$fleet$setWeekMonthNum()
#' 
#' 
#' #####################
#' ## Fishing Grounds ##
#' #####################
#' 
#' # Setup available data to identify fishing areas
#' yourSmartRstudy$setAvailData()
#' 
#' # Setup cluster analysis input
#' yourSmartRstudy$sampMap$setClusInpu()
#' 
#' # Run cluster analysis with the SKATER method
#' yourSmartRstudy$sampMap$calcFishGrou(numCuts = 3, minsize = 10,
#'  modeska = "S", skater_method = "manhattan", nei_queen = FALSE)
#' 
#' # Setup cluster plot with 3 clusters
#' yourSmartRstudy$sampMap$setCutResult(ind_clu = 3)
#' 
#' # Map of the clusters' configuration
#' print(yourSmartRstudy$sampMap$ggCutFGmap)
#' 



SmartProject <- R6Class("smartProject",
                        portable = FALSE,
                        class = TRUE,
                        public = list(
                          rawDataSurvey = NULL,
                          yearInSurvey = NULL,
                          specieInSurvey = NULL,
                          surveyBySpecie = NULL,
                          rawDataFishery = NULL,
                          yearInFishery = NULL,
                          specieInFishery = NULL,
                          fisheryBySpecie = NULL,
                          ggEffRaw = NULL,
                          ggEffFish = NULL,
                          ggEffGrid = NULL,
                          ggEff = NULL,
                          gooLstCoho = list(),
                          sampMap = NULL,
                          fleet = NULL,
                          simProd = list(),
                          simEffo = NULL,
                          simBanFG = NULL,
                          simSpatialCost = NULL,
                          simEffortCost = NULL,
                          simProdCost = NULL,
                          simTotalCost = NULL,
                          simRevenue = list(),
                          simTotalRevenue = NULL,
                          simCostRevenue = NULL,
                          simResPlot = NULL,
                          outGmat = NULL,
                          outOptimEffo = NULL,
                          outWeiProp = NULL,
                          outWeiPropQ = NULL,
                          assessData = list(),
                          assSingleRes = list(),
                          assSinglePlot = list(),
                          assMultiRes = list(),
                          assMultiPlot = list(),
                          assessInteract = list(),
                          setAssessInteract = function(intName, intType, intWho, intQty, intChi, intOm) {
                            assessInteract <<- list()
                            assessInteract$name <<- intName
                            assessInteract$type <<- intType
                            assessInteract$who <<- intWho
                            assessInteract$qty <<- intQty
                            assessInteract$chi <<- intChi
                            assessInteract$om <<- intOm
                            # assessInteract$per <<- intPer
                          },
                          setAssessData = function(species, forecast = FALSE) {
                            cat("\nSetup Assessment Data for", species, "...\n")
                            if (is.null(assessData[[species]])) {
                              assessData[[species]] <<- list()
                            }
                            indSpeFis <- which(specieInFishery == species)
                            indSpeSur <- which(specieInSurvey == species)
                            
                            cat("\n\tLoading time-scales... ")
                            assessData[[species]]$Amax <<- fisheryBySpecie[[indSpeFis]]$nCoho
                            assessData[[species]]$Yr1 <<- as.numeric(as.character(min(years(fisheryBySpecie[[indSpeFis]]$rawLFD$Date))))
                            
                            assessData[[species]]$Yr2 <<- as.numeric(as.character(max(years(fisheryBySpecie[[indSpeFis]]$rawLFD$Date))))
                            if (forecast) assessData[[species]]$Yr2 <<- assessData[[species]]$Yr2 + 1
                            
                            assessData[[species]]$Nyear <<- assessData[[species]]$Yr2 - assessData[[species]]$Yr1 + 1
                            
                            assessData[[species]]$Nlen <<- 2
                            assessData[[species]]$NCAL <<- 0
                            assessData[[species]]$Nsurvey <<- 1
                            assessData[[species]]$NSAA <<- assessData[[species]]$Nyear
                            assessData[[species]]$NSAL <<- 0
                            cat("Done!")
                            
                            ## Catch
                            cat("\n\tLoading Catch Data... ")
                            tmpDF <- aggregate(
                              Production ~ Year, data.frame(
                                Year = as.character(fleet$effoAllLoa$Year),
                                Production = apply(fleet$predProd[[species]], 1, sum),
                                stringsAsFactors = FALSE
                              ),
                              sum
                            )
                            assessData[[species]]$Catch <<- tmpDF$Production
                            if (forecast) {
                              if (is.null(simProd[[species]])) stop("\nMissing Simulated Production!\n")
                              assessData[[species]]$Catch <<- c(assessData[[species]]$Catch, sum(simProd[[species]]))
                            }
                            
                            ## Catch at Age
                            for (sex in 1:length(names(fisheryBySpecie[[indSpeFis]]$groMixout))) {
                              if (sex == 1) {
                                tmpGro <- fisheryBySpecie[[indSpeFis]]$groMixout[[sex]][, c("FG", "Age", "Date")]
                              } else {
                                tmpGro <- rbind(tmpGro, fisheryBySpecie[[indSpeFis]]$groMixout[[sex]][, c("FG", "Age", "Date")])
                              }
                            }
                            tmpGro$Season <- factor(quarters(tmpGro$Date + 30), levels = c("1Q", "2Q", "3Q", "4Q"), labels = c("winter", "spring", "summer", "fall"))
                            tmpAstat <- ddply(tmpGro, .(FG, Age, Season), summarise, Freq = length(Age))
                            tmpAstat <- tmpAstat[tmpAstat$Freq > 0, ]
                            outWeiQ <- list()
                            for (season in c("winter", "spring", "summer", "fall")) {
                              preCAA <- data.frame(FG = 1:(sampMap$cutFG + 1))
                              preCAA <- cbind(preCAA, setNames(lapply(0:(fisheryBySpecie[[indSpeFis]]$nCoho - 1), function(x) x <- NA), 0:(fisheryBySpecie[[indSpeFis]]$nCoho - 1)))
                              for (i in 1:nrow(preCAA)) {
                                tempRev <- tmpAstat[tmpAstat$FG == i & tmpAstat$Season == season, ]
                                if (nrow(tempRev) > 0) {
                                  tempRev$propAge <- tempRev$Freq / sum(tempRev$Freq)
                                  outClass <- merge(data.frame(Age = 0:(fisheryBySpecie[[indSpeFis]]$nCoho - 1)), aggregate(formula = propAge ~ Age, data = tempRev, FUN = sum), all.x = TRUE)
                                  preCAA[i, 2:(fisheryBySpecie[[indSpeFis]]$nCoho + 1)] <- outClass$propAge
                                }
                              }
                              outWeiQ[[season]] <- preCAA
                            }
                            outProp <- matrix(
                              data = 0,
                              nrow = nrow(fleet$predProd[[species]]),
                              ncol = fisheryBySpecie[[indSpeFis]]$nCoho
                            )
                            tmpSeason <- data.frame(
                              Month = 1:12,
                              Season = c(
                                "winter", "winter", "spring",
                                "spring", "spring", "summer",
                                "summer", "summer", "fall",
                                "fall", "fall", "winter"
                              )
                            )
                            for (season in c("winter", "spring", "summer", "fall")) {
                              tmpOutProp <- apply(fleet$predProd[[species]][fleet$effoAllLoa$MonthNum %in% tmpSeason$Month[tmpSeason$Season == season], ], 1, function(x) apply(outWeiQ[[season]][, -1] * t(x), 2, sum, na.rm = TRUE))
                              outProp[fleet$effoAllLoa$MonthNum %in% tmpSeason$Month[tmpSeason$Season == season], ] <- t(tmpOutProp)
                            }
                            outCAA <- aggregate(. ~ Year, data.frame(Year = fleet$effoAllLoa$Year, outProp), sum)
                            
                            assessData[[species]]$YCAA <<- as.numeric(as.character(outCAA[, 1]))
                            if (forecast) assessData[[species]]$YCAA <<- c(assessData[[species]]$YCAA, assessData[[species]]$Yr2)
                            
                            assessData[[species]]$CAA <<- outCAA[, -1] / apply(outCAA[, -1], 1, sum)
                            if (forecast) {
                              if (is.null(simProd[[species]])) stop("\nMissing Simulated Production!\n")
                              if (is.null(simEffo)) stop("\nMissing Simulated Effort!\n")
                              simProp <- matrix(
                                data = 0,
                                nrow = nrow(simProd[[species]]),
                                ncol = fisheryBySpecie[[indSpeFis]]$nCoho
                              )
                              for (season in c("winter", "spring", "summer", "fall")) {
                                tmpOutProp <- apply(simProd[[species]][simEffo$MonthNum %in% tmpSeason$Month[tmpSeason$Season == season], ], 1, function(x) apply(outWeiQ[[season]][, -1] * t(x), 2, sum, na.rm = TRUE))
                                simProp[simEffo$MonthNum %in% tmpSeason$Month[tmpSeason$Season == season], ] <- t(tmpOutProp)
                              }
                              simCAA <- aggregate(. ~ Year, data.frame(Year = simEffo$Year, simProp), sum)
                              assessData[[species]]$CAA <<- rbind(assessData[[species]]$CAA, simCAA[, -1] / apply(simCAA[, -1], 1, sum))
                            }
                            
                            assessData[[species]]$NCAA <<- nrow(assessData[[species]]$CAA)
                            assessData[[species]]$YCAL <<- 0
                            assessData[[species]]$CAL <<- matrix(0, ncol = assessData[[species]]$Nlen, nrow = assessData[[species]]$NCAL)
                            cat("Done!")
                            
                            ## Survey at Age
                            cat("\n\tLoading Survey Data... ")
                            for (sex in 1:length(names(surveyBySpecie[[indSpeSur]]$groMixout))) {
                              if (sex == 1) {
                                tmpAL <- table(
                                  round(surveyBySpecie[[indSpeSur]]$groMixout[[sex]]$Length),
                                  factor(surveyBySpecie[[indSpeSur]]$groMixout[[sex]]$Age, levels = 0:surveyBySpecie[[indSpeSur]]$nCoho)
                                )
                                tmpAL <- round(tmpAL / apply(tmpAL, 1, sum), 2)
                                tmpAlDf <- as.data.frame.matrix(tmpAL)
                                tmpAlDf$Class <- rownames(tmpAL)
                                tmpAbu <- aggregate(. ~ Class + Year, surveyBySpecie[[indSpeSur]]$abuAvg[, c("Class", "Year", paste0("wei", substr(names(surveyBySpecie[[indSpeSur]]$groMixout)[sex], 1, 3)))], sum)
                                # tmpAbu <- aggregate(. ~ Class + Year, surveyBySpecie[[indSpeSur]]$abuAvg[,c("Class", "Year", names(surveyBySpecie[[indSpeSur]]$groMixout)[sex])], sum)
                                tmpMem <- merge(tmpAbu, tmpAlDf)
                                for (numCoh in 1:ncol(tmpAL)) {
                                  tmpMem[, 3 + numCoh] <- tmpMem[, 3 + numCoh] * tmpMem[, 3]
                                }
                                colnames(tmpMem)[3] <- "Num"
                                outMem <- tmpMem
                              } else {
                                tmpAL <- table(
                                  round(surveyBySpecie[[indSpeSur]]$groMixout[[sex]]$Length),
                                  factor(surveyBySpecie[[indSpeSur]]$groMixout[[sex]]$Age, levels = 0:surveyBySpecie[[indSpeSur]]$nCoho)
                                )
                                tmpAL <- round(tmpAL / apply(tmpAL, 1, sum), 2)
                                tmpAlDf <- as.data.frame.matrix(tmpAL)
                                tmpAlDf$Class <- rownames(tmpAL)
                                tmpAbu <- aggregate(. ~ Class + Year, surveyBySpecie[[indSpeSur]]$abuAvg[, c("Class", "Year", names(surveyBySpecie[[indSpeSur]]$groMixout)[sex])], sum)
                                tmpMem <- merge(tmpAbu, tmpAlDf)
                                for (numCoh in 1:ncol(tmpAL)) {
                                  tmpMem[, 3 + numCoh] <- tmpMem[, 3 + numCoh] * tmpMem[, 3]
                                }
                                colnames(tmpMem)[3] <- "Num"
                                outMem <- rbind(outMem, tmpMem)
                              }
                            }
                            outSAA <- aggregate(. ~ Year, outMem[, c(2, 4:(4 + (surveyBySpecie[[indSpeSur]]$nCoho - 1)))], sum)
                            
                            assessData[[species]]$SAA <<- outSAA[, -1]
                            if (forecast) assessData[[species]]$SAA <<- rbind(assessData[[species]]$SAA, assessData[[species]]$SAA[nrow(assessData[[species]]$SAA), ])
                            
                            assessData[[species]]$YSAA <<- as.numeric(as.character(outSAA[, 1]))
                            if (forecast) assessData[[species]]$YSAA <<- c(assessData[[species]]$YSAA, assessData[[species]]$Yr2)
                            
                            assessData[[species]]$SSAA <<- rep(1, nrow(assessData[[species]]$SAA))
                            assessData[[species]]$SAL <<- matrix(0, ncol = assessData[[species]]$Nlen, nrow = assessData[[species]]$NSAL)
                            assessData[[species]]$YSAL <<- 0
                            assessData[[species]]$SSAL <<- 0
                            assessData[[species]]$SelsurvType <<- 1
                            assessData[[species]]$SelexType <<- 1
                            cat("Done!")
                            
                            ### Mean Weight
                            cat("\n\tLoading Weight and Growth Data... ")
                            for (sex in 1:length(names(fisheryBySpecie[[indSpeFis]]$groMixout))) {
                              if (sex == 1) {
                                tmpWei <- fisheryBySpecie[[indSpeFis]]$groMixout[[sex]][, c("Age", "Weight")]
                              } else {
                                tmpWei <- rbind(tmpWei, fisheryBySpecie[[indSpeFis]]$groMixout[[sex]][, c("Age", "Weight")])
                              }
                            }
                            tmpWei$Age <- factor(tmpWei$Age, levels = 0:(assessData[[species]]$Amax - 1))
                            
                            assessData[[species]]$WeightS <<- (ddply(tmpWei, .(Age), summarise, coh.mean = min(Weight) / 1000, .drop = FALSE))[, 2]
                            weightNaN <- which(is.nan(assessData[[species]]$WeightS))
                            if (length(weightNaN) > 0) {
                              for (i in 1:length(weightNaN)) {
                                if (weightNaN[i] == 1) {
                                  if (!is.na(assessData[[species]]$WeightS[2])) {
                                    assessData[[species]]$WeightS[1] <<- mean(c(0, assessData[[species]]$WeightS[2]))
                                  } else {
                                    assessData[[species]]$WeightS[1] <<- mean(assessData[[species]]$WeightS[-1], na.rm = TRUE)
                                  }
                                  next
                                }
                                if (weightNaN[i] == length(assessData[[species]]$WeightS)) {
                                  if (!is.na(assessData[[species]]$WeightS[weightNaN[i] - 1])) {
                                    assessData[[species]]$WeightS[weightNaN[i]] <<- assessData[[species]]$WeightS[weightNaN[i] - 1]
                                  } else {
                                    assessData[[species]]$WeightS[weightNaN[i]] <<- mean(assessData[[species]]$WeightS[-(weightNaN[i])], na.rm = TRUE)
                                  }
                                  next
                                }
                                if (!is.na(assessData[[species]]$WeightS[weightNaN[i] - 1] & assessData[[species]]$WeightS[weightNaN[i] + 1])) {
                                  assessData[[species]]$WeightS[weightNaN[i]] <<- mean(c(assessData[[species]]$WeightS[weightNaN[i] - 1], assessData[[species]]$WeightS[weightNaN[i] + 1]))
                                } else {
                                  assessData[[species]]$WeightS[weightNaN[i]] <<- mean(assessData[[species]]$WeightS[-(weightNaN[i])], na.rm = TRUE)
                                }
                              }
                            }
                            
                            assessData[[species]]$WeightH <<- (ddply(tmpWei, .(Age), summarise, coh.mean = mean(Weight) / 1000, .drop = FALSE))[, 2]
                            weightNaN <- which(is.nan(assessData[[species]]$WeightH))
                            if (length(weightNaN) > 0) {
                              for (i in 1:length(weightNaN)) {
                                if (weightNaN[i] == 1) {
                                  if (!is.na(assessData[[species]]$WeightH[2])) {
                                    assessData[[species]]$WeightH[1] <<- mean(c(0, assessData[[species]]$WeightH[2]))
                                  } else {
                                    assessData[[species]]$WeightH[1] <<- mean(assessData[[species]]$WeightH[-1], na.rm = TRUE)
                                  }
                                  next
                                }
                                if (weightNaN[i] == length(assessData[[species]]$WeightH)) {
                                  if (!is.na(assessData[[species]]$WeightH[weightNaN[i] - 1])) {
                                    assessData[[species]]$WeightH[weightNaN[i]] <<- assessData[[species]]$WeightH[weightNaN[i] - 1]
                                  } else {
                                    assessData[[species]]$WeightH[weightNaN[i]] <<- mean(assessData[[species]]$WeightH[-(weightNaN[i])], na.rm = TRUE)
                                  }
                                  next
                                }
                                if (!is.na(assessData[[species]]$WeightH[weightNaN[i] - 1] & assessData[[species]]$WeightH[weightNaN[i] + 1])) {
                                  assessData[[species]]$WeightH[weightNaN[i]] <<- mean(c(assessData[[species]]$WeightH[weightNaN[i] - 1], assessData[[species]]$WeightH[weightNaN[i] + 1]))
                                } else {
                                  assessData[[species]]$WeightH[weightNaN[i]] <<- mean(assessData[[species]]$WeightH[-(weightNaN[i])], na.rm = TRUE)
                                }
                              }
                            }
                            
                            
                            assessData[[species]]$Qinit <<- 0
                            assessData[[species]]$InitN <<- rep(0, fisheryBySpecie[[indSpeFis]]$nCoho)
                            assessData[[species]]$InitR0 <<- 20
                            assessData[[species]]$recDev <<- rep(0, assessData[[species]]$Nyear)
                            
                            ### Growth
                            for (sex in 1:length(names(fisheryBySpecie[[indSpeFis]]$groPars))) {
                              if (sex == 1) {
                                tmpLHat <- fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$LHat
                                tmpkHat <- fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$kHat
                                tmpt0Hat <- fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$t0Hat
                              } else {
                                tmpLHat <- cbind(tmpLHat, fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$LHat)
                                tmpkHat <- cbind(tmpkHat, fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$kHat)
                                tmpt0Hat <- cbind(tmpt0Hat, fisheryBySpecie[[indSpeFis]]$groPars[[sex]]$t0Hat)
                              }
                            }
                            GrowthAL <- c(
                              mean(tmpLHat),
                              mean(tmpkHat),
                              -mean(tmpt0Hat),
                              0.15,
                              0.1
                            )
                            
                            for (sex in 1:length(names(fisheryBySpecie[[indSpeFis]]$LWpar))) {
                              if (sex == 1) {
                                tmpAlpha <- fisheryBySpecie[[indSpeFis]]$LWpar[[sex]]$alpha
                                tmpBeta <- fisheryBySpecie[[indSpeFis]]$LWpar[[sex]]$beta
                              } else {
                                tmpAlpha <- cbind(tmpAlpha, fisheryBySpecie[[indSpeFis]]$LWpar[[sex]]$alpha)
                                tmpBeta <- cbind(tmpBeta, fisheryBySpecie[[indSpeFis]]$LWpar[[sex]]$beta)
                              }
                            }
                            GrowthLW <- c(
                              mean(tmpAlpha),
                              mean(tmpBeta)
                            )
                            LenClassMax <- seq(from = assessData[[species]]$Nlen, by = assessData[[species]]$Nlen, length = assessData[[species]]$Nlen)
                            GetALK <- GetALKMW(GrowthAL[1], GrowthAL[2], GrowthAL[3], GrowthAL[4], GrowthAL[5], GrowthLW[1], GrowthLW[2], assessData[[species]]$Amax, LenClassMax, 0.0)
                            ALK1 <- GetALK$ALK
                            GetALK <- GetALKMW(GrowthAL[1], GrowthAL[2], GrowthAL[3], GrowthAL[4], GrowthAL[5], GrowthLW[1], GrowthLW[2], assessData[[species]]$Amax, LenClassMax, 0.5)
                            ALK2 <- GetALK$ALK
                            WeightLen <- GrowthLW[1] * GetALK$Lengths^GrowthLW[2] / 1000
                            assessData[[species]]$ALK1 <<- ALK1
                            assessData[[species]]$ALK2 <<- ALK2
                            assessData[[species]]$WeightLen <<- WeightLen
                            SurvBio <- matrix(0, nrow = 10, ncol = assessData[[species]]$Nyear)
                            SurvBio[1, ] <- apply(t(assessData[[species]]$SAA) * assessData[[species]]$WeightH, 2, sum)
                            assessData[[species]]$SurvBio <<- SurvBio
                            cat("Done!")
                            
                            # from GUI
                            cat("\n\tLoading User Parameters... ")
                            # assessData[[species]]$M <<- c(2.3, 1.1, 0.8, 0.7)
                            if (is.null(assessData[[species]]$M)) {
                              assessData[[species]]$M <<- round(seq(from = 0.7, to = 0.3, length.out = assessData[[species]]$Amax), 2)
                            }
                            # assessData[[species]]$Mat <<- c(0.5, 1, 1, 1)
                            if (is.null(assessData[[species]]$Mat)) {
                              assessData[[species]]$Mat <<- round(seq(from = 0.2, to = 0.9, length.out = assessData[[species]]$Amax), 2)
                            }
                            # assessData[[species]]$Selex <<- c(0.1, 0.2, 0.6, 1)
                            if (is.null(assessData[[species]][["Selex"]])) {
                              assessData[[species]]$Selex <<- round(seq(from = 0.2, to = 0.9, length.out = assessData[[species]]$Amax), 2)
                            }
                            if (is.null(assessData[[species]]$SelexSurv)) {
                              assessData[[species]]$SelexSurv <<- matrix(0, nrow = 10, ncol = assessData[[species]]$Amax)
                            }
                            assessData[[species]]$SelexSurv[1, ] <<- round(seq(from = 0.5, to = 0.9, length.out = assessData[[species]]$Amax), 2)
                            # assessData[[species]]$SelexSurv[1,] <<- c(0.2, 0.5, 1, 1)
                            if (length(assessData[[species]]$PropZBeforeMat) == 0) {
                              assessData[[species]]$PropZBeforeMat <<- 0.3
                            }
                            cat("Done!\n\nSetup Assessment Data for", species, "Completed!\n")
                          },
                          assSingle = function(species = "") {
                            if (is.null(assessData[[species]])) {
                              stop("Missing Data! Run setAssesData() first.")
                            }
                            cat("\n\nAssessing ", species, "\n", sep = "")
                            assSingleRes[[species]] <<- list()
                            Nyear <- assessData[[species]]$Nyear
                            Amax <- assessData[[species]]$Amax
                            Nsurvey <- assessData[[species]]$Nsurvey
                            LogR0 <- assessData[[species]]$InitR0
                            RecDev <- rep(0, assessData[[species]]$Nyear)
                            Qinit <- rep(assessData[[species]]$Qinit, Nsurvey)
                            VecS <- NULL
                            VecC <- NULL
                            if (assessData[[species]]$SelsurvType == 1) {
                              VecS <- rep(0, Nsurvey * (Amax - 2))
                              for (Isurv in 1:Nsurvey) {
                                Offset <- (Amax - 2) * (Isurv - 1)
                                for (Iage in 1:(Amax - 2)) {
                                  VecS[Offset + Iage] <- log((1 - assessData[[species]]$SelexSurv[Isurv, Iage]) / assessData[[species]]$SelexSurv[Isurv, Iage])
                                }
                              }
                            }
                            if (assessData[[species]]$SelexType == 1) {
                              VecC <- rep(0, Amax - 2)
                              for (Iage in 1:(Amax - 2)) {
                                VecC[Iage] <- log((1 - assessData[[species]]$Selex[Iage]) / assessData[[species]]$Selex[Iage])
                              }
                            }
                            Fvals <- rep(0, Nyear)
                            InitF <- log(1)
                            Pars <- c(assessData[[species]]$InitN, RecDev, LogR0, VecS, VecC, Fvals, InitF)
                            Npar <- length(Pars)
                            Res <- fit1Pars(Pars,
                                            fun1opt,
                                            FullMin = TRUE,
                                            DoVarCo = TRUE,
                                            SpeciesData = assessData[[species]]
                            )
                            assSingleRes[[species]] <<- fun1opt(Res$par,
                                                                DoEst = FALSE,
                                                                SpeciesData = assessData[[species]]
                            )
                            assSingleRes[[species]]$par <<- Res$par
                            assSingleRes[[species]]$VarCo <<- Res$VarCo
                            assSingleRes[[species]]$SSBSD <<- Res$SSBSD
                            cat("\n\n", species, " Assessment Complete!\n", sep = "")
                          },
                          assMulti = function(varCov = TRUE) {
                            if (is.null(assessData)) {
                              stop("Missing Data! Run setAssesData() first.")
                            }
                            cat("\n\nMultiple Species Assessment\n", sep = "")
                            assMultiRes <<- list()
                            Nspecies <- length(assessData)
                            Pars <- NULL
                            ParsInit <- NULL
                            for (Ispec in 1:Nspecies) {
                              Qinit <- rep(assessData[[Ispec]]$Qinit, assessData[[Ispec]]$Nsurvey)
                              Amax <- assessData[[Ispec]]$Amax
                              Nsurvey <- assessData[[Ispec]]$Nsurvey
                              Nyear <- assessData[[Ispec]]$Nyear
                              InitN <- assessData[[Ispec]]$InitN
                              RecDev <- assessData[[Ispec]]$recDev
                              LogR0 <- assessData[[Ispec]]$InitR0
                              
                              VecS <- NULL
                              VecC <- NULL
                              if (assessData[[Ispec]]$SelsurvType == 1) {
                                VecS <- rep(0, Nsurvey * (Amax - 2))
                                for (Isurv in 1:Nsurvey) {
                                  Offset <- (Amax - 2) * (Isurv - 1)
                                  for (Iage in 1:(Amax - 2)) VecS[Offset + Iage] <- log((1 - assessData[[Ispec]]$SelexSurv[Isurv, Iage]) / assessData[[Ispec]]$SelexSurv[Isurv, Iage])
                                }
                              }
                              if (assessData[[Ispec]]$SelsurvType == 2) {
                                for (Isurv in 1:Nsurvey) {
                                  Offset <- 3*(Isurv-1)
                                  VecS[(Offset + 1):(Offset + 3)] <- log(assessData[[Ispec]]$SelexSurv[Isurv, 1:3])
                                }
                              }
                              
                              
                              if (assessData[[Ispec]]$SelexType == 1) {
                                VecC <- rep(0, Amax - 2)
                                for (Iage in 1:(Amax - 2)) {
                                  VecC[Iage] <- log((1 - assessData[[Ispec]]$Selex[Iage]) / assessData[[Ispec]]$Selex[Iage])
                                }
                              }
                              if (assessData[[Ispec]]$SelexType == 2) {
                                VecC[1:3] <- log(assessData[[Ispec]]$Selex[1:3])
                              }
                              
                              Fvals <- rep(0, Nyear)
                              InitF <- log(1)
                              Pars <- c(Pars, InitN, RecDev, LogR0, VecS, VecC, Fvals, InitF)
                              
                              # if(toOpt[1] == TRUE) Pars <- c(Pars, InitN) else ParsInit <- c(ParsInit, InitN)
                              # if(toOpt[2] == TRUE) Pars <- c(Pars, RecDev) else ParsInit <- c(ParsInit, RecDev)
                              # if(toOpt[3] == TRUE) Pars <- c(Pars, LogR0) else ParsInit <- c(ParsInit, LogR0)
                              # if(toOpt[4] == TRUE) Pars <- c(Pars, VecS) else ParsInit <- c(ParsInit, VecS)
                              # if(toOpt[5] == TRUE) Pars <- c(Pars, VecC) else ParsInit <- c(ParsInit, VecC)
                              # if(toOpt[6] == TRUE) Pars <- c(Pars, Fvals) else ParsInit <- c(ParsInit, Fvals)
                              # if(toOpt[7] == TRUE) Pars <- c(Pars, InitF) else ParsInit <- c(ParsInit, InitF)
                              
                            }
                            Npar <- length(Pars)
                            Res <- fitNPars(Pars, funNopt,
                                            FullMin = TRUE, DoVarCo = varCov,
                                            SpeciesData = assessData, Nspecies = Nspecies,
                                            PredationPars = assessInteract
                            )
                            
                            assMultiRes <<- funNopt(Res$par,
                                                    DoEst = FALSE, SpeciesData = assessData,
                                                    Nspecies = Nspecies, PredationPars = assessInteract
                            )
                            assMultiRes$par <<- Res$par
                            assMultiRes$VarCo <<- Res$VarCo
                            assMultiRes$SSBSD <<- Res$SSBSD
                            cat("\n\nAssessment Complete!\n", sep = "")
                          },
                          setPlotSingle = function(species = "") {
                            if (is.null(assSingleRes[[species]])) {
                              stop("\n\nMissing Results! Run Assessment First.")
                            }
                            assSinglePlot[[species]] <<- list()
                            
                            ssbData <- data.frame(
                              Year = 1:assessData[[species]]$Nyear + assessData[[species]]$Yr1 - 1,
                              SSB = assSingleRes[[species]]$SSB,
                              Lower = assSingleRes[[species]]$SSB - 1.96 * assSingleRes[[species]]$SSBSD,
                              Upper = assSingleRes[[species]]$SSB + 1.96 * assSingleRes[[species]]$SSBSD
                            )
                            
                            assSinglePlot[[species]]$SSB <<- ggplot_SSBsingle(choSpecie = species, assData = ssbData)
                            
                            
                            survData <- list()
                            tmpObs <- assSingleRes[[species]]$ObsSAA
                            tmpObs$Year <- assSingleRes[[species]]$YSAA
                            survData$obsSAA <- melt(tmpObs,
                                                    id.vars = "Year",
                                                    variable.name = "Age", value.name = "Index"
                            )
                            survData$obsSAA$Lower <- survData$obsSAA$Index * exp(-1.96 * assSingleRes[[species]]$CvIndex)
                            survData$obsSAA$Upper <- survData$obsSAA$Index * exp(+1.96 * assSingleRes[[species]]$CvIndex)
                            levels(survData$obsSAA$Age) <- paste("Age", levels(survData$obsSAA$Age))
                            
                            tmpPred <- as.data.frame(assSingleRes[[species]]$PredSAA)
                            colnames(tmpPred) <- colnames(assSingleRes[[species]]$ObsSAA)
                            tmpPred$Year <- assSingleRes[[species]]$YSAA
                            survData$predSAA <- melt(tmpPred,
                                                     id.vars = "Year",
                                                     variable.name = "Age", value.name = "Index"
                            )
                            levels(survData$predSAA$Age) <- paste("Age", levels(survData$predSAA$Age))
                            
                            assSinglePlot[[species]]$ObsPredSurv <<- ggplot_OPSsingle(choSpecie = species, assData = survData)
                            
                            
                            caaData <- data.frame(
                              Age = 0:(assSingleRes[[species]]$Amax - 1),
                              obsCAA = apply(assSingleRes[[species]]$ObsCAA, 2, sum),
                              predCAA = apply(assSingleRes[[species]]$PredCAA, 2, sum)
                            )
                            
                            assSinglePlot[[species]]$ObsPredCAA <<- ggplot_OPCsingle(choSpecie = species, assData = caaData)
                            
                            
                            catcData <- data.frame(
                              Year = 1:assessData[[species]]$Nyear + assessData[[species]]$Yr1 - 1,
                              Catch = assSingleRes[[species]]$ObsCatch,
                              Lower = assSingleRes[[species]]$ObsCatch * exp(-1.96 * assSingleRes[[species]]$CvCatch),
                              Upper = assSingleRes[[species]]$ObsCatch * exp(+1.96 * assSingleRes[[species]]$CvCatch)
                            )
                            
                            assSinglePlot[[species]]$totCatc <<- ggplot_TCsingle(choSpecie = species, assData = catcData)
                          },
                          setPlotMulti = function() {
                            if (is.null(assMultiRes)) {
                              stop("\n\nMissing Results! Run Assessment First.")
                            }
                            
                            for (species in 1:length(assessData)) {
                              assMultiPlot[[names(assessData)[species]]] <<- list()
                              
                              ssbData <- data.frame(
                                Year = 1:assessData[[species]]$Nyear + assessData[[species]]$Yr1 - 1,
                                SSB = assMultiRes$SSB[species, ],
                                Lower = assMultiRes$SSB[species, ] - 1.96 * assMultiRes$SSBSD[seq(from = 1, to = length(assMultiRes$SSBSD), by = length(assessData)) + species - 1],
                                Upper = assMultiRes$SSB[species, ] + 1.96 * assMultiRes$SSBSD[seq(from = 1, to = length(assMultiRes$SSBSD), by = length(assessData)) + species - 1]
                              )
                              
                              assMultiPlot[[names(assessData)[species]]]$SSB <<- ggplot_SSBsingle(choSpecie = names(assessData)[species], assData = ssbData)
                              
                              
                              survData <- list()
                              tmpObs <- as.data.frame(assMultiRes$ObsSAA[species, 1:assessData[[species]]$Nyear, 1:assessData[[species]]$Amax])
                              tmpObs$Year <- assMultiRes$YSAA[species, 1:assessData[[species]]$Nyear]
                              survData$obsSAA <- melt(tmpObs,
                                                      id.vars = "Year",
                                                      variable.name = "Age", value.name = "Index"
                              )
                              survData$obsSAA$Lower <- survData$obsSAA$Index * exp(-1.96 * assMultiRes$CvIndex)
                              survData$obsSAA$Upper <- survData$obsSAA$Index * exp(+1.96 * assMultiRes$CvIndex)
                              levels(survData$obsSAA$Age) <- paste("Age", substr(levels(survData$obsSAA$Age), start = 2, stop = 2))
                              
                              tmpPred <- as.data.frame(assMultiRes$PredSAA[species, 1:assessData[[species]]$Nyear, 1:assessData[[species]]$Amax])
                              # colnames(tmpPred) <- colnames(assSingleRes[[species]]$ObsSAA)
                              tmpPred$Year <- assMultiRes$YSAA[species, 1:assessData[[species]]$Nyear]
                              survData$predSAA <- melt(tmpPred,
                                                       id.vars = "Year",
                                                       variable.name = "Age", value.name = "Index"
                              )
                              levels(survData$predSAA$Age) <- paste("Age", substr(levels(survData$predSAA$Age), start = 2, stop = 2))
                              
                              assMultiPlot[[names(assessData)[species]]]$ObsPredSurv <<- ggplot_OPSsingle(choSpecie = names(assessData)[species], assData = survData)
                              
                              
                              caaData <- data.frame(
                                Age = 0:(assMultiRes$Amax[species] - 1),
                                obsCAA = apply(assMultiRes$ObsCAA[species, 1:assessData[[species]]$Nyear, 1:assessData[[species]]$Amax], 2, sum),
                                predCAA = apply(assMultiRes$PredCAA[species, 1:assessData[[species]]$Nyear, 1:assessData[[species]]$Amax], 2, sum)
                              )
                              
                              assMultiPlot[[names(assessData)[species]]]$ObsPredCAA <<- ggplot_OPCsingle(choSpecie = names(assessData)[species], assData = caaData)
                              
                              catcData <- data.frame(
                                Year = 1:assessData[[species]]$Nyear + assessData[[species]]$Yr1 - 1,
                                Catch = assMultiRes$ObsCatch[species, ],
                                Lower = assMultiRes$ObsCatch[species, ] * exp(-1.96 * assMultiRes$CvCatch),
                                Upper = assMultiRes$ObsCatch[species, ] * exp(+1.96 * assMultiRes$CvCatch)
                              )
                              
                              assMultiPlot[[names(assessData)[species]]]$totCatc <<- ggplot_TCsingle(choSpecie = names(assessData)[species], assData = catcData)
                            }
                          },
                          setCostInput = function() {
                            if (is.null(fleet$effortIndex)) stop("Missing Effort Index")
                            if (is.null(fleet$daysAtSea)) stop("Missing Days at Sea Index")
                            if (is.null(fleet$effoAllLoa)) stop("Missing Production Index")
                            fleet$setInSpatial()
                            fleet$setInEffort()
                            setInProduction()
                          },
                          setInProduction = function() {
                            tmp_Prod <- data.frame(
                              Year = fleet$effoProdAllLoa$Year,
                              I_NCEE = fleet$effoProdAllLoa$I_NCEE,
                              MonthNum = fleet$effoProdAllLoa$MonthNum,
                              Production = apply(fleet$effoProdAllLoa[, (4 + sampMap$cutFG + 1):ncol(fleet$effoProdAllLoa)], 1, sum)
                            )
                            agg_Prod <- aggregate(Production ~ I_NCEE + Year, tmp_Prod, sum)
                            tmp_effoCost <- fleet$rawEconomy[, c("VessID", "Year", "ProductionCost")]
                            fleet$inProductionReg <<- merge(
                              x = tmp_effoCost, y = agg_Prod,
                              by.x = c("VessID", "Year"), by.y = c("I_NCEE", "Year")
                            )
                          },
                          setDaysAtSea = function() {
                            cat("\nProcessing year: ", sep = "")
                            for (year in names(fleet$rawEffort)) {
                              cat(year, "... ", sep = "")
                              tmp_vmsY <- fleet$rawEffort[[year]][, c("I_NCEE", "DATE", "MonthNum")]
                              tmp_vmsY$DATE <- floor(tmp_vmsY$DATE)
                              tmp_vmsY <- tmp_vmsY[!duplicated(tmp_vmsY), ]
                              out_tbl <- table(I_NCEE = tmp_vmsY$I_NCEE, Month = tmp_vmsY$MonthNum)
                              if (year == names(fleet$rawEffort)[1]) {
                                tmp_days <- data.frame(Year = year, out_tbl)
                              } else {
                                tmp_days <- rbind(
                                  tmp_days,
                                  data.frame(Year = year, out_tbl)
                                )
                              }
                            }
                            fleet$daysAtSea <<- suppressWarnings(merge(tmp_days, data.frame(
                              I_NCEE = as.numeric(substr(fleet$rawRegister$CFR, 4, nchar(fleet$rawRegister$CFR))),
                              Loa = fleet$rawRegister$Loa,
                              Kw = fleet$rawRegister$Power.Main
                            ), all.x = TRUE))
                            cat(" Completed!", sep = "")
                          },
                          setEffortIndex = function() {
                            cat("\n\nComputing Effort Index... ", sep = "")
                            if (is.null(sampMap$fgWeigDist)) stop("Missing Harbours Weighted Distance")
                            if (is.null(sampMap$cutFG)) stop("Missing Fishing Grounds")
                            if (is.null(fleet$effoAllLoa)) stop("Missing Effort-LOA data")
                            tmp_ei <- apply(data.frame(mapply(
                              `*`,
                              fleet$effoAllLoa[, 4:(sampMap$cutFG + 4)],
                              sampMap$fgWeigDist
                            )), 1, sum)
                            fleet$effortIndex <<- data.frame(fleet$effoAllLoa[, c(1:3, ncol(fleet$effoAllLoa))],
                                                             EffInd = tmp_ei
                            )
                            cat("Completed!", sep = "")
                          },
                          setProductionIndex = function() {
                            fleet$prodIndex <<- data.frame(
                              Year = fleet$effoProdAllLoa$Year,
                              I_NCEE = fleet$effoProdAllLoa$I_NCEE,
                              MonthNum = fleet$effoProdAllLoa$MonthNum,
                              Production = apply(fleet$effoProdAllLoa[, (4 + sampMap$cutFG + 1):ncol(fleet$effoProdAllLoa)], 1, sum)
                            )
                          },
                          getHarbFgDist = function() {
                            if (is.null(fleet$regHarbsUni)) stop("Missing Harbours Coordinates")
                            cat("\nRunnnig Fishing ground average distances routine... ", cat = "")
                            setRegHarbBox()
                            setFgWeigDist()
                            cat("\nFishing ground average distances routine completed!", cat = "")
                          },
                          setFgWeigDist = function() {
                            cat("\n\tSetting Fishing ground weighted distances... ", cat = "")
                            harb_fg_dist <- spDists(
                              y = as.matrix(sampMap$cutResShpCent[, 1:2]),
                              x = as.matrix(fleet$regHarbsBox[, 2:3]), longlat = TRUE
                            )
                            dimnames(harb_fg_dist) <- list(
                              as.character(fleet$regHarbsBox$Name),
                              paste("FG", sampMap$cutResShpCent$FG, sep = "")
                            )
                            harb_fg_dist <- data.frame(harb_fg_dist)
                            harb_wei_dist <- numeric(length = ncol(harb_fg_dist))
                            names(harb_wei_dist) <- names(harb_fg_dist)
                            for (i in 1:ncol(harb_fg_dist)) {
                              harb_wei_dist[i] <- weighted.mean(harb_fg_dist[, i], fleet$regHarbsBox$relFreq)
                            }
                            sampMap$fgWeigDist <<- harb_wei_dist[order(as.numeric(substr(names(harb_wei_dist), 3, nchar(names(harb_wei_dist)))))]
                            cat(" OK!", cat = "")
                          },
                          setRegHarbBox = function() {
                            cat("\n\tComputing Harbours-FishingGround distances...", cat = "")
                            tmp_dist <- gDistance(sampMap$gridShp,
                                                  SpatialPoints(fleet$regHarbsUni[, 2:3], proj4string = CRS(proj4string(sampMap$gridShp))),
                                                  byid = TRUE
                            )
                            fleet$regHarbsUni$shpDist <<- apply(tmp_dist, 1, min)
                            fleet$regHarbsBox <<- fleet$regHarbsUni[fleet$regHarbsUni$shpDist < 0.5, ]
                            harb_cur_box <- as.data.frame(table(fleet$vmsRegister[fleet$vmsRegister$Port.Name %in% fleet$regHarbsBox$Name, ]$Port.Name))
                            colnames(harb_cur_box) <- c("Name", "absFreq")
                            harb_cur_box$relFreq <- harb_cur_box$absFreq / sum(harb_cur_box$absFreq)
                            fleet$regHarbsBox <<- merge(fleet$regHarbsBox, harb_cur_box)
                            cat(" OK!", cat = "")
                          },
                          loadSurveyLFD = function(csv_path) {
                            cat("\nLoading survey data...", sep = "")
                            rawDataSurvey <<- read.csv(file = csv_path, stringsAsFactors = FALSE, header = TRUE)
                            surveyBySpecie <<- list()
                            
                            cat("\nSetting Years... ", sep = "")
                            setYearSurvey()
                            cat(" from ", min(yearInSurvey), " to ", max(yearInSurvey), "\nSetting Species... ", sep = "")
                            setSpecieSurvey()
                            cat(" found: ", paste(specieInSurvey, collapse = " - "), "\nSplitting Species...", sep = "")
                            splitSpecieSurvey()
                            
                            if (!is.null(sampMap$cutResShp)) {
                              addFg2Survey()
                              setSpreaSurvey()
                              setSpatSurvey()
                              sampMap$set_ggMapFgSurvey(rawDataSurvey)
                            }
                            
                            cat(" completed!", sep = "")
                          },
                          loadFisheryLFD = function(csv_path) {
                            cat("\nLoading fishery data...", sep = "")
                            rawDataFishery <<- read.csv(file = csv_path, stringsAsFactors = FALSE, header = TRUE)
                            if (is.null(rawDataFishery$Unsex)) rawDataFishery$Unsex <<- rawDataFishery$Female + rawDataFishery$Male
                            fisheryBySpecie <<- list()
                            cat("\nSetting Years... ", sep = "")
                            setYearFishery()
                            cat(" from ", min(levels(yearInFishery)[as.numeric(yearInFishery)]), " to ", max(levels(yearInFishery)[as.numeric(yearInFishery)]), "\nSetting Species... ", sep = "")
                            setSpecieFishery()
                            cat(" found: ", paste(specieInFishery, collapse = " - "), "\nSplitting Species...", sep = "")
                            splitSpecieFishery()
                            
                            if (!is.null(sampMap$cutResShp)) {
                              addFg2Fishery()
                              setSpreaFishery()
                              setSpatFishery()
                              sampMap$set_ggMapFgFishery(rawDataFishery)
                            }
                            
                            cat(" completed!", sep = "")
                          },
                          setYearSurvey = function() {
                            yearInSurvey <<- sort(unique(years(rawDataSurvey[, "Date"])), decreasing = FALSE)
                          },
                          setYearFishery = function() {
                            yearInFishery <<- sort(unique(years(rawDataFishery[, "Date"])), decreasing = FALSE)
                          },
                          loadMap = function(map_path) {
                            sampMap <<- SampleMap$new(map_path)
                          },
                          createFleet = function() {
                            fleet <<- FishFleet$new()
                          },
                          importEnv = function(envLst) {
                            sampMap <<- SampleMap$new()
                            sampMap$gridName <<- envLst$gridName
                            sampMap$gridShp <<- envLst$gridShp
                            sampMap$nCells <<- envLst$nCells
                            sampMap$gridPolySet <<- envLst$gridPolySet
                            sampMap$gridFortify <<- envLst$gridFortify
                            sampMap$griCent <<- envLst$griCent
                            sampMap$gridBbox <<- envLst$gridBbox
                            sampMap$gridBboxExt <<- envLst$gridBboxExt
                            sampMap$gridBboxSP <<- envLst$gridBboxSP
                            sampMap$gooMap <<- envLst$gooMap
                            sampMap$gooMapPlot <<- envLst$gooMapPlot
                            sampMap$plotRange <<- envLst$plotRange
                            sampMap$gooGrid <<- envLst$gooGrid
                            sampMap$gooBbox <<- envLst$gooBbox
                            sampMap$gridBathy <<- envLst$gridBathy
                            sampMap$centDept <<- envLst$centDept
                            sampMap$ggDepth <<- envLst$ggDepth
                            sampMap$bioDF <<- envLst$bioDF
                            sampMap$ggBioDF <<- envLst$ggBioDF
                          },
                          setSpecieSurvey = function() {
                            specieInSurvey <<- sort(unique(rawDataSurvey[, "Species"]))
                          },
                          setSpecieFishery = function() {
                            specieInFishery <<- sort(unique(rawDataFishery[, "Species"]))
                          },
                          splitSpecieSurvey = function() {
                            if (length(specieInSurvey) == 1) {
                              addSpecieSurvey(rawDataSurvey)
                            } else {
                              for (i in 1:length(specieInSurvey)) {
                                addSpecieSurvey(rawDataSurvey[rawDataSurvey[, "Species"] == specieInSurvey[i], ])
                              }
                            }
                          },
                          splitSpecieFishery = function() {
                            if (length(specieInFishery) == 1) {
                              addSpecieFishery(rawDataFishery)
                            } else {
                              for (i in 1:length(specieInFishery)) {
                                addSpecieFishery(rawDataFishery[rawDataFishery[, "Species"] == specieInFishery[i], ])
                              }
                            }
                          },
                          addSpecieSurvey = function(sing_spe) {
                            surveyBySpecie <<- c(surveyBySpecie, SurveyBySpecie$new(sing_spe))
                          },
                          addSpecieFishery = function(sing_spe) {
                            fisheryBySpecie <<- c(fisheryBySpecie, FisheryBySpecie$new(sing_spe))
                          },
                          setSpreaFishery = function() {
                            for (i in 1:length(fisheryBySpecie)) {
                              fisheryBySpecie[[i]]$setSpreDistSing()
                              fisheryBySpecie[[i]]$setAvailSex()
                            }
                          },
                          setSpatFishery = function() {
                            for (i in 1:length(fisheryBySpecie)) {
                              fisheryBySpecie[[i]]$setSpatDistSing()
                            }
                          },
                          setSpreaSurvey = function() {
                            for (i in 1:length(surveyBySpecie)) {
                              surveyBySpecie[[i]]$setSpreDistSing()
                              surveyBySpecie[[i]]$setAvailSex()
                            }
                          },
                          setSpatSurvey = function() {
                            for (i in 1:length(surveyBySpecie)) {
                              surveyBySpecie[[i]]$setSpatDistSing()
                            }
                          },
                          setDepthSurvey = function() {
                            cat("\n\nSetting depth of survey data:", sep = "")
                            for (i in 1:length(surveyBySpecie)) {
                              cat("\n\t", surveyBySpecie[[i]]$species, "... ", sep = "")
                              surveyBySpecie[[i]]$setDepth(bathyMatrix = sampMap$gridBathy)
                              cat("Done!", sep = "")
                            }
                            cat("\nCompleted!\n", sep = "")
                          },
                          setStratumSurvey = function(vectorStrata = c(0, 10, 100, 1000, Inf)) {
                            cat("\nSetting stratum of survey data:", sep = "")
                            for (i in 1:length(surveyBySpecie)) {
                              cat("\n\t", surveyBySpecie[[i]]$species, "... ", sep = "")
                              surveyBySpecie[[i]]$setStratum(vecStrata = vectorStrata)
                              cat("Done!", sep = "")
                            }
                            cat("\nCompleted!\n", sep = "")
                          },
                          setAbuAvgAll = function() {
                            cat("\nComputing average Number of individuals x Size x Stratum: ", sep = "")
                            for (i in 1:length(surveyBySpecie)) {
                              cat("\n\t", surveyBySpecie[[i]]$species, "... ", sep = "")
                              surveyBySpecie[[i]]$setAbuAvg()
                              cat("Done!", sep = "")
                            }
                            cat("\nCompleted!\n", sep = "")
                          },
                          setMeditsIndex = function() {
                            cat("\nComputing MEDITS index: ", sep = "")
                            for (i in 1:length(surveyBySpecie)) {
                              cat("\n\t", surveyBySpecie[[i]]$species, "... ", sep = "")
                              surveyBySpecie[[i]]$setIndSpe()
                              cat("Done!", sep = "")
                            }
                            cat("\nCompleted!\n", sep = "")
                          },
                          setStrataAbu = function() {
                            cat("\nComputing weighted Number of individuals x Size x Stratum: ", sep = "")
                            for (i in 1:length(surveyBySpecie)) {
                              cat("\n\t", surveyBySpecie[[i]]$species, "... ", sep = "")
                              surveyBySpecie[[i]]$abuAvg$weiFem <<- surveyBySpecie[[i]]$abuAvg$Female * sampMap$weightStrata[surveyBySpecie[[i]]$abuAvg$Stratum]
                              surveyBySpecie[[i]]$abuAvg$weiMal <<- surveyBySpecie[[i]]$abuAvg$Male * sampMap$weightStrata[surveyBySpecie[[i]]$abuAvg$Stratum]
                              surveyBySpecie[[i]]$abuAvg$weiUns <<- surveyBySpecie[[i]]$abuAvg$Unsex * sampMap$weightStrata[surveyBySpecie[[i]]$abuAvg$Stratum]
                              cat("Done!", sep = "")
                            }
                            cat("\nCompleted!\n", sep = "")
                          },
                          # setLFDPopSurvey = function(){
                          #   if(length(specieInSurvey) == 1){
                          #     calcLFDPopSurvey(1)
                          #   }else{
                          #     for(i in 1:length(specieInSurvey)){
                          #       calcLFDPopSurvey(i)
                          #     }}
                          #   # speDisPlot("All")
                          # },
                          # setLFDPopFishery = function(){
                          #   if(length(specieInFishery) == 1){
                          #     calcLFDPopFishery(1)
                          #   }else{
                          #     for(i in 1:length(specieInFishery)){
                          #       calcLFDPopFishery(i)
                          #     }}
                          #   # speDisPlot("All")
                          # },
                          loadFleeEffoDbs = function(effort_path, met_nam, onBox = TRUE, perOnBox = 1) {
                            cat("\n   ---   Extracting Effort data   ---", sep = "")
                            sort_files <- sort(effort_path)
                            fleet$rawEffort <<- list()
                            for (i in sort_files) {
                              cat("\n\nLoading db: ", i, sep = "")
                              # cat("\nSelecting tracks in box...", sep = "")
                              tmp_eff <- fn$sqldf("select * from (select * from (select *, rowid as i_id from intrp) join (select distinct I_NCEE, T_NUM from intrp where I_NCEE in (select distinct I_NCEE from nn_clas where met_des = '`met_nam`') and LON > `sampMap$gridBboxSP@bbox[1,1]` and LON < `sampMap$gridBboxSP@bbox[1,2]` and LAT > `sampMap$gridBboxSP@bbox[2,1]` and LAT < `sampMap$gridBboxSP@bbox[2,2]`) using (I_NCEE, T_NUM)) join (select * from p_depth) using (i_id)", dbname = i) ### Over in B-Box
                              if (onBox) {
                                in_box <- over(SpatialPoints(tmp_eff[, c("LON", "LAT")]), sampMap$gridBboxSP)
                              } else {
                                in_box <- over(
                                  SpatialPoints(tmp_eff[, c("LON", "LAT")], proj4string = CRS(proj4string(sampMap$gridShp))),
                                  unionSpatialPolygons(sampMap$gridShp,
                                                       IDs = rep(
                                                         1,
                                                         length(sampMap$gridShp@polygons)
                                                       )
                                  )
                                )
                              }
                              cat("   -   Completed!", sep = "")
                              
                              in_box[is.na(in_box)] <- 0
                              tmp_eff$in_box <- in_box
                              in_box_ping <- sqldf("select I_NCEE, T_NUM, sum(in_box) from tmp_eff group by I_NCEE, T_NUM")
                              all_ping <- sqldf("select I_NCEE, T_NUM, count(*) from tmp_eff group by I_NCEE, T_NUM")
                              
                              ### Selecting tracks with at least X points in the bounding box
                              if (perOnBox > 100) perOnBox <- 1
                              if (perOnBox > 1) perOnBox <- perOnBox / 100
                              cat("\nLoading tracks with at least ", perOnBox * 100, "% of points in the bounding box...", sep = "")
                              perOnInd <- in_box_ping[, 3] / all_ping[, 3] >= perOnBox
                              if (sum(perOnInd) == 0) {
                                cat("\nNo tracks with ", perOnBox * 100, "% of points in the bounding box.\nNo tracking data loaded!", sep = "")
                              } else {
                                all_in_box <- in_box_ping[perOnInd, 1:2]
                                all_sos <- sqldf("select * from tmp_eff join (select * from all_in_box) using (I_NCEE, T_NUM)")
                                cat("\nSaving Data", sep = "")
                                numYea <- as.character(sort(unique(years(all_sos$DATE))))
                                for (yea in 1:length(numYea)) {
                                  if (is.null(fleet$rawEffort[[numYea[yea]]])) {
                                    fleet$rawEffort[[numYea[yea]]] <<- all_sos[years(all_sos$DATE) == numYea[yea], ]
                                  } else {
                                    fleet$rawEffort[[numYea[yea]]] <<- rbind(fleet$rawEffort[[numYea[yea]]], all_sos[years(all_sos$DATE) == numYea[yea], ])
                                  }
                                }
                              }
                            }
                          },
                          # effPlot = function(whichYear){
                          #   if(whichYear == "All"){
                          #     all_sum <- apply(fleet$rawEffort, 1, sum)
                          #     round_perc <- 1+round(100*all_sum/max(all_sum))
                          #     # col_palette <- rev(heat.colors(100))
                          #     # cell_colors <- col_palette[round_perc]
                          #     distrPlotCols(cols = rev(heat.colors(101)), vals = round_perc,
                          #                   maxVal = ceiling(max(all_sum)),
                          #                   plotTitle = "Effort all years",
                          #                   legendUnits = "Hours")
                          #
                          #   }else{
                          #     num_col <- which(colnames(fleet$rawEffort) == whichYear)
                          #
                          #     yea_eff <- round(fleet$rawEffort[,num_col])
                          #     round_yea <- 1+100*yea_eff/max(yea_eff)
                          #
                          #     distrPlotCols(cols = rev(heat.colors(101)), vals = round_yea,
                          #                   maxVal = ceiling(max(yea_eff)),
                          #                   plotTitle = paste("Effort ", whichYear, sep = ""),
                          #                   legendUnits = "Hours")
                          #   }
                          # },
                          speDisPlot = function(whoPlo) {
                            if (whoPlo == "All") {
                              sampMap$plotSamMap("All species")
                              for (i in 1:length(specieInSurvey)) {
                                points(surveyBySpecie[[i]]$rawLFD[, c("Lon", "Lat")], pch = 20, col = 1 + i, cex = 0.4)
                              }
                            } else {
                              sampMap$plotSamMap(whoPlo)
                              points(surveyBySpecie[[which(specieInSurvey == whoPlo)]]$rawLFD[, c("Lon", "Lat")], pch = 20, col = 1 + which(specieInSurvey == whoPlo), cex = 0.4)
                            }
                          },
                          plotGooSpe = function(whiSpe, whiSou) {
                            if (whiSou == "Survey") {
                              if (whiSpe == "All") {
                                tmp_data <- unique(rawDataSurvey[, c("Species", "Lat", "Lon")])
                              } else {
                                tmp_data <- unique(rawDataSurvey[which(rawDataSurvey$Species == whiSpe), c("Species", "Lat", "Lon")])
                              }
                              # levels(tmp_data[,1]) <- unique(rawDataSurvey[,"Species"])
                              
                              sampMap$plotGooSpeSur(tmp_data)
                            } else {
                              if (whiSpe == "All") {
                                tmp_data <- unique(rawDataFishery[, c("Species", "Lat", "Lon")])
                              } else {
                                tmp_data <- unique(rawDataFishery[which(rawDataFishery$Species == whiSpe), c("Species", "Lat", "Lon")])
                              }
                              # levels(tmp_data[,1]) <- unique(rawDataFishery[,"Species"])
                              
                              sampMap$plotGooSpeFis(tmp_data)
                            }
                          },
                          setGooPlotCohoFish = function(species = "", sex = "Female", speCol = "", smooPoi = 500, smooBin = 0.5) {
                            cat("\n\nProcessing ", species, " - ", sex, "... Cohort ", sep = "")
                            gooLstCoho[[species]] <<- list()
                            gooLstCoho[[species]][[sex]] <<- list()
                            
                            tmpMix <- fisheryBySpecie[[which(specieInFishery == species)]]$groMixout[[sex]]
                            ageFGtbl <- table(tmpMix$FG, tmpMix$Age)
                            cohAbuFG <- as.data.frame(cbind(FG = as.numeric(rownames(ageFGtbl)), ageFGtbl))
                            outPalette <- rainbow(ncol(cohAbuFG) - 1)
                            
                            for (coh_i in 2:ncol(cohAbuFG)) {
                              cat(colnames(ageFGtbl)[coh_i - 1], "... ", sep = "")
                              
                              if (speCol == "") {
                                cohFillPal <- outPalette[coh_i - 1]
                              } else {
                                cohFillPal <- speCol
                              }
                              all_cell <- merge(x = sampMap$cutResShpFort$id, data.frame(x = cohAbuFG$FG, y = round(100 * cohAbuFG[, coh_i] / max(cohAbuFG[, coh_i]))), all = TRUE)
                              grid_data <- cbind(sampMap$cutResShpFort, NumInd = all_cell[, 2])
                              
                              
                              if (length(sampMap$cutFG) > 0) {
                                if (length(sampMap$gridShp@polygons) == (sampMap$cutFG + 1)) {
                                  tmp_coo <- data.frame(coordinates(sampMap$gridShp), cell_id = 1:length(sampMap$gridShp))
                                  colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                                } else {
                                  tmp_coo <- sampMap$cutResShpCent
                                }
                              } else {
                                tmp_coo <- sampMap$cutResShpCent
                              }
                              
                              tmp_dens <- data.frame(
                                lon = rep(grid_data[!is.na(grid_data$NumInd), ]$long, grid_data[!is.na(grid_data$NumInd), ]$NumInd),
                                lat = rep(grid_data[!is.na(grid_data$NumInd), ]$lat, grid_data[!is.na(grid_data$NumInd), ]$NumInd)
                              )
                              
                              mapCoho <- suppressMessages(sampMap$gooMapPlot +
                                                            geom_polygon(aes(x = long, y = lat, group = group, fill = NumInd),
                                                                         colour = "black", size = 0.1, data = grid_data, alpha = 0.8
                                                            ) +
                                                            scale_fill_gradient(paste0("Max abundance\n", max(cohAbuFG[, coh_i]), " specimens"),
                                                                                low = "Grey85", high = cohFillPal,
                                                                                # trans = "log10",
                                                                                breaks = pretty(1:100, 5), limits = c(1, 100)
                                                            ) +
                                                            ggtitle(paste0("Spatial Distribution of ", species, " - Cohort ", coh_i - 2)) +
                                                            geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                      data = tmp_coo, size = 2
                                                            ) +
                                                            theme(
                                                              legend.position = c(0.1, 0.22),
                                                              legend.text = element_text(size = 10, colour = "grey19"),
                                                              # legend.title = element_blank(),
                                                              legend.background = element_rect(fill = rgb(1, 1, 1, 0.5)),
                                                              axis.text.x = element_text(size = 10),
                                                              axis.title.x = element_text(size = 12),
                                                              axis.text.y = element_text(size = 10),
                                                              legend.key.size = unit(0.75, "cm"),
                                                              axis.title.y = element_text(size = 12),
                                                              plot.title = element_text(size = 20)
                                                            ) +
                                                            stat_density_2d(
                                                              data = tmp_dens,
                                                              mapping = aes(x = lon, y = lat, alpha = ..level..),
                                                              fill = cohFillPal,
                                                              geom = "polygon",
                                                              n = smooPoi,
                                                              h = smooBin,
                                                              show.legend = FALSE
                                                            ) +
                                                            scale_alpha_continuous(limits = c(0, 0.4), breaks = seq(0, 0.4, by = 0.025)) +
                                                            geom_polygon(aes(x = long, y = lat, group = group, alpha = 0.1),
                                                                         colour = "black", size = 0.1, data = grid_data, alpha = 0.8, fill = NA
                                                            ) +
                                                            geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                      data = tmp_coo, size = 2
                                                            ))
                              
                              gooLstCoho[[species]][[sex]][[colnames(ageFGtbl)[coh_i - 1]]] <<- mapCoho
                            }
                            cat(" Completed!", sep = "")
                          },
                          distrPlotCols = function(cols = NULL, vals = NULL, maxVal = 100,
                                                   plotTitle = "NoTitle", legendUnits = "NoUnits") {
                            def.par <- par(no.readonly = TRUE)
                            par(mar = c(2.5, 2.5, 3, 1))
                            layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(6, 1))
                            sampMap$plotSamMap(title = plotTitle, celCol = cols[vals])
                            par(mar = c(5, 3, 5, 1))
                            plot(NULL, xlim = c(0, 1), ylim = c(0, 1), bty = "n", axes = FALSE, ann = FALSE, main = "Hours")
                            rect(0.25, seq(0.2, 0.79, length.out = 100),
                                 0.55, seq(0.21, 0.80, length.out = 100),
                                 col = cols, border = rainbow(1, alpha = 0.01)
                            )
                            mtext(ceiling(seq(from = 0, to = maxVal, length.out = 10)),
                                  side = 2,
                                  at = seq(0.21, 0.80, length.out = 10), las = 2, cex = 1
                            )
                            text(legendUnits, x = 0.25, y = 0.15)
                            par(def.par)
                          },
                          ggplotRawPoints = function(year) {
                            tmp_dat <- fleet$rawEffort[[year]][
                              sample(1:nrow(fleet$rawEffort[[year]]), min(c(50000, nrow(fleet$rawEffort[[year]])))),
                              c("LON", "LAT", "W_HARB")
                              ]
                            tmp_dat$Status <- factor(tmp_dat$W_HARB,
                                                     levels = c("0", "1"),
                                                     labels = c("At sea", "In harbour")
                            )
                            ggEffRaw <<- suppressMessages(sampMap$gooMapPlot +
                                                            geom_point(
                                                              data = tmp_dat,
                                                              aes(x = LON, y = LAT, shape = Status, color = Status),
                                                              size = 0.6, alpha = 0.3
                                                            ) +
                                                            geom_point(
                                                              data = subset(tmp_dat, Status == "In harbour"),
                                                              aes(x = LON, y = LAT, shape = Status, color = Status),
                                                              size = 0.6, alpha = 0.3
                                                            ) +
                                                            scale_colour_manual(values = c("coral", "darkseagreen1")) +
                                                            guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
                                                            ggtitle(paste("Sample raw points - ", year, sep = "")) +
                                                            theme_tufte(base_size = 14, ticks = T) +
                                                            theme(
                                                              legend.position = "bottom",
                                                              axis.text.x = element_text(size = 8),
                                                              axis.title.x = element_text(size = 10),
                                                              panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                                              axis.text.y = element_text(size = 8),
                                                              axis.title.y = element_text(size = 10),
                                                              legend.text = element_text(size = 8),
                                                              legend.title = element_text(size = 10)
                                                            ))
                          },
                          setGgEff = function() {
                            if (is.null(ggEffRaw)) {
                              tmp_raw <- ggplot() + geom_blank() + ggtitle("Raw Effort Not Loaded Yet")
                            } else {
                              tmp_raw <- ggEffRaw
                            }
                            if (is.null(ggEffFish)) {
                              tmp_fish <- ggplot() + geom_blank() + ggtitle("Fishing Point Not Loaded Yet")
                            } else {
                              tmp_fish <- ggEffFish
                            }
                            if (is.null(ggEffGrid)) {
                              tmp_grid <- ggplot() + geom_blank() + ggtitle("Gridded Effort Not Loaded Yet")
                            } else {
                              tmp_grid <- ggEffGrid
                            }
                            ggEff <<- suppressWarnings(grid.arrange(tmp_raw,
                                                                    tmp_fish,
                                                                    tmp_grid,
                                                                    layout_matrix = matrix(1:3, 1, 3)
                            ))
                          },
                          plotGgEff = function() {
                            suppressWarnings(grid.draw(ggEff))
                          },
                          ggplotFgWeigDists = function() {
                            all_cell <- merge(
                              x = sampMap$cutResShpFort$id,
                              data.frame(
                                x = as.numeric(substr(
                                  names(sampMap$fgWeigDist), 3,
                                  nchar(names(sampMap$fgWeigDist))
                                )),
                                y = sampMap$fgWeigDist
                              ), all = TRUE
                            )
                            all_cell[is.na(all_cell)] <- 0
                            grid_data <- cbind(sampMap$cutResShpFort, DistAvg = all_cell[, 2])
                            
                            if (length(sampMap$cutFG) > 0) {
                              if (length(sampMap$gridShp@polygons) == (sampMap$cutFG + 1)) {
                                tmp_coo <- data.frame(coordinates(sampMap$gridShp), cell_id = 1:length(sampMap$gridShp))
                                colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                              } else {
                                tmp_coo <- sampMap$cutResShpCent
                              }
                            } else {
                              tmp_coo <- sampMap$cutResShpCent
                            }
                            
                            suppressWarnings(
                              print(
                                suppressMessages(
                                  sampMap$gooMapPlot +
                                    geom_polygon(aes(
                                      x = long, y = lat,
                                      group = group, fill = DistAvg
                                    ),
                                    colour = "black", size = 0.1,
                                    data = grid_data, alpha = 0.8
                                    ) +
                                    scale_fill_gradient("Weighted\nDistance", low = "snow1", high = "orange1") +
                                    geom_text(aes(label = FG, x = Lon, y = Lat),
                                              data = tmp_coo, size = 2
                                    ) +
                                    ggtitle("Average Distance x Fishing Ground") +
                                    xlab("Longitude") + ylab("Latitude") +
                                    geom_point(
                                      data = fleet$regHarbsBox,
                                      mapping = aes(x = Lon, y = Lat, size = absFreq),
                                      fill = NA, color = "tomato3", shape = 21
                                    ) +
                                    scale_size_continuous("Number of\nVessels",
                                                          breaks = pretty(fleet$regHarbsBox$absFreq, 5),
                                                          range = c(1, 15)
                                    ) +
                                    geom_label_repel(
                                      data = fleet$regHarbsBox, mapping = aes(x = Lon, y = Lat, label = Name),
                                      size = 3, nudge_x = 0.1, nudge_y = 0.1, color = "grey3", fill = "grey89"
                                    ) +
                                    theme_tufte(base_size = 14, ticks = T) +
                                    theme(
                                      legend.position = "bottom",
                                      axis.text.x = element_text(size = 8),
                                      axis.title.x = element_text(size = 10),
                                      panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                      axis.text.y = element_text(size = 8),
                                      axis.title.y = element_text(size = 10)
                                    )
                                )
                              )
                            )
                          },
                          setAvailData = function() {
                            sampMap$availData <<- character(0)
                            sampMap$rawInpu <<- list()
                            
                            cat("\n\nLoading available data:\n")
                            if (is.null(sampMap$bioDF)) {
                              stop("\nMissing Seabed Data!\n")
                            } else {
                              sampMap$availData <<- c(sampMap$availData, "Seabed")
                              cat("\n   -   Seabed Category")
                              sampMap$rawInpu <<- c(sampMap$rawInpu, Seabed = list(data.frame(sampMap$bioDF)))
                              cat("\t\t-   Loaded!")
                            }
                            
                            if (is.null(fleet$rawEffort)) {
                              stop("\nMissing Effort Data!\n")
                            } else {
                              sampMap$availData <<- c(sampMap$availData, "Effort")
                              cat("\n   -   Effort Distribution")
                              raw_effort <- numeric(length = sampMap$nCells)
                              for (i in names(fleet$rawEffort)) {
                                tmp_effo <- as.data.frame(table(fleet$rawEffort[[i]]$Cell[which(fleet$rawEffort[[i]]$FishPoint)]))
                                names(tmp_effo) <- c("Cell", "Freq")
                                tmp_effo$Cell <- as.numeric(as.character(tmp_effo$Cell))
                                miss_rows <- as.numeric(setdiff(as.character(sampMap$gridShp@plotOrder), as.character(tmp_effo$Cell)))
                                if (length(miss_rows) > 0) {
                                  tmp_effo <- rbind(tmp_effo, data.frame(Cell = miss_rows, Freq = 0))
                                  tmp_effo <- tmp_effo[order(tmp_effo[, 1]), ]
                                }
                                raw_effort <- cbind(raw_effort, tmp_effo[, 2])
                                colnames(raw_effort)[ncol(raw_effort)] <- paste("Year_", i, sep = "")
                              }
                              sampMap$rawInpu <<- c(sampMap$rawInpu, Effort = list(raw_effort[, -1]))
                              cat("\t-   Loaded!")
                            }
                            
                            if (is.null(sampMap$griCent)) {
                              stop("\nMissing Depth Data!\n")
                            } else {
                              sampMap$availData <<- c(sampMap$availData, "Depth")
                              cat("\n   -   Cell Depth")
                              sampMap$rawInpu <<- c(sampMap$rawInpu, Depth = list(sampMap$centDept[, 3]))
                              cat("\t\t-   Loaded!\n")
                            }
                          },
                          predictProduction = function(species) {
                            Prod <- matrix(data = NA, nrow(fleet$effoAllLoa), ncol = sampMap$cutFG + 1)
                            lyears <- sort(as.numeric(as.character(unique(fleet$effoAllLoa$Year))))
                            thrZero <- mean(fleet$effoProdAllLoa[, species][fleet$effoProdAllLoa[, species] < fleet$specSett[[species]]$threshold & fleet$effoProdAllLoa[, species] > 0])
                            fgClms <- which(colnames(fleet$effoAllLoa) %in% as.character(seq(1, sampMap$cutFG + 1)))
                            datalog <- fleet$effoAllLoa
                            datalog$MonthNum <- as.factor(datalog$MonthNum)
                            datalog$Year <- as.factor(datalog$Year)
                            if (fleet$specLogit[[species]]$logit$Name == "GLM") {
                              infish <- which(predict(fleet$specLogit[[species]]$logit$Model, datalog, type = "response") > fleet$specLogit[[species]]$logit$Cut)
                            } else {
                              infish <- which(predict(fleet$specLogit[[species]]$logit$Model, datalog, type = "prob")[, 2] > fleet$specLogit[[species]]$logit$Cut)
                            }
                            # infish <- which(predict(fleet$specLogit[[species]]$logit$logit_f, datalog, type="response") > fleet$specLogit[[species]]$optCut)
                            for (i in 1:length(infish)) {
                              idata <- as.numeric(fleet$effoAllLoa[infish[i], fgClms])
                              iloa <- as.numeric(fleet$effoAllLoa[infish[i], "Loa"])
                              iy <- which(lyears == fleet$effoAllLoa[infish[i], "Year"])
                              im <- as.numeric(as.character(fleet$effoAllLoa[infish[i], "MonthNum"]))
                              ib <- fleet$resNNLS[[species]]$bmat[which((fleet$resNNLS[[species]]$SceMat$YEAR == iy) & (fleet$resNNLS[[species]]$SceMat$MONTH == im)), ]
                              # Prod[infish[i]] <- sum(ib * idata * iloa) + mean(fleet$effoProdAllLoa[,species][fleet$effoProdAllLoa[,species] < fleet$specSett[[species]]$threshold & fleet$effoProdAllLoa[,species] > 0])
                              if (sum(ib * idata) > 0) {
                                Prod[infish[i], ] <- (ib * idata * iloa) + ((ib * idata) / sum(ib * idata)) * thrZero
                              }
                            }
                            Prod[is.na(Prod)] <- 0
                            colnames(Prod) <- paste("PR_", as.character(seq(1, ncol(Prod))), sep = "")
                            fleet$predProd[[species]] <<- Prod
                          },
                          simProdAll = function(selRow = numeric(0)) {
                            if (length(selRow) == 0) {
                              selRow <- 1:nrow(simEffo)
                            }
                            
                            lyears <- sort(as.numeric(as.character(unique(fleet$effoAllLoa$Year))))
                            datalog <- simEffo[selRow, ]
                            datalog$MonthNum <- as.factor(datalog$MonthNum)
                            datalog$Year <- as.factor(datalog$Year)
                            fgClms <- which(colnames(simEffo) %in% as.character(seq(1, sampMap$cutFG + 1)))
                            
                            for (species in names(fleet$specLogit)) {
                              Prod <- matrix(data = NA, length(selRow), ncol = sampMap$cutFG + 1)
                              thrZero <- mean(fleet$effoProdAllLoa[, species][fleet$effoProdAllLoa[, species] < fleet$specSett[[species]]$threshold & fleet$effoProdAllLoa[, species] > 0])
                              if (fleet$specLogit[[species]]$logit$Name == "GLM") {
                                infish <- which(predict(fleet$specLogit[[species]]$logit$Model, datalog, type = "response") > fleet$specLogit[[species]]$logit$Cut)
                              } else {
                                infish <- which(predict(fleet$specLogit[[species]]$logit$Model, datalog, type = "prob")[, 2] > fleet$specLogit[[species]]$logit$Cut)
                              }
                              
                              Prod[infish, ] <- do.call(rbind, lapply(split(simEffo[selRow[infish],], seq(nrow(simEffo[selRow[infish],]))),
                                                                      FUN = function(x) getSingleProd(oneEff = as.numeric(x),
                                                                                                      betas = fleet$resNNLS[[species]]$bmat,
                                                                                                      scenarios = fleet$resNNLS[[species]]$SceMat,
                                                                                                      thres = thrZero,
                                                                                                      fgs = fgClms)))
                              
                              Prod[is.na(Prod)] <- 0
                              colnames(Prod) <- paste("PR_", as.character(seq(1, ncol(Prod))), sep = "")
                              if (length(selRow) == nrow(simEffo)) {
                                simProd[[species]] <<- Prod
                              } else {
                                simProd[[species]][selRow, ] <<- Prod
                              }
                            }
                          },
                          genSimEffo = function(selRow = numeric(0), areaBan = numeric(0)) {
                            if (is.null(simEffo)) {
                              simEffo <<- fleet$effoAllLoa[fleet$effoAllLoa$Year == max(as.numeric(as.character(unique(fleet$effoAllLoa$Year)))), ]
                            } else {
                              
                              if (length(selRow) == 0) {
                                selRow <- 1:nrow(simEffo)
                              }
                              
                              obsZero <- apply(simEffo[selRow, 4:(ncol(simEffo) - 1)], 2, sum)
                              posZero <- which(obsZero == 0)
                              lenPosZero <- length(posZero)
                              lenAreaBan <- length(areaBan)
                              
                              if (lenPosZero == 0 & lenAreaBan == 0) {
                                selMode <- "flat"
                              } else if (lenPosZero == 0 & lenAreaBan > 0) {
                                selMode <- "ban"
                              } else if (lenPosZero > 0 & lenAreaBan == 0) {
                                selMode <- "zero"
                                areaBan <- posZero
                              } else if (lenPosZero > 0 & lenAreaBan > 0) {
                                selMode <- "zeroBan"
                                areaBan <- union(areaBan, posZero)
                              }
                              
                              simEffo[selRow, 4:(ncol(simEffo) - 1)] <<- switch(selMode,
                                                                                flat = {
                                                                                  t(apply(simEffo[selRow, 4:(ncol(simEffo) - 1)], 1, function(x) genFlatEffo(effoPatt = x)))
                                                                                },
                                                                                zero = {
                                                                                  t(apply(simEffo[selRow, 4:(ncol(simEffo) - 1)], 1, function(x) genBanEffo(effoPatt = x, set0 = areaBan)))
                                                                                },
                                                                                zeroBan = {
                                                                                  t(apply(simEffo[selRow, 4:(ncol(simEffo) - 1)], 1, function(x) genBanEffo(effoPatt = x, set0 = areaBan)))
                                                                                },
                                                                                ban = {
                                                                                  t(apply(simEffo[selRow, 4:(ncol(simEffo) - 1)], 1, function(x) genBanEffo(effoPatt = x, set0 = areaBan)))
                                                                                }
                              )
                            }
                          },
                          getSimSpatialCost = function() {
                            tmp_ei <- apply(data.frame(mapply(`*`, simEffo[, 4:(ncol(simEffo) - 1)], sampMap$fgWeigDist)), 1, sum)
                            tmpIndex <- data.frame(simEffo[, c(1:3, ncol(simEffo))], EffInd = tmp_ei)
                            simSpatialIndex <- aggregate(EffInd ~ I_NCEE + Year + Loa, tmpIndex, sum)
                            predSpatCost <- predict(fleet$outSpatialReg, simSpatialIndex)
                            simSpatialCost <<- cbind(simSpatialIndex, predSpatCost)
                          },
                          getSimEffortCost = function() {
                            effortIndex <- aggregate(Freq ~ I_NCEE + Year + Loa + Kw, fleet$daysAtSea[fleet$daysAtSea$Year == max(as.numeric(as.character(fleet$daysAtSea$Year))), ], sum)
                            predEffoCost <- predict(fleet$outEffortReg, effortIndex)
                            simEffortCost <<- cbind(effortIndex, predEffoCost)
                          },
                          getSimProdCost = function() {
                            outProd <- numeric(nrow(simProd[[1]]))
                            for (species in names(simProd)) {
                              outProd <- outProd + apply(simProd[[species]], 1, sum)
                            }
                            tmp_Prod <- data.frame(
                              Year = simEffo$Year,
                              I_NCEE = simEffo$I_NCEE,
                              MonthNum = simEffo$MonthNum,
                              Production = outProd
                            )
                            agg_ProdInd <- aggregate(Production ~ I_NCEE + Year, tmp_Prod, sum)
                            predProdCost <- predict(fleet$outProductionReg, agg_ProdInd)
                            simProdCost <<- cbind(agg_ProdInd, predProdCost = predProdCost)
                          },
                          getSimTotalCost = function() {
                            getSimSpatialCost()
                            getSimEffortCost()
                            getSimProdCost()
                            
                            tmp_out_costs <- merge(merge(simEffortCost, simSpatialCost), simProdCost)
                            out_costs <- tmp_out_costs[, c("I_NCEE", "Year", "Loa", "predEffoCost", "predSpatCost", "predProdCost")]
                            out_costs$totCost <- out_costs$predEffoCost + out_costs$predSpatCost + out_costs$predProdCost
                            
                            simTotalCost <<- out_costs
                          },
                          getSimRevenue = function(selRow = numeric(0), timeScale = "Year") {
                            if (length(selRow) == 0) {
                              selRow <- 1:nrow(simEffo)
                            }
                            
                            # simTotalRevenue <- data.frame(matrix(NA, nrow = nrow(simEffo), ncol = length(names(simProd))))
                            speNam <- names(simProd)
                            tmp_Revenue <- cbind(simEffo[, 1:3], (matrix(NA, nrow = nrow(simEffo[, ]), ncol = length(speNam))))
                            names(tmp_Revenue)[4:(4 + length(speNam) - 1)] <- speNam
                            for (species in speNam) {
                              if (timeScale == "Year") {
                                tmpRev <- getFleetRevenue(
                                  predProd = simProd[[species]][selRow, ],
                                  lwStat = outWeiProp[[species]][, -1],
                                  priceVec = fleet$ecoPrice[[species]]$Price
                                )
                              } else {
                                tmpRev <- getFleetRevSeason(
                                  predProd = simProd[[species]][selRow, ],
                                  monthVec = tmp_Revenue$MonthNum[selRow],
                                  lwStat = outWeiPropQ[[species]],
                                  priceVec = fleet$ecoPrice[[species]]$Price
                                )
                              }
                              if (length(selRow) == nrow(simEffo)) {
                                simRevenue[[species]] <<- tmpRev
                              } else {
                                simRevenue[[species]][selRow, ] <<- tmpRev
                              }
                              tmp_Revenue[, species] <- apply(simRevenue[[species]], 1, sum, na.rm = TRUE)
                            }
                            if (ncol(tmp_Revenue) == 4) {
                              tmp_Revenue$totRevenue <- tmp_Revenue[, 4]
                            } else {
                              tmp_Revenue$totRevenue <- apply(tmp_Revenue[, 4:(4 + length(speNam) - 1)], 1, sum)
                            }
                            simTotalRevenue <<- aggregate(totRevenue ~ I_NCEE + Year, tmp_Revenue, sum)
                          },
                          getLWstat = function(fill = FALSE) {
                            if (is.null(fisheryBySpecie)) {
                              stop("No mcmc output found")
                            }
                            specList <- intersect(
                              specieInFishery,
                              names(fleet$ecoPrice)
                            )
                            
                            for (species in specList) {
                              priIdx <- which(names(fleet$ecoPrice) == species)
                              fisIdx <- which(specieInFishery == species)
                              if (is.null(fleet$ecoPrice[[priIdx]])) {
                                stop(paste0("No size/price data found for species ", names(fleet$ecoPrice)[priIdx]))
                              }
                              
                              vecSize <- sort(unique(c(fleet$ecoPrice[[priIdx]]$LowerBound, fleet$ecoPrice[[priIdx]]$UpperBound)))
                              curUnit <- unique(fleet$ecoPrice[[priIdx]]$Units)[1]
                              
                              fisheryBySpecie[[fisIdx]]$setLWstat(lwUnit = curUnit)
                              fgNames <- paste0("LW_", 1:(sampMap$cutFG + 1))
                              
                              ## yearly
                              preRevenue <- data.frame(FG = 1:length(fgNames))
                              preRevenue <- cbind(preRevenue, setNames(lapply(1:(length(vecSize) - 1), function(x) x <- NA), 1:(length(vecSize) - 1)))
                              for (i in 1:nrow(preRevenue)) {
                                tempRev <- fisheryBySpecie[[fisIdx]]$LWstat[fisheryBySpecie[[fisIdx]]$LWstat$FG == i, ]
                                if (nrow(tempRev) > 0) {
                                  if (curUnit == "Length") {
                                    tempRev$propWei <- tempRev$relAbb / sum(tempRev$relAbb)
                                    tempRev$SizeClass <- factor(findInterval(x = tempRev$avgLen, vec = vecSize, all.inside = TRUE), levels = 1:(length(vecSize) - 1))
                                    outClass <- merge(data.frame(SizeClass = levels(tempRev$SizeClass)), aggregate(formula = propWei ~ SizeClass, data = tempRev, FUN = sum), all.x = TRUE)
                                    preRevenue[i, 2:length(vecSize)] <- outClass$propWei
                                  } else {
                                    tempRev$propWei <- tempRev$Freq / sum(tempRev$Freq)
                                    tempRev$SizeClass <- factor(findInterval(x = tempRev$Weight, vec = vecSize), levels = 1:length(vecSize))
                                    outClass <- merge(data.frame(SizeClass = levels(tempRev$SizeClass)), aggregate(formula = propWei ~ SizeClass, data = tempRev, FUN = sum), all.x = TRUE)
                                    preRevenue[i, 2:length(vecSize)] <- outClass$propWei
                                  }
                                }
                              }
                              outWeiProp[[fisheryBySpecie[[fisIdx]]$species]] <<- preRevenue
                              
                              
                              ## seasonal
                              outWeiPropQ[[fisheryBySpecie[[fisIdx]]$species]] <<- list()
                              for (season in c("winter", "spring", "summer", "fall")) {
                                preReveSea <- data.frame(FG = 1:length(fgNames))
                                preReveSea <- cbind(preReveSea, setNames(lapply(1:(length(vecSize) - 1), function(x) x <- NA), 1:(length(vecSize) - 1)))
                                for (i in 1:nrow(preReveSea)) {
                                  tempRev <- fisheryBySpecie[[fisIdx]]$LWstatQ[fisheryBySpecie[[fisIdx]]$LWstatQ$FG == i & fisheryBySpecie[[fisIdx]]$LWstatQ$Season == season, ]
                                  if (nrow(tempRev) > 0) {
                                    if (curUnit == "Length") {
                                      tempRev$propWei <- tempRev$relAbb / sum(tempRev$relAbb)
                                      tempRev$SizeClass <- factor(findInterval(x = tempRev$avgLen, vec = vecSize, all.inside = TRUE), levels = 1:(length(vecSize) - 1))
                                      outClass <- merge(data.frame(SizeClass = levels(tempRev$SizeClass)), aggregate(formula = propWei ~ SizeClass, data = tempRev, FUN = sum), all.x = TRUE)
                                      preReveSea[i, 2:length(vecSize)] <- outClass$propWei
                                    } else {
                                      tempRev$propWei <- tempRev$Freq / sum(tempRev$Freq)
                                      tempRev$SizeClass <- factor(findInterval(x = tempRev$Weight, vec = vecSize), levels = 1:length(vecSize))
                                      outClass <- merge(data.frame(SizeClass = levels(tempRev$SizeClass)), aggregate(formula = propWei ~ SizeClass, data = tempRev, FUN = sum), all.x = TRUE)
                                      preReveSea[i, 2:length(vecSize)] <- outClass$propWei
                                    }
                                  }
                                }
                                if ( fill == TRUE) {
                                  avgLW <- apply(preReveSea[,2:ncol(preReveSea)], 2, mean, na.rm = TRUE)/sum(apply(preReveSea[,2:ncol(preReveSea)], 2, mean, na.rm = TRUE))
                                  idxNArows <- which(apply(is.na(preReveSea[,2:ncol(preReveSea)]), 1, all))
                                  if (length(idxNArows) > 0) {
                                    preReveSea[idxNArows, 2:ncol(preReveSea)] <- rep(avgLW, each = length(idxNArows))
                                  } 
                                }
                                outWeiPropQ[[fisheryBySpecie[[fisIdx]]$species]][[season]] <<- preReveSea
                              }
                            }
                          },
                          getCostRevenue = function() {
                            simCostRevenue <<- merge(
                              simTotalCost[, c("I_NCEE", "Year", "totCost")],
                              simTotalRevenue[, c("I_NCEE", "Year", "totRevenue")]
                            )
                          },
                          simulateFishery = function(thr0 = 100, effoBan = numeric(0), timeStep = "Year", maxEffo = 0, fillLW = FALSE) {
                            cat("\nGetting length-weight statistics...", sep = "")
                            getLWstat(fill = fillLW)
                            cat("Done!", sep = "")
                            
                            cat("\nSetup initial parameters...", sep = "")
                            
                            genSimEffo()
                            simProdAll()
                            getSimTotalCost()
                            getSimRevenue(timeScale = timeStep)
                            getCostRevenue()
                            cat("Done!\n", sep = "")
                            
                            nFG <- sampMap$cutFG + 1
                            Esim <- Etemp <- simEffo
                            nVessels <- length(unique(simEffo$I_NCEE))
                            Gmat <- Pmat <- matrix(0, nVessels, 1)
                            Gmat[, 1] <- simCostRevenue$totRevenue - simCostRevenue$totCost
                            noV <- numeric(0)
                            nRec <- nrow(Esim)
                            nVproc <- nVessels
                            nIter <- 0
                            toOpt <- numeric(0)
                            
                            while (length(noV) < nVessels) {
                              nIter <- nIter + 1
                              cat("\nIteration", nIter)
                              
                              cat("\n\tOptimising effort... ", sep = "")
                              genSimEffo(selRow = toOpt, areaBan = effoBan)
                              cat("Done!", sep = "")
                              
                              if (maxEffo > 0) {
                                cat("\n\tScaling max effort... ", sep = "")
                                simEffo[, -c(1:3, ncol(simEffo))] <<- lvlSimEffo(simuEffo = simEffo[, -c(1:3, ncol(simEffo))], maxEff = maxEffo)
                                cat("Done!", sep = "")
                              }
                              cat("\n\tComputing production...", sep = "")
                              simProdAll(selRow = toOpt)
                              cat("Done!", sep = "")
                              
                              cat("\n\tComputing cost-revenues...", sep = "")
                              getSimTotalCost()
                              getSimRevenue(selRow = toOpt, timeScale = timeStep)
                              getCostRevenue()
                              cat("Done!", sep = "")
                              
                              EsimG <- simCostRevenue$totRevenue - simCostRevenue$totCost
                              Gmat <- cbind(Gmat, EsimG)
                              set_plus <- which(Gmat[, ncol(Gmat)] > Gmat[, ncol(Gmat) - 1])
                              set_minus <- setdiff(1:nrow(Gmat), set_plus)
                              pvec <- rep(1, nrow(Gmat))
                              if (length(set_plus) > 0) pvec[set_plus] <- 0
                              Pmat <- cbind(Pmat, pvec)
                              if (ncol(Pmat) >= thr0) {
                                noV <- unique(c(noV, which(apply(Pmat[, (ncol(Pmat) - thr0 + 1):ncol(Pmat)], 1, sum, na.rm = T) >= thr0)))
                                toOpt <- which(!(simEffo$I_NCEE %in% simCostRevenue$I_NCEE[noV]))
                              }
                              nVproc <- c(nVproc, nVessels - length(noV))
                              rec_minus <- which(simEffo$I_NCEE %in% simCostRevenue$I_NCEE[set_minus])
                              simEffo[rec_minus, ] <<- Etemp[rec_minus, ]
                              Etemp <- simEffo
                              Gmat[set_minus, ncol(Gmat)] <- Gmat[set_minus, ncol(Gmat) - 1]
                              
                              par(mfrow = c(1, 2), las = 2)
                              plot(1:ncol(Gmat), apply(Gmat, 2, sum, na.rm = TRUE) / 1000000,
                                   type = "l",
                                   xlab = "Iteration", ylab = "10^6 Euros", lwd = 3, col = 2
                              )
                              title(main = "Gains")
                              plot(1:ncol(Gmat), nVproc,
                                   type = "l",
                                   xlab = "Iteration", ylab = "", lwd = 3, col = 4
                              )
                              title(main = "Vessels to optimize")
                            }
                            
                            cat("\nSaving results...", sep = "")
                            outGmat <<- Gmat
                            # outNVlst <<- nVproc
                            outOptimEffo <<- Etemp
                            cat("Done!\n", sep = "")
                          },
                          setSimResults = function() {
                            simResPlot <<- list()
                            
                            outObsEffo <- fleet$effoAllLoa[fleet$effoAllLoa$Year == max(as.numeric(as.character(unique(fleet$effoAllLoa$Year)))), ]
                            cumObsEffo <- data.frame(
                              FG = colnames(outObsEffo[, 4:(ncol(outObsEffo) - 1)]),
                              Hours = apply(outObsEffo[, 4:(ncol(outObsEffo) - 1)], 2, sum)
                            )
                            rm(outObsEffo)
                            cumOptEffo <- data.frame(
                              FG = colnames(simEffo[, 4:(ncol(simEffo) - 1)]),
                              Hours = apply(simEffo[, 4:(ncol(simEffo) - 1)], 2, sum)
                            )
                            
                            deltaObsOpt <- data.frame(
                              x = as.numeric(colnames(simEffo[, 4:(ncol(simEffo) - 1)])),
                              obs = cumObsEffo$Hours,
                              opt = cumOptEffo$Hours,
                              delta = cumOptEffo$Hours - cumObsEffo$Hours,
                              deltaPerc = 100 * (cumOptEffo$Hours - cumObsEffo$Hours) / cumObsEffo$Hours
                            )
                            
                            all_cell <- merge(x = sampMap$cutResShpFort$id, deltaObsOpt, all = TRUE)
                            all_cell[is.na(all_cell)] <- 0
                            
                            grid_data <- cbind(sampMap$cutResShpFort, all_cell[, 2:5])
                            
                            if (length(sampMap$cutFG) > 0) {
                              if (length(sampMap$gridShp@polygons) == (sampMap$cutFG + 1)) {
                                tmp_coo <- data.frame(coordinates(sampMap$gridShp), cell_id = 1:length(sampMap$gridShp))
                                colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                              } else {
                                tmp_coo <- sampMap$cutResShpCent
                              }
                            } else {
                              tmp_coo <- sampMap$cutResShpCent
                            }
                            
                            grid_data$deltaPerc[grid_data$deltaPerc > 100] <- 101
                            grid_data$deltaPerc[grid_data$deltaPerc < -100] <- -101
                            grid_data$deltaPerc[grid_data$deltaPerc == 0] <- NA
                            grid_data$delta[grid_data$delta == 0] <- NA
                            
                            simResPlot[["obsEffort"]] <<- suppressMessages(sampMap$gooMapPlot +
                                                                             geom_polygon(aes(x = long, y = lat, group = group, fill = obs),
                                                                                          colour = "grey20", size = 0.1, data = grid_data, alpha = 0.8
                                                                             ) +
                                                                             scale_fill_gradient("Observed\nlog10(Hours)",
                                                                                                 low = "snow1",
                                                                                                 high = "#fc8d59", trans = "log10"
                                                                             ) +
                                                                             geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                                       data = tmp_coo, size = 2
                                                                             ) +
                                                                             ggtitle("Map of observed effort pattern") +
                                                                             xlab("Longitude") + ylab("Latitude") +
                                                                             theme(legend.position = "left"))
                            
                            simResPlot[["optEffort"]] <<- suppressMessages(sampMap$gooMapPlot +
                                                                             geom_polygon(aes(x = long, y = lat, group = group, fill = opt),
                                                                                          colour = "grey20", size = 0.1, data = grid_data, alpha = 0.8
                                                                             ) +
                                                                             scale_fill_gradient("Optimized\nlog10(Hours)",
                                                                                                 low = "snow1",
                                                                                                 high = "#fc8d59", trans = "log10"
                                                                             ) +
                                                                             geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                                       data = tmp_coo, size = 2
                                                                             ) +
                                                                             ggtitle("Map of optimized effort pattern") +
                                                                             xlab("Longitude") + ylab("Latitude"))
                            
                            simResPlot[["absChange"]] <<- suppressMessages(sampMap$gooMapPlot +
                                                                             geom_polygon(aes(x = long, y = lat, group = group, fill = delta),
                                                                                          colour = "black", size = 0.1, data = grid_data, alpha = 1
                                                                             ) +
                                                                             scale_fill_gradient2("Effort Delta\n(Hours)",
                                                                                                  low = "#91bfdb",
                                                                                                  high = "#fc8d59", mid = "#ffffbf", na.value = "grey20"
                                                                             ) +
                                                                             geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                                       data = tmp_coo, size = 2
                                                                             ) +
                                                                             ggtitle("Map of Absolute Change") +
                                                                             xlab("Longitude") + ylab("Latitude") +
                                                                             theme(legend.position = "left"))
                            
                            simResPlot[["relChange"]] <<- suppressMessages(sampMap$gooMapPlot +
                                                                             geom_polygon(aes(x = long, y = lat, group = group, fill = deltaPerc),
                                                                                          colour = "black", size = 0.1, data = grid_data, alpha = 1
                                                                             ) +
                                                                             scale_fill_gradient2("Effort Delta\n(%)",
                                                                                                  low = "#91bfdb",
                                                                                                  high = "#fc8d59", mid = "#ffffbf", na.value = "grey20"
                                                                             ) +
                                                                             geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                                       data = tmp_coo, size = 2
                                                                             ) +
                                                                             ggtitle("Map of Relative Change") +
                                                                             xlab("Longitude") + ylab("Latitude"))
                          },
                          ggplotFishingPoints = function(year) {
                            tmp_dat <- fleet$rawEffort[[year]][sample(1:nrow(fleet$rawEffort[[year]]), min(c(50000, nrow(fleet$rawEffort[[year]])))), c("LON", "LAT", "FishPoint")]
                            tmp_dat$Status <- factor(tmp_dat$FishPoint, levels = c("FALSE", "TRUE"), labels = c("Not fishing", "Fishing"))
                            ggEffFish <<- suppressMessages(sampMap$gooMapPlot +
                                                             geom_point(
                                                               data = tmp_dat,
                                                               aes(x = LON, y = LAT, color = Status), size = 0.25, alpha = 0.2
                                                             ) +
                                                             scale_colour_manual(values = c("coral", "darkseagreen1")) +
                                                             guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
                                                             ggtitle(paste("Sample fishing points - ", year, sep = "")) +
                                                             theme_tufte(base_size = 14, ticks = T) +
                                                             theme(
                                                               legend.position = "bottom",
                                                               axis.text.x = element_text(size = 8),
                                                               axis.title.x = element_text(size = 10),
                                                               panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                                               axis.text.y = element_text(size = 8),
                                                               axis.title.y = element_text(size = 10),
                                                               legend.text = element_text(size = 8),
                                                               legend.title = element_text(size = 10)
                                                             ))
                          },
                          setPlotBetaMeltYear = function(species, year) {
                            tmp_melt_sub <- subset(fleet$betaMeltYear[[species]], Year == year)
                            all_cell <- merge(
                              x = sampMap$cutResShpFort$id,
                              data.frame(
                                x = as.numeric(substr(
                                  as.character(tmp_melt_sub$FishGround), 4,
                                  nchar(as.character(tmp_melt_sub$FishGround))
                                )),
                                y = tmp_melt_sub$Productivity
                              ), all = TRUE
                            )
                            all_cell[is.na(all_cell)] <- 0
                            grid_data <- cbind(sampMap$cutResShpFort, Beta = all_cell[, 2])
                            
                            if (length(sampMap$cutFG) > 0) {
                              if (length(sampMap$gridShp@polygons) == (sampMap$cutFG + 1)) {
                                tmp_coo <- data.frame(coordinates(sampMap$gridShp), cell_id = 1:length(sampMap$gridShp))
                                colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                              } else {
                                tmp_coo <- sampMap$cutResShpCent
                              }
                            } else {
                              tmp_coo <- sampMap$cutResShpCent
                            }
                            
                            sampMap$ggBetaFGmap <<- suppressMessages(sampMap$gooMapPlot + geom_polygon(aes(x = long, y = lat, group = group, fill = Beta),
                                                                                                       colour = "black", size = 0.1,
                                                                                                       data = grid_data, alpha = 0.8
                            ) +
                              scale_fill_gradient("Beta\nValues", low = "lightyellow", high = "mediumseagreen") +
                              geom_text(aes(label = FG, x = Lon, y = Lat),
                                        data = tmp_coo, size = 2
                              ) +
                              ggtitle(paste("Betas x Fishing Ground - ", year, sep = "")) +
                              xlab("Longitude") + ylab("Latitude") +
                              theme_tufte(base_size = 14, ticks = F) +
                              theme(
                                legend.position = "right",
                                plot.title = element_text(size = 14),
                                axis.text.x = element_text(size = 8),
                                axis.title = element_blank(),
                                panel.grid = element_line(size = 0.1, linetype = 2, colour = "grey20"),
                                axis.text.y = element_text(size = 10),
                                axis.ticks.y = element_blank()
                              ))
                            sampMap$ggBetaFGbox <<- suppressMessages(ggplot(
                              fleet$betaMeltYear[[species]],
                              aes(
                                x = FishGround, y = Productivity,
                                group = FishGround
                              )
                            ) +
                              geom_boxplot() +
                              coord_flip() +
                              geom_point(
                                data = tmp_melt_sub,
                                aes(
                                  x = FishGround, y = Productivity,
                                  fill = Productivity, group = FishGround
                                ),
                                size = 2, shape = 21, color = "grey40"
                              ) +
                              geom_line(
                                data = tmp_melt_sub,
                                aes(x = FishGround, y = Productivity, group = Year),
                                color = "grey40"
                              ) +
                              scale_fill_gradient(low = "lightyellow", high = "mediumseagreen") +
                              xlab("Fishing Ground") +
                              theme_tufte(base_size = 14, ticks = F) +
                              theme(
                                legend.position = "none",
                                plot.title = element_text(size = 14),
                                axis.text.x = element_text(size = 8),
                                axis.title = element_blank(),
                                panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                axis.text.y = element_text(size = 10),
                                axis.ticks.y = element_blank()
                              ))
                          },
                          setPlotProdMeltYear = function(species, year) {
                            tmp_melt_sub <- subset(fleet$prodMeltYear[[species]], Year == year)
                            all_cell <- merge(
                              x = sampMap$cutResShpFort$id,
                              data.frame(
                                x = substr(as.character(tmp_melt_sub$FishGround), 4, nchar(as.character(tmp_melt_sub$FishGround))),
                                y = tmp_melt_sub$Production
                              ), all = TRUE
                            )
                            all_cell[is.na(all_cell)] <- 0
                            grid_data <- cbind(sampMap$cutResShpFort, Hours = all_cell[, 2])
                            
                            if (length(sampMap$cutFG) > 0) {
                              if (length(sampMap$gridShp@polygons) == (sampMap$cutFG + 1)) {
                                tmp_coo <- data.frame(coordinates(sampMap$gridShp), cell_id = 1:length(sampMap$gridShp))
                                colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                              } else {
                                tmp_coo <- sampMap$cutResShpCent
                              }
                            } else {
                              tmp_coo <- sampMap$cutResShpCent
                            }
                            
                            sampMap$ggProdFGmap <- suppressMessages(sampMap$gooMapPlot +
                                                                      geom_polygon(aes(x = long, y = lat, group = group, fill = Hours),
                                                                                   colour = "black", size = 0.1,
                                                                                   data = grid_data, alpha = 0.8
                                                                      ) +
                                                                      scale_fill_gradient("Production\nValues",
                                                                                          low = "lightyellow", high = "slateblue1"
                                                                      ) +
                                                                      geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                                data = tmp_coo, size = 2
                                                                      ) +
                                                                      ggtitle(paste("Production x Fishing Ground - ", year, sep = "")) +
                                                                      xlab("Longitude") + ylab("Latitude") +
                                                                      theme_tufte(base_size = 14, ticks = F) +
                                                                      theme(
                                                                        legend.position = "right",
                                                                        plot.title = element_text(size = 14),
                                                                        axis.text.x = element_text(size = 8),
                                                                        axis.title = element_blank(),
                                                                        panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                                        axis.text.y = element_text(size = 10),
                                                                        axis.ticks.y = element_blank()
                                                                      ))
                            sampMap$ggProdFGbox <<- suppressMessages(ggplot(fleet$prodMeltYear[[species]], aes(x = FishGround, y = Production, group = FishGround)) +
                                                                       geom_boxplot() +
                                                                       coord_flip() +
                                                                       geom_point(
                                                                         data = tmp_melt_sub, aes(
                                                                           x = FishGround, y = Production,
                                                                           fill = Production, group = FishGround
                                                                         ),
                                                                         size = 2, shape = 21, color = "grey40"
                                                                       ) +
                                                                       geom_line(
                                                                         data = tmp_melt_sub, aes(x = FishGround, y = Production, group = Year),
                                                                         color = "grey40"
                                                                       ) +
                                                                       scale_fill_gradient(low = "lightyellow", high = "slateblue1") +
                                                                       xlab("Fishing Ground") +
                                                                       theme_tufte(base_size = 14, ticks = F) +
                                                                       theme(
                                                                         legend.position = "none",
                                                                         plot.title = element_text(size = 14),
                                                                         axis.text.x = element_text(size = 8),
                                                                         axis.title = element_blank(),
                                                                         panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                                         axis.text.y = element_text(size = 10),
                                                                         axis.ticks.y = element_blank()
                                                                       ))
                          },
                          setCellPoin = function() {
                            sampMap$gridShp@plotOrder <- 1:length(sampMap$gridShp@plotOrder)
                            tmp_polygons <- SpatialPolygons(sampMap$gridShp@polygons)
                            cat("\nGridding year ", sep = "")
                            for (j in names(fleet$rawEffort)) {
                              cat(j, "... ", sep = "")
                              fleet$rawEffort[[j]]$Cell <<- over(SpatialPoints(fleet$rawEffort[[j]][, c("LON", "LAT")]), tmp_polygons)
                            }
                            cat("Done!", sep = "")
                          },
                          setTrackHarb = function() {
                            fleet$trackHarbs <<- list()
                            for (i in names(fleet$rawEffort)) {
                              cat("\nLoading effort year: ", i, "... ", sep = "")
                              tmp_eff <- fleet$rawEffort[[i]]
                              tmp_harbs <- sqldf("select xCFR I_NCEE, xTnum T_NUM, LON_S, LAT_S, DATE_S, LON_E, LAT_E, DATE_E from
                                                 (select x.I_NCEE xCFR, LAT LAT_S, LON LON_S, DATE DATE_S, T_NUM xTnum from tmp_eff x where W_HARB = 1)
                                                 join
                                                 (select y.I_NCEE yCFR, LAT LAT_E, LON LON_E, DATE DATE_E, T_NUM yTnum from tmp_eff y where W_HARB = 1)
                                                 on xCFR = yCFR and xTnum = yTnum and DATE_S != DATE_E and DATE_S < DATE_E order by xCFR, xTnum, DATE_S")
                              cat("Done!", sep = "")
                              uni_harbs <- as.data.frame(unique(rbind(as.matrix(tmp_harbs[, c("LON_S", "LAT_S")]), as.matrix(tmp_harbs[, c("LON_E", "LAT_E")]))))
                              uni_harbs$Name <- ""
                              cat("\nSetting nearest harbor name... ", sep = "")
                              tmp_dist <- apply(spDists(
                                x = as.matrix(uni_harbs[, 1:2]),
                                y = as.matrix(sampMap$harbDbf[, 1:2]), longlat = TRUE
                              ), 1, which.min)
                              uni_harbs$Name <- as.character(sampMap$harbDbf[tmp_dist, 3])
                              cat("Done!", sep = "")
                              cat("\nSaving... ", sep = "")
                              fleet$trackHarbs[[i]] <<- sqldf("select I_NCEE, T_NUM, k.LON_S LON_S, k.LAT_S LAT_S, DATE_S, HARB_S, LON_E, LAT_E, DATE_E, Name HARB_E from (select I_NCEE, T_NUM, LON_S, LAT_S, DATE_S, Name HARB_S, LON_E, LAT_E, DATE_E from tmp_harbs join uni_harbs using (LON_S, LAT_S)) k join uni_harbs x on x.LON_S = LON_E and x.LAT_S = LAT_E")
                              cat("Done!", sep = "")
                            }
                          },
                          setFishGround = function(numCut) {
                            sampMap$cutFG <<- numCut
                            sampMap$setCutResult(ind_clu = numCut)
                            simBanFG <<- data.frame(FG = sampMap$cutResShpCent$FG, Banned = "0", stringsAsFactors = FALSE)
                            tmp_clust <- cbind(
                              Cell = 1:sampMap$nCells,
                              FishGround = sampMap$clusMat[, numCut]
                            )
                            cat("\n\nSetting Fishing Ground year\n", sep = "")
                            for (j in names(fleet$rawEffort)) {
                              cat(j, "... ", sep = "")
                              fleet$rawEffort[[j]]$FishGround <<- tmp_clust[fleet$rawEffort[[j]]$Cell, 2]
                            }
                            cat("Done!\n", sep = "")
                          },
                          addFg2Fishery = function() {
                            cat("\n\nConnecting coordinates to fishing ground...", sep = "")
                            
                            rawDataFishery$numFG <<- names(sampMap$cutResShp)[over(
                              SpatialPoints(data.frame(
                                Lon = rawDataFishery$Lon,
                                Lat = rawDataFishery$Lat
                              ), proj4string = CRS(proj4string(sampMap$gridShp))),
                              sampMap$cutResShp
                            )]
                            
                            for (i in 1:length(fisheryBySpecie)) {
                              fisheryBySpecie[[i]]$rawLFD$numFG <<- names(sampMap$cutResShp)[over(
                                SpatialPoints(data.frame(
                                  Lon = fisheryBySpecie[[i]]$rawLFD$Lon,
                                  Lat = fisheryBySpecie[[i]]$rawLFD$Lat
                                ), proj4string = CRS(proj4string(sampMap$gridShp))),
                                sampMap$cutResShp
                              )]
                            }
                            
                            cat(" Done!\n", sep = "")
                          },
                          addFg2Survey = function() {
                            cat("\n\nConnecting coordinates to fishing ground...", sep = "")
                            rawDataSurvey$numFG <<- names(sampMap$cutResShp)[over(
                              SpatialPoints(data.frame(
                                Lon = rawDataSurvey$Lon,
                                Lat = rawDataSurvey$Lat
                              ), proj4string = CRS(proj4string(sampMap$gridShp))),
                              sampMap$cutResShp
                            )]
                            
                            for (i in 1:length(surveyBySpecie)) {
                              surveyBySpecie[[i]]$rawLFD$numFG <<- names(sampMap$cutResShp)[over(
                                SpatialPoints(data.frame(
                                  Lon = surveyBySpecie[[i]]$rawLFD$Lon,
                                  Lat = surveyBySpecie[[i]]$rawLFD$Lat
                                ), proj4string = CRS(proj4string(sampMap$gridShp))),
                                sampMap$cutResShp
                              )]
                            }
                            cat(" Done!\n", sep = "")
                          },
                          setWeekEffoMatrCell = function() {
                            fleet$weekEffoMatr <<- list()
                            for (j in names(fleet$rawEffort)) {
                              cat("\n\nLoading year ", j, " ... ", sep = "")
                              tmp_dat <- fleet$rawEffort[[j]][fleet$rawEffort[[j]]$FishPoint == TRUE & !is.na(fleet$rawEffort[[j]]$Cell), c("I_NCEE", "T_NUM", "WeekNum", "Cell", "FishPoint")]
                              cat("Done!", sep = "")
                              tmp_dat$Cell <- as.factor(tmp_dat$Cell)
                              cat("\nCreating weekly fishing effort matrix... ", sep = "")
                              tmp_matrix <- dcast(tmp_dat,
                                                  I_NCEE + T_NUM + WeekNum ~ Cell,
                                                  fun.aggregate = sum,
                                                  na.rm = TRUE, value.var = "FishPoint"
                              )
                              cat("Done!", sep = "")
                              cat("\nChecking... ", sep = "")
                              miss_cols <- setdiff(as.character(sampMap$gridShp@plotOrder), names(tmp_matrix)[4:ncol(tmp_matrix)])
                              if (length(miss_cols) > 0) {
                                cat(length(miss_cols), " cells with no points... ", sep = "")
                                tmp_matrix[, miss_cols] <- 0
                                tmp_matrix <- tmp_matrix[, c(1:3, 3 + order(as.numeric(names(tmp_matrix)[4:ncol(tmp_matrix)])))]
                              }
                              cat(" Done!", sep = "")
                              fleet$weekEffoMatr[[j]] <<- tmp_matrix
                            }
                          },
                          setWeekEffoMatrGround = function() {
                            fleet$weekEffoMatr <<- list()
                            for (j in names(fleet$rawEffort)) {
                              cat("\n\nLoading year ", j, " ... ", sep = "")
                              tmp_dat <- fleet$rawEffort[[j]][fleet$rawEffort[[j]]$FishPoint == TRUE & !is.na(fleet$rawEffort[[j]]$Cell), c("I_NCEE", "T_NUM", "WeekNum", "MonthNum", "FishGround", "FishPoint")]
                              cat("Done!", sep = "")
                              tmp_dat$FishGround <- as.factor(tmp_dat$FishGround)
                              cat("\nCreating weekly fishing effort matrix... ", sep = "")
                              tmp_matrix <- dcast(tmp_dat,
                                                  I_NCEE + T_NUM + WeekNum + MonthNum ~ FishGround,
                                                  fun.aggregate = sum,
                                                  na.rm = TRUE, value.var = "FishPoint"
                              )
                              cat("Done!", sep = "")
                              cat("\nChecking... ", sep = "")
                              miss_cols <- setdiff(
                                as.character(unique(fleet$rawEffort[[j]]$FishGround[!is.na(fleet$rawEffort[[j]]$FishGround)])),
                                names(tmp_matrix)[5:ncol(tmp_matrix)]
                              )
                              if (length(miss_cols) > 0) {
                                cat(length(miss_cols), " cells with no points... ", sep = "")
                                tmp_matrix[, miss_cols] <- 0
                                tmp_matrix <- tmp_matrix[, c(1:4, 4 + order(as.numeric(names(tmp_matrix)[5:ncol(tmp_matrix)])))]
                              }
                              tmp_matrix <- sqldf("select * from tmp_matrix join (select I_NCEE, T_NUM, min(DATE) DATE_S, max(DATE) DATE_E from tmp_dat group by I_NCEE, T_NUM) using (I_NCEE, T_NUM)")
                              cat(" Done!", sep = "")
                              fleet$weekEffoMatr[[j]] <<- tmp_matrix
                            }
                          },
                          ggplotGridEffort = function(year) {
                            tmp_dat <- table(fleet$rawEffort[[year]][fleet$rawEffort[[year]]$FishPoint == TRUE & !is.na(fleet$rawEffort[[year]]$Cell), c("Cell")])
                            all_cell <- merge(
                              x = sampMap$gridPolySet$PID,
                              data.frame(x = as.numeric(names(tmp_dat)), y = tmp_dat), all = TRUE
                            )[, c(1, 3)]
                            all_cell[is.na(all_cell)] <- 0
                            # grid_data <- cbind(sampMap$gridPolySet, LogCount = log10(all_cell[,2] + 1))
                            grid_data <- cbind(sampMap$gridPolySet, Count = all_cell[, 2])
                            ggEffGrid <<- suppressMessages(sampMap$gooMapPlot + geom_polygon(aes(x = X, y = Y, group = PID, fill = Count),
                                                                                             size = 0.2,
                                                                                             data = grid_data, alpha = 0.8
                            ) +
                              scale_fill_gradient(
                                low = "Yellow", high = "coral",
                                trans = "log10",
                                breaks = trans_breaks("log10", function(x) 10^x),
                                labels = trans_format("log10", math_format(10^.x))
                              ) +
                              ggtitle(paste("Fishing Effort - ", year, sep = "")) +
                              theme_tufte(base_size = 14, ticks = T) +
                              theme(
                                legend.position = "bottom",
                                axis.text.x = element_text(size = 8),
                                axis.title.x = element_text(size = 10),
                                panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                axis.text.y = element_text(size = 8),
                                axis.title.y = element_text(size = 10),
                                legend.text = element_text(size = 8),
                                legend.title = element_text(size = 10)
                              ))
                          },
                          getNnlsModel = function(species, minobs, thr_r2) {
                            data <- fleet$effoProdAllLoa
                            nFG <- sampMap$cutFG + 1
                            thrB <- fleet$specSett[[species]]$threshold
                            thr0 <- mean(data[, species][data[, species] < thrB & data[, species] > 0])
                            naFind <- which(is.na(data) == TRUE, arr.ind = TRUE)
                            if (length(naFind) > 0) data <- data[-naFind[, 1], ]
                            norow <- which(data[, which(colnames(data) == species)] <= thrB)
                            X0 <- data[-norow, which(colnames(data) %in% c("Year", "MonthNum", "Loa", as.character(seq(1:nFG))))]
                            Y0 <- data[-norow, which(colnames(data) == species)]
                            
                            # Gen Matrix of scenaria and fill by membership
                            nY <- length(unique(X0$Year))
                            nM <- 12
                            membY <- as.numeric(as.factor(X0$Year)) # Membership delle osservazioni x anno
                            membM <- as.numeric(X0$MonthNum) # Membership delle osservazioni x mese
                            SceMat <- expand.grid(unique(membY), unique(membM))
                            colnames(SceMat) <- c("YEAR", "MONTH")
                            SceList <- vector(mode = "list", length = nrow(SceMat))
                            il <- 0
                            for (i in 1:nrow(SceMat)) {
                              il <- il + 1
                              SceList[[il]] <- which((membY == SceMat[il, 1]) & (membM == SceMat[il, 2]))
                            }
                            lsce <- as.numeric(lapply(SceList, length))
                            wsce <- which(lsce >= minobs)
                            nSce <- length(wsce)
                            
                            # Gen and fill the matrix of betas
                            bmat <- matrix(NA, length(unique(membY)) * length(unique(membM)), nFG)
                            obsY <- fittedY <- vector(mode = "list", length = nSce)
                            nnls_r2 <- rep(NA, nSce)
                            nfitted <- 0
                            nno <- 0
                            wcol <- which(colnames(X0) %in% as.character(seq(1, nFG)))
                            
                            for (iS in wsce) {
                              subX <- X0[SceList[[iS]], wcol] * (X0$Loa[SceList[[iS]]])
                              subY <- Y0[SceList[[iS]]]
                              zeroFG <- which(apply(subX, 2, sum) == 0)
                              eFG <- setdiff(1:nFG, zeroFG)
                              if (length(zeroFG) > 0) subX <- subX[, -zeroFG]
                              if (min(dim(subX)) > 0) {
                                nnls_m <- getNNLS(subX, subY, zeroFG)
                                if (abs(nnls_m$r2) > thr_r2) {
                                  nnls_r2[iS] <- nnls_m$r2
                                  obsY[[iS]] <- nnls_m$obs
                                  fittedY[[iS]] <- nnls_m$fitted
                                  tcoo <- 12 * (SceMat[iS, "YEAR"] - 1) + SceMat[iS, "MONTH"]
                                  bmat[tcoo, eFG] <- nnls_m$betas
                                  if (length(zeroFG) > 0) bmat[tcoo, zeroFG] <- 0
                                  nfitted <- nfitted + 1
                                } else {
                                  nno <- nno + 1
                                }
                              } else {
                                nno <- nno + 1
                              }
                            }
                            cat("\nNNLS: ", nSce, " actual scenarios - ", nfitted, " fitted", "(", floor(100 * (nSce - nno) / nSce), "%)", sep = "")
                            blist <- vector(mode = "list", length = 4)
                            colnames(bmat) <- paste("BE_", formatC(1:ncol(bmat), width = nchar(ncol(bmat)), format = "d", flag = "0"), sep = "")
                            
                            if (anyNA(bmat)) {
                              zero_chk <- which(apply(bmat, 2, sum, na.rm = TRUE) == 0)
                              if (length(zero_chk) > 0) {
                                par_betas <- fillbetas(bmat[, -zero_chk])
                                zero_beta <- matrix(0,
                                                    nrow = nrow(data.frame(bmat[, zero_chk])),
                                                    ncol = ncol(data.frame(bmat[, zero_chk]))
                                )
                                colnames(zero_beta) <- colnames(bmat)[zero_chk]
                                colnames(par_betas) <- colnames(bmat)[-zero_chk]
                                full_betas <- cbind(par_betas, zero_beta)
                                bmat <- full_betas[, order(colnames(full_betas))]
                              } else {
                                bmat <- fillbetas(bmat)
                              }
                            }
                            
                            if (any(bmat < 0)) bmat[which(bmat < 0)] <- 0
                            
                            blist[[1]] <- bmat
                            blist[[2]] <- unlist(obsY)
                            blist[[3]] <- unlist(fittedY)
                            blist[[4]] <- unlist(nnls_r2)
                            blist[[5]] <- SceMat
                            blist[[6]] <- nfitted
                            blist[[7]] <- nSce
                            names(blist) <- c("bmat", "obsY", "fittedY", "nnls_r2", "SceMat", "nfitted", "nSce")
                            fleet$resNNLS[[species]] <<- blist
                            fleet$setBetaMeltYear(species)
                          },
                          cohoDisPlot = function(whoSpe, whoCoh, whiYea, interp) {
                            if (interp == FALSE) {
                              if (whoCoh == "All") {
                                if (whiYea == "All") {
                                  # 1+round(apply(surveyBySpecie[[whoSpe]]$Coh_A[,,,],1,sum)/max(apply(surveyBySpecie[[whoSpe]]$Coh_A[,,,],1,sum)), 2)*100
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A[, , , ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - All cohorts - All years", sep = ""), legendUnits = "N."
                                  )
                                } else {
                                  # 1+round(apply(surveyBySpecie[[whoSpe]]$Coh_A[,,whiYea,],1,sum)/max(apply(surveyBySpecie[[whoSpe]]$Coh_A[,,whiYea,],1,sum)), 2)*100
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A[, , whiYea, ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - All cohorts - Year: ", whiYea, sep = ""),
                                    legendUnits = "N."
                                  )
                                }
                              } else {
                                if (whiYea == "All") {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A[, whoCoh, , ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - Cohort: ", whoCoh, "- All years", sep = ""),
                                    legendUnits = "N."
                                  )
                                } else {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A[, whoCoh, whiYea, ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - Cohort: ", whoCoh, " - Year: ", whiYea, sep = ""),
                                    legendUnits = "N."
                                  )
                                }
                              }
                            } else {
                              if (whoCoh == "All") {
                                if (whiYea == "All") {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A_Int[, , , ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - All cohorts - All years", sep = ""),
                                    legendUnits = "N."
                                  )
                                } else {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A_Int[, , whiYea, ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - All cohorts - Year: ", whiYea, sep = ""),
                                    legendUnits = "N."
                                  )
                                }
                              } else {
                                if (whiYea == "All") {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A_Int[, whoCoh, , ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - Cohort: ", whoCoh, "- All years", sep = ""),
                                    legendUnits = "N."
                                  )
                                } else {
                                  yea_abb <- round(apply(surveyBySpecie[[whoSpe]]$Coh_A_Int[, whoCoh, whiYea, ], 1, sum))
                                  round_yea <- 1 + 100 * yea_abb / max(yea_abb)
                                  
                                  distrPlotCols(
                                    cols = rev(topo.colors(101)), vals = round_yea,
                                    maxVal = ceiling(max(yea_abb)),
                                    plotTitle = paste("Species: ", specieInSurvey[whoSpe], " - Cohort: ", whoCoh, " - Year: ", whiYea, sep = ""),
                                    legendUnits = "N."
                                  )
                                }
                              }
                            }
                          },
                          # setCoh_A_Survey = function(){
                          #   if(length(specieInSurvey) == 1){
                          #     calcCoh_A_Survey(1)
                          #   }else{
                          #     for(i in 1:length(specieInSurvey)){
                          #       calcCoh_A_Survey(i)
                          #     }}
                          # },
                          # setCoh_A_Fishery = function(){
                          #   if(length(specieInFishery) == 1){
                          #     calcCoh_A_Fishery(1)
                          #   }else{
                          #     for(i in 1:length(specieInFishery)){
                          #       calcCoh_A_Fishery(i)
                          #     }}
                          # },
                          calcCoh_A_Survey = function(ind_num) {
                            Pop <- surveyBySpecie[[ind_num]]$LFDPop
                            LC <- surveyBySpecie[[ind_num]]$lengClas[-length(surveyBySpecie[[ind_num]]$lengClas)]
                            sp <- surveyBySpecie[[ind_num]]$species
                            nc <- surveyBySpecie[[ind_num]]$nCoho
                            surveyBySpecie[[ind_num]]$Coh_A <<- array(dim = c(sampMap$nCells, nc, length(surveyBySpecie[[ind_num]]$year), 2))
                            for (y in 1:length(surveyBySpecie[[ind_num]]$year)) {
                              for (sex in c(1:2)) {
                                mms <- surveyBySpecie[[ind_num]]$mixPar[[sex]][[1]][y, ]
                                sds <- surveyBySpecie[[ind_num]]$mixPar[[sex]][[2]][y, ]
                                opt <- matrix(0, length(LC), nc)
                                for (ij in 1:sampMap$nCells) {
                                  vv <- Pop[ij, , y, sex]
                                  coh.abb <- numeric(nc)
                                  if (sum(vv) > 0) {
                                    for (coh in c(1:nc)) opt[, coh] <- dnorm(LC, mms[coh], sds[coh])
                                    opt.ass <- apply(opt, 1, which.max)
                                    for (coh in c(1:nc)) coh.abb[coh] <- sum(vv[which(opt.ass == coh)])
                                  }
                                  surveyBySpecie[[ind_num]]$Coh_A[ij, 1:nc, y, sex] <- as.numeric(coh.abb)
                                }
                              }
                            }
                          },
                          calcCoh_A_Fishery = function(ind_num) {
                            Pop <- fisheryBySpecie[[ind_num]]$LFDPop
                            LC <- fisheryBySpecie[[ind_num]]$lengClas[-length(fisheryBySpecie[[ind_num]]$lengClas)]
                            sp <- fisheryBySpecie[[ind_num]]$species
                            nc <- fisheryBySpecie[[ind_num]]$nCoho
                            fisheryBySpecie[[ind_num]]$Coh_A <<- array(dim = c(sampMap$nCells, nc, length(fisheryBySpecie[[ind_num]]$year), 2))
                            for (y in 1:length(fisheryBySpecie[[ind_num]]$year)) {
                              for (sex in c(1:2)) {
                                mms <- fisheryBySpecie[[ind_num]]$mixPar[[sex]][[1]][y, ]
                                sds <- fisheryBySpecie[[ind_num]]$mixPar[[sex]][[2]][y, ]
                                opt <- matrix(0, length(LC), nc)
                                for (ij in 1:sampMap$nCells) {
                                  vv <- Pop[ij, , y, sex]
                                  coh.abb <- numeric(nc)
                                  if (sum(vv) > 0) {
                                    for (coh in c(1:nc)) opt[, coh] <- dnorm(LC, mms[coh], sds[coh])
                                    opt.ass <- apply(opt, 1, which.max)
                                    for (coh in c(1:nc)) coh.abb[coh] <- sum(vv[which(opt.ass == coh)])
                                  }
                                  fisheryBySpecie[[ind_num]]$Coh_A[ij, 1:nc, y, sex] <- as.numeric(coh.abb)
                                }
                              }
                            }
                          },
                          intrpCoh_A_Survey = function(ind_num) {
                            surveyBySpecie[[ind_num]]$Coh_A_Int <<- array(dim = c(sampMap$nCells, surveyBySpecie[[ind_num]]$nCoho, length(surveyBySpecie[[ind_num]]$year), 2))
                            for (y in 1:length(surveyBySpecie[[ind_num]]$year)) {
                              for (sex in 1:2) {
                                for (coh in 1:surveyBySpecie[[ind_num]]$nCoho) {
                                  xdata <- cbind(sampMap$griCent, surveyBySpecie[[ind_num]]$Coh_A[, coh, y, sex])
                                  colnames(xdata) <- c("Lon", "Lat", "Coh")
                                  xdata <- as.data.frame(xdata)
                                  yea_poi <- surveyBySpecie[[ind_num]]$rawLFD[which(surveyBySpecie[[ind_num]]$rawLFD$Year == surveyBySpecie[[ind_num]]$year[y]), c("Lon", "Lat")]
                                  cMEDITS <- which(!is.na(over(sampMap$gridShp, SpatialPoints(unique(yea_poi)))))
                                  noMEDITS <- setdiff(c(1:sampMap$nCells), cMEDITS)
                                  Areacell <- 9.091279 * 11.112
                                  RateArea <- Areacell / 100
                                  surveyBySpecie[[ind_num]]$Coh_A_Int[, coh, y, sex] <- IntInvDis(RateArea * xdata, cMEDITS, noMEDITS,
                                                                                                  Refmax = 5, Refmin = 3,
                                                                                                  sampMap$nCells,
                                                                                                  sampMap$gridShp, graph = T, logplot = F
                                  )[, 3]
                                }
                              }
                            }
                          },
                          intrpCoh_A_Fishery = function(ind_num) {
                            fisheryBySpecie[[ind_num]]$Coh_A_Int <<- array(dim = c(sampMap$nCells, fisheryBySpecie[[ind_num]]$nCoho, length(fisheryBySpecie[[ind_num]]$year), 2))
                            for (y in 1:length(fisheryBySpecie[[ind_num]]$year)) {
                              for (sex in 1:2) {
                                for (coh in 1:fisheryBySpecie[[ind_num]]$nCoho) {
                                  xdata <- cbind(sampMap$griCent, fisheryBySpecie[[ind_num]]$Coh_A[, coh, y, sex])
                                  colnames(xdata) <- c("Lon", "Lat", "Coh")
                                  xdata <- as.data.frame(xdata)
                                  yea_poi <- fisheryBySpecie[[ind_num]]$rawLFD[which(fisheryBySpecie[[ind_num]]$rawLFD$Year == fisheryBySpecie[[ind_num]]$year[y]), c("Lon", "Lat")]
                                  cMEDITS <- which(!is.na(over(sampMap$gridShp, SpatialPoints(unique(yea_poi)))))
                                  noMEDITS <- setdiff(c(1:sampMap$nCells), cMEDITS)
                                  Areacell <- 9.091279 * 11.112
                                  RateArea <- Areacell / 100
                                  fisheryBySpecie[[ind_num]]$Coh_A_Int[, coh, y, sex] <- IntInvDis(RateArea * xdata, cMEDITS, noMEDITS,
                                                                                                   Refmax = 5, Refmin = 3,
                                                                                                   sampMap$nCells,
                                                                                                   sampMap$gridShp, graph = T, logplot = F
                                  )[, 3]
                                }
                              }
                            }
                          }
                        )
)



#### SurveyBySpecie#################################################
#' SurveyBySpecie
#'
#' The \code{SurveyBySpecie} class implements the class of SMART to
#'  handle species samplings.
#'
#' @docType class
#' @usage NULL
#' @keywords data
#' @return Object of \code{\link{R6Class}} with attributes and methods for the survey data.
#'
#' @format \code{\link{R6Class}} object.
#'
#' @field species Name of the species.
#' @field year Years in the time-serie.
#' @field rawLFD data.frame, raw length frequency distribution.
#' @field abuAvg data.frame, average abundances by depth' stratum.
#' @field meditsIndex data.frame, medits index by depth' stratum.
#' @field lengClas numeric, length classes.
#' @field nCoho numeric, number of cohorts.
#' @field spreDist list of DF, lfd by sex.
#' @field sprePlot plots of LFD statistics.
#' @field spreSpat list of DF, spatial distribution by sex.
#' @field sampMcmc list, mcmc output chains.
#' @field groMixout list of DF, aged individuals by sex.
#' @field groPars list of DF, growth parameters by sex.
#' @field LWpar list of DF, length/weight parameters by sex.
#'
#' @section Methods:
#' \describe{
#'   \item{\code{initialize(sing_spe)}}{Automatic initialization made by the
#'   SmartProject class}
#'   \item{\code{setRawData(raw_data)}}{This method is used load the initial
#'   raw dataset}
#'   \item{\code{setYears()}}{This method is used to store the years in the
#'   provided time-serie}
#'   \item{\code{setSpecie()}}{This method is used to store the name of
#'   the species of the initial raw data}
#'   \item{\code{setLClass()}}{This method is used to store the unique
#'   length values of the sampled species}
#'   \item{\code{setDepth(bathyMatrix)}}{This method is used to assign the
#'   depth value corresponding to each sampling location}
#'   \item{\code{setStratum(vecStrata)}}{This method is used to set the
#'   depth strata of each sampling location}
#'   \item{\code{setIndSpe()}}{This method is used to aggregate the abundance
#'   data into the medits index}
#'   \item{\code{setAbuAvg()}}{This method is used to standardize the
#'   spatial abundances by depth strata}
#'   \item{\code{setNCoho(num_coh)}}{This method is used to setup the number
#'   of cohorts for the ageing module}
#'   \item{\code{setLWpar(alphaVal, betaVal, sex)}}{This method is used to
#'   store the alpha and beta values for the length/weight relationship}
#'   \item{\code{setWeight(sexVal = "Female")}}{This method is used to
#'   compute the fish weight given their length and the LWrelationship}
#'   \item{\code{setSpreDistSing()}}{This method is used to spread the
#'   aggregated LFD abundances into single individuals}
#'   \item{\code{setSprePlot(sampSex)}}{This method is used to setup the
#'   plots of the LFD statistics}
#'   \item{\code{setSpatDistSing()}}{This method is used to setup the spatial
#'   distribution of the single specimens}
#'   \item{\code{setSpatPlot(sampSex)}}{This method is used to store the
#'   spatial plots of the population}
#'   \item{\code{getMCsamps(numSamp, numAdap, numIter, sexDrop, curveSel)}}{
#'   This method is used to get a sample of the population to feed the mcmc
#'   module}
#'   \item{\code{getGrowPar(sexDrop)}}{This method is used to extract the
#'   growth parameters from the mcmc results}
#'   \item{\code{getMCage(sexDrop)}}{This method is used to assign an age to
#'   each fish}
#'   \item{\code{setMCplot(sexDrop, selCurve)}}{This method is used to setup
#'   the plot of the mcmc results}
#'   \item{\code{calcMixDate(nAdap, nSamp, nIter, sexDrop, curveSel)}}{
#'   This method is used to estimate the growth parameters of a population}
#'   \item{\code{ggplotMcmcOut(selCompo, selSex)}}{This method is used to
#'   output the stored plots of mcmc results}
#'   }


SurveyBySpecie <- R6Class("SurveyBySpecie",
                          portable = FALSE,
                          class = TRUE,
                          public = list(
                            species = NULL,
                            year = NULL,
                            rawLFD = NULL,
                            abuAvg = NULL,
                            meditsIndex = NULL,
                            lengClas = NULL,
                            LFDPop = NULL,
                            mixPar = NULL,
                            nCoho = NULL,
                            speSex = NULL,
                            spreDist = list(),
                            sprePlot = list(),
                            spreSpat = list(),
                            sampMcmc = list(),
                            groMixout = list(),
                            groPars = list(),
                            Coh_A = NULL,
                            Coh_A_Int = NULL,
                            LWpar = NULL,
                            initialize = function(sing_spe) {
                              setRawData(sing_spe)
                              setYears()
                              setSpecie()
                              setLClass()
                            },
                            setRawData = function(raw_data) {
                              rawLFD <<- raw_data
                            },
                            plotLFD = function() {
                              plotSpeAllYea(rawLFD)
                            },
                            setYears = function() {
                              year <<- sort(unique(rawLFD[, "Date"]), decreasing = FALSE)
                            },
                            setSpecie = function() {
                              species <<- unique(rawLFD[, "Species"])
                            },
                            setLClass = function() {
                              lengClas <<- seq(from = min(rawLFD[, "Class"]), to = max(rawLFD[, "Class"]), by = 1)
                            },
                            setDepth = function(bathyMatrix) {
                              rawLFD$Depth <<- get.depth(bathyMatrix, x = rawLFD$Lon, y = rawLFD$Lat, locator = FALSE)[, 3]
                            },
                            setStratum = function(vecStrata = c(0, 10, 100, 1000, Inf)) {
                              tmp_mem <- findInterval(x = -rawLFD$Depth, vec = vecStrata)
                              rawLFD$Stratum <<- factor(tmp_mem, levels = 1:(length(vecStrata) - 1), labels = paste(vecStrata[-length(vecStrata)], vecStrata[-1], sep = " - "))
                            },
                            setIndSpe = function() {
                              meditsIndex <<- aggregate(weiFem ~ Class + Year, data = abuAvg, sum)
                              meditsIndex <<- merge(x = meditsIndex, y = aggregate(weiMal ~ Class + Year, data = abuAvg, sum), all = TRUE)
                              meditsIndex <<- merge(x = meditsIndex, y = aggregate(weiUns ~ Class + Year, data = abuAvg, sum), all = TRUE)
                            },
                            setAbuAvg = function() {
                              tmp_lfd <- rawLFD
                              tmp_lfd$Year <- years(tmp_lfd$Date)
                              if ("Area" %in% colnames(tmp_lfd)) {
                                tmp_lfd$Female <- tmp_lfd$Female / tmp_lfd$Area
                                tmp_lfd$Male <- tmp_lfd$Male / tmp_lfd$Area
                                tmp_lfd$Unsex <- tmp_lfd$Unsex / tmp_lfd$Area
                              }
                              
                              tmp_aggFem <- aggregate(Female ~ Class + Stratum + Year, data = tmp_lfd, mean)
                              tmp_aggMal <- aggregate(Male ~ Class + Stratum + Year, data = tmp_lfd, mean)
                              tmp_aggUns <- aggregate(Unsex ~ Class + Stratum + Year, data = tmp_lfd, mean)
                              tmp_all <- merge(x = tmp_aggFem, y = tmp_aggMal, all = TRUE)
                              tmp_all <- merge(x = tmp_all, y = tmp_aggUns, all = TRUE)
                              abuAvg <<- tmp_all
                            },
                            setNCoho = function(num_coh) {
                              nCoho <<- num_coh
                            },
                            setLWpar = function(alphaVal, betaVal, sex) {
                              LWpar[[sex]] <<- list(alpha = as.numeric(alphaVal), beta = as.numeric(betaVal))
                            },
                            setWeight = function(sexVal = "Female") {
                              groMixout[[sexVal]]$Weight <<- LWpar[[sexVal]][["alpha"]] * groMixout[[sexVal]]$Length^LWpar[[sexVal]][["beta"]]
                            },
                            setSpreDistSing = function() {
                              for (sex in c("Female", "Male", "Unsex")) {
                                tmp_spre <- rawLFD[!is.na(rawLFD$numFG), c("Date", "Class", "numFG", sex)]
                                
                                num_sex <- sum(tmp_spre[, 4])
                                cat("Found", num_sex, sex, as.character(species), "samples\n", sep = " ")
                                
                                spreDist <- data.frame(
                                  UTC = rep(tmp_spre$Date, tmp_spre[, 4]),
                                  Length = rep(tmp_spre$Class, tmp_spre[, 4]) + runif(num_sex, -0.5, 0.5),
                                  NumFG = rep(tmp_spre$numFG, tmp_spre[, 4])
                                )
                                
                                spreDist$Year <- years(spreDist$UTC)
                                spreDist$Month <- months(as.chron(spreDist$UTC))
                                
                                spreDist[[sex]] <<- spreDist
                                
                                setSprePlot(sampSex = sex)
                              }
                            },
                            setAvailSex = function() {
                              speSex <<- sort(names(which(lapply(spreDist, nrow) > 0)))
                            },
                            setSprePlot = function(sampSex) {
                              sprePlot[[sampSex]] <<- list(
                                histLfdTot = set_ggHistLfdTot(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                histUtcTot = set_ggHistUtcTot(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                dotUtcSplit = set_ggDotUtcSplit(spreDist[[sampSex]]) + scale_color_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                histUtcLfd = set_ggHistUtcLfd(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE")))
                              )
                            },
                            setSpatDistSing = function() {
                              for (sex in c("Female", "Male", "Unsex")) {
                                tmp_fishSpat <- rawLFD[!is.na(rawLFD$numFG) & rawLFD[, sex] > 0, c("Lon", "Lat", "numFG", sex)]
                                if (nrow(tmp_fishSpat) > 0) {
                                  barploFgAll <- data.frame(table(tmp_fishSpat$numFG))
                                  barploFgAll <- barploFgAll[order(as.numeric(as.character(barploFgAll[, 1]))), ]
                                  barploFgAll$FG <- factor(barploFgAll$Var1, levels = barploFgAll$Var1)
                                  barploFgAll$relFreq <- round(100 * barploFgAll$Freq / sum(barploFgAll$Freq), 1)
                                  spreSpat[[sex]] <<- barploFgAll
                                  setSpatPlot(sampSex = sex)
                                }
                              }
                            },
                            setSpatPlot = function(sampSex) {
                              sprePlot[[sampSex]][["spatAbbTbl"]] <<- set_spatAbbTbl(spreSpat[[sampSex]])
                              sprePlot[[sampSex]][["spatAbsFreq"]] <<- set_spatAbsFreq(spreSpat[[sampSex]])
                              sprePlot[[sampSex]][["spatRelFreq"]] <<- set_spatRelFreq(spreSpat[[sampSex]])
                            },
                            getMCsamps = function(numSamp = 2000, numAdap = 100, numIter = 500, sexDrop = "Female", curveSel = "von Bertalanffy") {
                              sub_idx <- sample(1:nrow(spreDist[[sexDrop]]), size = numSamp, replace = ifelse(numSamp < nrow(spreDist[[sexDrop]]), FALSE, TRUE))
                              sub_data <- spreDist[[sexDrop]][sub_idx, ]
                              
                              N <- length(sub_data$Length)
                              alpha <- rep(1, nCoho)
                              Z <- rep(NA, N)
                              Z[which.min(sub_data$Length)] <- 1
                              Z[which.max(sub_data$Length)] <- nCoho
                              
                              dataList <- list(
                                y = sub_data$Length,
                                maxLeng = max(sub_data$Length), ## !!!
                                alpha = alpha,
                                Z = Z,
                                N = N,
                                Nclust = nCoho
                              )
                              
                              inits <- list(
                                list(Linf = min(sub_data$Length), k = 0.5, t0 = 0.0),
                                list(Linf = mean(sub_data$Length), k = 0.5, t0 = 0.0),
                                list(Linf = max(sub_data$Length), k = 0.5, t0 = 0.0)
                              )
                              
                              modelGrow <- ifelse(curveSel == "von Bertalanffy",
                                                  system.file("model/bertGrow.jags", package = "smartR"),
                                                  system.file("model/gompGrow.jags", package = "smartR")
                              )
                              
                              jags.m <- jags.model(modelGrow,
                                                   data = dataList,
                                                   inits = inits,
                                                   n.chains = 3,
                                                   n.adapt = numAdap
                              )
                              
                              ### MCMC chain sampling
                              # n.iter <- 500
                              obsNode <- c("Linf", "k", "t0", "tau", "p")
                              samps <- coda.samples(jags.m, obsNode, n.iter = numIter)
                              
                              sampMcmc[[sexDrop]] <<- samps
                            },
                            getGrowPar = function(sexDrop = "Female") {
                              groPars[[sexDrop]]$LHat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "Linf"]))
                              groPars[[sexDrop]]$kHat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "k"]))
                              groPars[[sexDrop]]$t0Hat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "t0"]))
                              groPars[[sexDrop]]$taus <<- as.matrix(sampMcmc[[sexDrop]][, grep("tau", varnames(sampMcmc[[sexDrop]]))])
                              groPars[[sexDrop]]$sigma2s <<- 1 / groPars[[sexDrop]]$taus
                              groPars[[sexDrop]]$sigma2Hat <<- apply(groPars[[sexDrop]]$sigma2s, 2, mean)
                            },
                            getMCage = function(sexDrop = "Female") {
                              tt <- as.POSIXlt(chron(spreDist[[sexDrop]]$UTC))$yday / 366
                              
                              zHat <- apply(cbind(spreDist[[sexDrop]]$Length, tt), 1, FUN = function(x) length2age(
                                numCoh = nCoho,
                                Linf = groPars[[sexDrop]]$LHat,
                                kappa = groPars[[sexDrop]]$kHat,
                                tZero = groPars[[sexDrop]]$t0Hat,
                                lengthIn = x[1],
                                timeIn = x[2],
                                sqrtSigma = sqrt(groPars[[sexDrop]]$sigma2Hat)
                              ))
                              
                              ages.f <- zHat - 1 + tt
                              AA <- floor(ages.f)
                              
                              FGlabels <- as.numeric(as.character(spreDist[[sexDrop]]$NumFG))
                              FGnames <- unique(FGlabels)
                              FG <- numeric(length(FGlabels))
                              for (FGname in 1:length(FGnames)) {
                                idx_FG <- which(FGlabels == FGnames[FGname])
                                FG[idx_FG] <- rep(FGname, length(idx_FG))
                              }
                              
                              mix_out <- data.frame(
                                Length = spreDist[[sexDrop]]$Length,
                                Date = spreDist[[sexDrop]]$UTC,
                                Day = tt,
                                Age = AA,
                                AgeNF = ages.f,
                                FG = FGlabels
                              )
                              
                              mix_out$Year <- years(mix_out$Date)
                              mix_out$Month <- as.numeric(months(as.chron(mix_out$Date)))
                              mix_out$MonthChar <- spreDist[[sexDrop]]$Month
                              mix_out$Quarter <- as.numeric(quarters(mix_out$Date))
                              mix_out$Birth <- as.numeric(as.character(mix_out$Year)) - mix_out$Age
                              
                              zeroedMonth <- ifelse(nchar(mix_out$Month) == 2, mix_out$Month, paste("0", mix_out$Month, sep = ""))
                              mix_out$CatcDate <- factor(paste(mix_out$Year,
                                                               zeroedMonth,
                                                               sep = "-"
                              ),
                              levels = paste(rep(sort(unique(mix_out$Year)), each = 12),
                                             ifelse(nchar(1:12) == 2, 1:12, paste("0", 1:12, sep = "")),
                                             sep = "-"
                              )
                              )
                              
                              groMixout[[sexDrop]] <<- mix_out
                            },
                            setMCplot = function(sexDrop = "Female", selCurve = "von Bertalanffy") {
                              nIter <- nrow(as.matrix(sampMcmc[[sexDrop]][[1]]))
                              dfLinf <- data.frame(
                                Parameter = "Linf",
                                Iter = 1:nIter,
                                Chain = as.matrix(sampMcmc[[sexDrop]][, "Linf"], chains = TRUE)[, 1],
                                Value = as.matrix(sampMcmc[[sexDrop]][, "Linf"], chains = TRUE)[, 2]
                              )
                              dfKapp <- data.frame(
                                Parameter = "Kappa",
                                Iter = 1:nIter,
                                Chain = as.matrix(sampMcmc[[sexDrop]][, "k"], chains = TRUE)[, 1],
                                Value = as.matrix(sampMcmc[[sexDrop]][, "k"], chains = TRUE)[, 2]
                              )
                              
                              ggdataSamps <- rbind(dfLinf, dfKapp)
                              ggdataSampScat <- cbind(dfLinf[, 2:3],
                                                      Linf = dfLinf[, 4],
                                                      Kappa = dfKapp[, 4]
                              )
                              
                              outPalette <- rainbow(nCoho)
                              
                              cat("\n\tSetting mcmc diagnostic plots... ", sep = "")
                              ### MCMC chain Traceplot
                              sprePlot[[sexDrop]][["traceChain"]] <<- set_ggChainTrace(ggdataSamps)
                              ### MCMC chain scatterplot
                              sprePlot[[sexDrop]][["scatLK"]] <<- set_ggChainScatter(gg_DFscat = ggdataSampScat, meanL = groPars[[sexDrop]]$LHat, meanK = groPars[[sexDrop]]$kHat)
                              ### MCMC chain Boxplot Tau
                              sprePlot[[sexDrop]][["cohoPreciGG"]] <<- set_ggTausBox(df_taus = groPars[[sexDrop]]$taus[, 1:(max(groMixout[[sexDrop]]$Age) + 1)], tauPalette = outPalette, numCoho = nCoho)
                              ### MCMC Boxplot Sigma
                              sprePlot[[sexDrop]][["cohoVariGG"]] <<- set_ggSigmaBox(df_sigma = groPars[[sexDrop]]$sigma2s[, 1:(max(groMixout[[sexDrop]]$Age) + 1)], sigPalette = outPalette, numCoho = nCoho)
                              cat("Done!", sep = "")
                              
                              coho_AL <- ddply(groMixout[[sexDrop]], .(Age), summarise,
                                               coh.mean = mean(Length), coh.var = var(Length), coh.num = length(Length)
                              )
                              
                              cat("\n\tSetting Age-Length plots... ", sep = "")
                              ### MCMC Plot Age-Length
                              sprePlot[[sexDrop]][["ageLength"]] <<- set_ggAgeLength(df_mix = groMixout[[sexDrop]], mixPalette = outPalette)
                              ### MCMC Age-Length Key
                              sprePlot[[sexDrop]][["ageLengthTbl"]] <<- set_tblAgeLength(df_mix = groMixout[[sexDrop]])
                              ### MCMC output cohort stats
                              sprePlot[[sexDrop]][["cohoStatTbl"]] <<- set_tblCohoStat(df_coho = coho_AL)
                              cat("Done!", sep = "")
                              
                              growPath <- data.frame(
                                Birth = rep(min(groMixout[[sexDrop]]$Birth):(min(groMixout[[sexDrop]]$Birth) + 11), each = length(levels(groMixout[[sexDrop]]$CatcDate))),
                                Date = rep(levels(groMixout[[sexDrop]]$CatcDate), times = length(min(groMixout[[sexDrop]]$Birth):(min(groMixout[[sexDrop]]$Birth) + 11))),
                                Length = NA
                              )
                              growPath$Age <- as.numeric(strtrim(growPath$Date, 4)) - growPath$Birth + as.numeric(substr(growPath$Date, 6, 7)) / 12
                              
                              if (selCurve == "von Bertalanffy") {
                                growPath$Length <- calcVonBert(groPars[[sexDrop]]$LHat, groPars[[sexDrop]]$kHat, growPath$Age)
                              } else {
                                growPath$Length <- calcGomp(groPars[[sexDrop]]$LHat, groPars[[sexDrop]]$kHat, growPath$Age)
                              }
                              
                              growPath$Date <- factor(growPath$Date, levels = levels(groMixout[[sexDrop]]$CatcDate))
                              growPath <- growPath[growPath$Length > floor(min(groMixout[[sexDrop]]$Length)), ]
                              
                              cat("\n\tSetting Population plots... ", sep = "")
                              ### MCMC quarter vertical hist
                              sprePlot[[sexDrop]][["histBirth"]] <<- set_ggHistBirth(df_mix = groMixout[[sexDrop]], df_grow = growPath)
                              
                              ### MCMC calc birth
                              out_birth <- table(paste(groMixout[[sexDrop]]$Year, groMixout[[sexDrop]]$Quarter, sep = "_"), groMixout[[sexDrop]]$Birth)
                              birth_melt <- melt(out_birth)
                              names(birth_melt) <- c("Catch", "Birth", "Qty")
                              birth_melt$Catch <- factor(birth_melt$Catch, levels = paste(rep(levels(groMixout[[sexDrop]]$Year), each = 4),
                                                                                          rep(1:4, times = length(levels(groMixout[[sexDrop]]$Year))),
                                                                                          sep = "_"
                              ))
                              birth_melt$Birth <- as.factor(birth_melt$Birth)
                              birth_melt <- birth_melt[birth_melt$Qty != 0, ]
                              
                              ### MCMC Catch * Quarters
                              sprePlot[[sexDrop]][["lineCatch"]] <<- set_ggCatchLine(df_birth = birth_melt)
                              
                              ### MCMC Calc Survivors
                              tot_count <- apply(out_birth, 2, sum)
                              surv_tbl <- out_birth
                              for (i in 1:nrow(out_birth)) {
                                surv_tbl[i, ] <- tot_count
                                tot_count <- tot_count - out_birth[i, ]
                              }
                              
                              surv_melt <- melt(surv_tbl)
                              names(surv_melt) <- c("Catch", "Birth", "Qty")
                              surv_melt$Catch <- factor(surv_melt$Catch, levels = paste(rep(levels(groMixout[[sexDrop]]$Year), each = 4),
                                                                                        rep(1:4, times = length(levels(groMixout[[sexDrop]]$Year))),
                                                                                        sep = "_"
                              ))
                              surv_melt <- surv_melt[!duplicated(surv_melt[, 2:3], fromLast = TRUE), ]
                              surv_melt <- surv_melt[surv_melt$Qty != 0, ]
                              surv_melt$Age <- as.numeric(strtrim(surv_melt$Catch, 4)) - surv_melt$Birth + as.numeric(substr(surv_melt$Catch, 6, 7)) / 4
                              surv_melt$Birth <- as.factor(surv_melt$Birth)
                              surv_melt$QtyNorm <- 100 * round(as.numeric(surv_melt$Qty / apply(surv_tbl, 2, max)[surv_melt$Birth]), 2)
                              # surv_melt$QtyNorm <- 100*round(as.numeric(surv_melt$Qty/max(surv_tbl)), 1)
                              
                              surv_melt$Zeta <- 0
                              for (i in unique(surv_melt$Birth)) {
                                tmp_surv_i <- surv_melt[surv_melt$Birth == i, ]
                                surv_melt$Zeta[surv_melt$Birth == i] <- c(0, -diff(tmp_surv_i$Qty) / diff(tmp_surv_i$Age) / tmp_surv_i$Qty[1])
                                # surv_melt$Zeta[surv_melt$Birth == i] <- c(0,1/diff(tmp_surv_i$Age)*log(tmp_surv_i$Qty[-nrow(tmp_surv_i)]/tmp_surv_i$Qty[-1]))
                                # surv_melt$Zeta[surv_melt$Birth == i] <- c(-diff(tmp_surv_i$Qty)/tmp_surv_i$Qty[-nrow(tmp_surv_i)], 0)
                              }
                              # surv_melt$Zeta <- 0.2*(surv_melt$Zeta)/(1/surv_melt$Zeta)
                              
                              ### MCMC Survivors * quarter
                              sprePlot[[sexDrop]][["lineSurv"]] <<- set_ggSurvLine(df_surv = surv_melt)
                              cat("Done!\n", sep = "")
                            },
                            calcMixDate = function(nAdap = 100, nSamp = 2000, nIter = 500, sexDrop = "Female", curveSel = "von Bertalanffy") {
                              cat("\n\tGetting mcmc samples... ", sep = "")
                              getMCsamps(numAdap = nAdap, numSamp = nSamp, numIter = nIter, sexDrop = sexDrop, curveSel = curveSel)
                              cat("Done!", sep = "")
                              cat("\n\tGetting growth parameters... ", sep = "")
                              getGrowPar(sexDrop = sexDrop)
                              cat("Done!", sep = "")
                              cat("\n\tGetting age estimates... ", sep = "")
                              getMCage(sexDrop = sexDrop)
                              cat("Done!", sep = "")
                              setMCplot(sexDrop = sexDrop, selCurve = curveSel)
                            },
                            ggplotMcmcOut = function(selCompo = "MCMC", selSex = "Female") {
                              switch(selCompo,
                                     MCMC = suppressWarnings(grid.arrange(sprePlot[[selSex]][["traceChain"]],
                                                                          sprePlot[[selSex]][["scatLK"]],
                                                                          sprePlot[[selSex]][["cohoPreciGG"]],
                                                                          sprePlot[[selSex]][["cohoVariGG"]],
                                                                          layout_matrix = rbind(
                                                                            c(1, 1, 1, 2),
                                                                            c(1, 1, 1, 2),
                                                                            c(4, 4, 5, 5)
                                                                          )
                                     )),
                                     Key = suppressWarnings(grid.arrange(sprePlot[[selSex]][["ageLength"]],
                                                                         sprePlot[[selSex]][["ageLengthTbl"]],
                                                                         sprePlot[[selSex]][["cohoStatTbl"]],
                                                                         layout_matrix = rbind(
                                                                           c(1, 1, 2),
                                                                           c(1, 1, 2),
                                                                           c(1, 1, 3)
                                                                         )
                                     )),
                                     Birth = suppressWarnings(grid.arrange(sprePlot[[selSex]][["histBirth"]],
                                                                           sprePlot[[selSex]][["lineCatch"]],
                                                                           sprePlot[[selSex]][["lineSurv"]],
                                                                           layout_matrix = rbind(
                                                                             c(1, 1),
                                                                             c(1, 1),
                                                                             c(2, 3)
                                                                           )
                                     ))
                              )
                            }
                          )
)



#### FisheryBySpecie#################################################
#' FisheryBySpecie
#'
#' The \code{FisheryBySpecie} class implements the class of SMART to
#'  handle species samplings.
#'
#' @docType class
#' @usage NULL
#' @keywords data
#' @return Object of \code{\link{R6Class}} with attributes and methods for the fishery data.
#'
#' @format \code{\link{R6Class}} object.
#'
#' @field species Name of the species.
#' @field year Years of the time serie.
#' @field rawLFD data.frame, raw length frequency distribution.
#' @field abuAvg data.frame, average abundances by depth' stratum.
#' @field meditsIndex data.frame, medits index by depth' stratum.
#' @field lengClas numeric, length classes.
#' @field nCoho numeric, number of cohorts.
#' @field spreDist list of DF, lfd by sex.
#' @field sprePlot plots of LFD statistics.
#' @field spreSpat list of DF, spatial distribution by sex.
#' @field sampMcmc list, mcmc output chains.
#' @field groMixout list of DF, aged individuals by sex.
#' @field groPars list of DF, growth parameters by sex.
#' @field LWpar list of DF, length/weight parameters by sex.
#'
#'   #' @section Methods:
#' \describe{
#'   \item{\code{initialize(sing_spe)}}{Automatic initialization made by the
#'   SmartProject class}
#'   \item{\code{setRawData(raw_data)}}{This method is used load the initial
#'   raw dataset}
#'   \item{\code{setYears()}}{This method is used to store the years in the
#'   provided time-serie}
#'   \item{\code{setSpecie()}}{This method is used to store the name of
#'   the species of the initial raw data}
#'   \item{\code{setLClass()}}{This method is used to store the unique
#'   length values of the sampled species}
#'   \item{\code{setDepth(bathyMatrix)}}{This method is used to assign the
#'   depth value corresponding to each sampling location}
#'   \item{\code{setStratum(vecStrata)}}{This method is used to set the
#'   depth strata of each sampling location}
#'   \item{\code{setIndSpe()}}{This method is used to aggregate the abundance
#'   data into the medits index}
#'   \item{\code{setAbuAvg()}}{This method is used to standardize the
#'   spatial abundances by depth strata}
#'   \item{\code{setNCoho(num_coh)}}{This method is used to setup the number
#'   of cohorts for the ageing module}
#'   \item{\code{setLWpar(alphaVal, betaVal, sex)}}{This method is used to
#'   store the alpha and beta values for the length/weight relationship}
#'   \item{\code{setWeight(sexVal = "Female")}}{This method is used to
#'   compute the fish weight given their length and the LWrelationship}
#'   \item{\code{setSpreDistSing()}}{This method is used to spread the
#'   aggregated LFD abundances into single individuals}
#'   \item{\code{setSprePlot(sampSex)}}{This method is used to setup the
#'   plots of the LFD statistics}
#'   \item{\code{setSpatDistSing()}}{This method is used to setup the spatial
#'   distribution of the single specimens}
#'   \item{\code{setSpatPlot(sampSex)}}{This method is used to store the
#'   spatial plots of the population}
#'   \item{\code{getMCsamps(numSamp, numAdap, numIter, sexDrop, curveSel)}}{
#'   This method is used to get a sample of the population to feed the mcmc
#'   module}
#'   \item{\code{getGrowPar(sexDrop)}}{This method is used to extract the
#'   growth parameters from the mcmc results}
#'   \item{\code{getMCage(sexDrop)}}{This method is used to assign an age to
#'   each fish}
#'   \item{\code{setMCplot(sexDrop, selCurve)}}{This method is used to setup
#'   the plot of the mcmc results}
#'   \item{\code{calcMixDate(nAdap, nSamp, nIter, sexDrop, curveSel)}}{
#'   This method is used to estimate the growth parameters of a population}
#'   \item{\code{ggplotMcmcOut(selCompo, selSex)}}{This method is used to
#'   output the stored plots of mcmc results}
#'   }

FisheryBySpecie <- R6Class("FisheryBySpecie",
                           portable = FALSE,
                           class = TRUE,
                           public = list(
                             species = NULL,
                             year = NULL,
                             rawLFD = NULL,
                             lengClas = NULL,
                             LFDPop = NULL,
                             mixPar = NULL,
                             nCoho = NULL,
                             speSex = NULL,
                             spreDist = list(),
                             spreSpat = list(),
                             sprePlot = list(),
                             sampMcmc = list(),
                             groMixout = list(),
                             groPars = list(),
                             Coh_A = NULL,
                             Coh_A_Int = NULL,
                             LWpar = list(),
                             LWstat = NULL,
                             LWstatQ = NULL,
                             selPar = NULL,
                             initialize = function(sing_spe) {
                               setRawData(sing_spe)
                               setYears()
                               setSpecie()
                               setLClass()
                             },
                             setRawData = function(raw_data) {
                               rawLFD <<- raw_data
                             },
                             plotLFD = function() {
                               plotSpeAllYea(rawLFD)
                             },
                             setYears = function() {
                               year <<- sort(unique(years(rawLFD[, "Date"])), decreasing = FALSE)
                             },
                             setSpecie = function() {
                               species <<- unique(rawLFD[, "Species"])
                             },
                             setLClass = function() {
                               lengClas <<- seq(from = min(rawLFD[, "Class"]), to = max(rawLFD[, "Class"]), by = 1)
                             },
                             setNCoho = function(num_coh) {
                               nCoho <<- num_coh
                             },
                             setLWpar = function(alphaVal, betaVal, sex) {
                               LWpar[[sex]] <<- list(alpha = as.numeric(alphaVal), beta = as.numeric(betaVal))
                             },
                             setAvailSex = function() {
                               speSex <<- sort(names(which(lapply(spreDist, nrow) > 0)))
                             },
                             setSprePlot = function(sampSex) {
                               sprePlot[[sampSex]] <<- list(
                                 histLfdTot = set_ggHistLfdTot(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                 histUtcTot = set_ggHistUtcTot(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                 dotUtcSplit = set_ggDotUtcSplit(spreDist[[sampSex]]) + scale_color_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE"))),
                                 histUtcLfd = set_ggHistUtcLfd(spreDist[[sampSex]]) + scale_fill_manual(values = ifelse(sampSex == "Female", "#FF6A6A", ifelse(sampSex == "Male", "#63B8FF", "#63FFAE")))
                               )
                             },
                             setSpreDistSing = function() {
                               for (sex in c("Female", "Male", "Unsex")) {
                                 tmp_spre <- rawLFD[!is.na(rawLFD$numFG), c("Date", "Class", "numFG", sex)]
                                 
                                 num_sex <- sum(tmp_spre[, 4])
                                 cat("Found", num_sex, sex, as.character(species), "samples\n", sep = " ")
                                 
                                 spreDist <- data.frame(
                                   UTC = rep(tmp_spre$Date, tmp_spre[, 4]),
                                   Length = rep(tmp_spre$Class, tmp_spre[, 4]) + runif(num_sex, -0.5, 0.5),
                                   NumFG = rep(tmp_spre$numFG, tmp_spre[, 4])
                                 )
                                 
                                 spreDist$Year <- years(spreDist$UTC)
                                 spreDist$Month <- months(as.chron(spreDist$UTC))
                                 
                                 spreDist[[sex]] <<- spreDist
                                 setSprePlot(sampSex = sex)
                               }
                             },
                             setSpatDistSing = function() {
                               for (sex in c("Female", "Male", "Unsex")) {
                                 tmp_fishSpat <- rawLFD[!is.na(rawLFD$numFG) & rawLFD[, sex] > 0, c("Lon", "Lat", "numFG", sex)]
                                 if (nrow(tmp_fishSpat) > 0) {
                                   barploFgAll <- data.frame(table(tmp_fishSpat$numFG))
                                   barploFgAll <- barploFgAll[order(as.numeric(as.character(barploFgAll[, 1]))), ]
                                   barploFgAll$FG <- factor(barploFgAll$Var1, levels = barploFgAll$Var1)
                                   barploFgAll$relFreq <- round(100 * barploFgAll$Freq / sum(barploFgAll$Freq), 1)
                                   spreSpat[[sex]] <<- barploFgAll
                                   setSpatPlot(sampSex = sex)
                                 }
                               }
                             },
                             setSpatPlot = function(sampSex) {
                               sprePlot[[sampSex]][["spatAbbTbl"]] <<- set_spatAbbTbl(spreSpat[[sampSex]])
                               sprePlot[[sampSex]][["spatAbsFreq"]] <<- set_spatAbsFreq(spreSpat[[sampSex]])
                               sprePlot[[sampSex]][["spatRelFreq"]] <<- set_spatRelFreq(spreSpat[[sampSex]])
                             },
                             # setPreMix = function(){
                             #   # preMix <<- weight2number(rawLFD[,c("DATE","Class","Unsex")])
                             # },
                             setWeight = function(sexVal = "Female") {
                               groMixout[[sexVal]]$Weight <<- LWpar[[sexVal]][["alpha"]] * groMixout[[sexVal]]$Length^LWpar[[sexVal]][["beta"]]
                             },
                             setLWstat = function(lwUnit = "Length") {
                               if (length(groMixout) > 1) {
                                 tmp_out <- do.call(rbind, groMixout)
                               } else {
                                 tmp_out <- groMixout[[1]]
                               }
                               tmp_out$Weight <- ceiling(tmp_out$Weight)
                               tmp_out$Season <- factor(quarters(tmp_out$Date + 30), levels = c("1Q", "2Q", "3Q", "4Q"), labels = c("winter", "spring", "summer", "fall"))
                               
                               if (lwUnit == "Length") {
                                 tmp_lw <- ddply(tmp_out, .(Weight, FG), summarise,
                                                 avgLen = mean(Length), sdLen = sd(Length), relAbb = length(Length)
                                 )
                                 tmp_lw <- tmp_lw[-which(is.na(tmp_lw), arr.ind = TRUE)[, 1], ]
                               } else {
                                 tmp_lw <- as.data.frame(table(Weight = tmp_out$Weight, FG = tmp_out$FG), stringsAsFactors = FALSE)
                                 tmp_lw$Weight <- as.numeric(tmp_lw$Weight)
                                 tmp_lw <- tmp_lw[tmp_lw$Freq > 0, ]
                               }
                               LWstat <<- tmp_lw
                               if (lwUnit == "Length") {
                                 tmp_lwq <- ddply(tmp_out, .(Weight, FG, Season), summarise,
                                                  avgLen = mean(Length), sdLen = sd(Length), relAbb = length(Length)
                                 )
                                 tmp_lwq <- tmp_lwq[-which(is.na(tmp_lwq), arr.ind = TRUE)[, 1], ]
                               } else {
                                 tmp_lwq <- ddply(tmp_out, .(Weight, FG, Season), summarise, Freq = length(Weight))
                                 tmp_lwq <- tmp_lwq[tmp_lwq$Freq > 0, ]
                               }
                               LWstatQ <<- tmp_lwq
                             },
                             getMCsamps = function(numSamp = 2000, numAdap = 100, numIter = 500, sexDrop = "Female", curveSel = "von Bertalanffy") {
                               sub_idx <- sample(1:nrow(spreDist[[sexDrop]]), size = numSamp, replace = ifelse(numSamp < nrow(spreDist[[sexDrop]]), FALSE, TRUE))
                               sub_data <- spreDist[[sexDrop]][sub_idx, ]
                               
                               N <- length(sub_data$Length)
                               alpha <- rep(1, nCoho)
                               Z <- rep(NA, N)
                               Z[which.min(sub_data$Length)] <- 1
                               Z[which.max(sub_data$Length)] <- nCoho
                               
                               dataList <- list(
                                 y = sub_data$Length,
                                 maxLeng = max(sub_data$Length), ## !!!
                                 alpha = alpha,
                                 Z = Z,
                                 N = N,
                                 Nclust = nCoho
                               )
                               
                               inits <- list(
                                 list(Linf = min(sub_data$Length), k = 0.5, t0 = 0.0),
                                 list(Linf = mean(sub_data$Length), k = 0.5, t0 = 0.0),
                                 list(Linf = max(sub_data$Length), k = 0.5, t0 = 0.0)
                               )
                               
                               modelGrow <- ifelse(curveSel == "von Bertalanffy",
                                                   system.file("model/bertGrow.jags", package = "smartR"),
                                                   system.file("model/gompGrow.jags", package = "smartR")
                               )
                               
                               jags.m <- jags.model(modelGrow,
                                                    data = dataList,
                                                    inits = inits,
                                                    n.chains = 3,
                                                    n.adapt = numAdap
                               )
                               
                               ### MCMC chain sampling
                               # n.iter <- 500
                               obsNode <- c("Linf", "k", "t0", "tau", "p")
                               samps <- coda.samples(jags.m, obsNode, n.iter = numIter)
                               
                               sampMcmc[[sexDrop]] <<- samps
                             },
                             getGrowPar = function(sexDrop = "Female") {
                               groPars[[sexDrop]]$LHat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "Linf"]))
                               groPars[[sexDrop]]$kHat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "k"]))
                               groPars[[sexDrop]]$t0Hat <<- mean(as.matrix(sampMcmc[[sexDrop]][, "t0"]))
                               groPars[[sexDrop]]$taus <<- as.matrix(sampMcmc[[sexDrop]][, grep("tau", varnames(sampMcmc[[sexDrop]]))])
                               groPars[[sexDrop]]$sigma2s <<- 1 / groPars[[sexDrop]]$taus
                               groPars[[sexDrop]]$sigma2Hat <<- apply(groPars[[sexDrop]]$sigma2s, 2, mean)
                             },
                             getMCage = function(sexDrop = "Female") {
                               tt <- as.POSIXlt(chron(spreDist[[sexDrop]]$UTC))$yday / 366
                               
                               zHat <- apply(cbind(spreDist[[sexDrop]]$Length, tt), 1, FUN = function(x) length2age(
                                 numCoh = nCoho,
                                 Linf = groPars[[sexDrop]]$LHat,
                                 kappa = groPars[[sexDrop]]$kHat,
                                 tZero = groPars[[sexDrop]]$t0Hat,
                                 lengthIn = x[1],
                                 timeIn = x[2],
                                 sqrtSigma = sqrt(groPars[[sexDrop]]$sigma2Hat)
                               ))
                               
                               ages.f <- zHat - 1 + tt
                               AA <- floor(ages.f)
                               
                               FGlabels <- as.numeric(as.character(spreDist[[sexDrop]]$NumFG))
                               FGnames <- unique(FGlabels)
                               FG <- numeric(length(FGlabels))
                               for (FGname in 1:length(FGnames)) {
                                 idx_FG <- which(FGlabels == FGnames[FGname])
                                 FG[idx_FG] <- rep(FGname, length(idx_FG))
                               }
                               
                               mix_out <- data.frame(
                                 Length = spreDist[[sexDrop]]$Length,
                                 Date = spreDist[[sexDrop]]$UTC,
                                 Day = tt,
                                 Age = AA,
                                 AgeNF = ages.f,
                                 FG = FGlabels
                               )
                               
                               mix_out$Year <- years(mix_out$Date)
                               mix_out$Month <- as.numeric(months(as.chron(mix_out$Date)))
                               mix_out$MonthChar <- spreDist[[sexDrop]]$Month
                               mix_out$Quarter <- as.numeric(quarters(mix_out$Date))
                               mix_out$Birth <- as.numeric(as.character(mix_out$Year)) - mix_out$Age
                               
                               zeroedMonth <- ifelse(nchar(mix_out$Month) == 2, mix_out$Month, paste("0", mix_out$Month, sep = ""))
                               mix_out$CatcDate <- factor(paste(mix_out$Year,
                                                                zeroedMonth,
                                                                sep = "-"
                               ),
                               levels = paste(rep(sort(unique(mix_out$Year)), each = 12),
                                              ifelse(nchar(1:12) == 2, 1:12, paste("0", 1:12, sep = "")),
                                              sep = "-"
                               )
                               )
                               
                               groMixout[[sexDrop]] <<- mix_out
                             },
                             setMCplot = function(sexDrop = "Female", selCurve = "von Bertalanffy") {
                               nIter <- nrow(as.matrix(sampMcmc[[sexDrop]][[1]]))
                               dfLinf <- data.frame(
                                 Parameter = "Linf",
                                 Iter = 1:nIter,
                                 Chain = as.matrix(sampMcmc[[sexDrop]][, "Linf"], chains = TRUE)[, 1],
                                 Value = as.matrix(sampMcmc[[sexDrop]][, "Linf"], chains = TRUE)[, 2]
                               )
                               dfKapp <- data.frame(
                                 Parameter = "Kappa",
                                 Iter = 1:nIter,
                                 Chain = as.matrix(sampMcmc[[sexDrop]][, "k"], chains = TRUE)[, 1],
                                 Value = as.matrix(sampMcmc[[sexDrop]][, "k"], chains = TRUE)[, 2]
                               )
                               
                               ggdataSamps <- rbind(dfLinf, dfKapp)
                               ggdataSampScat <- cbind(dfLinf[, 2:3],
                                                       Linf = dfLinf[, 4],
                                                       Kappa = dfKapp[, 4]
                               )
                               
                               outPalette <- rainbow(nCoho)
                               
                               cat("\n\tSetting mcmc diagnostic plots... ", sep = "")
                               ### MCMC chain Traceplot
                               sprePlot[[sexDrop]][["traceChain"]] <<- set_ggChainTrace(ggdataSamps)
                               ### MCMC chain scatterplot
                               sprePlot[[sexDrop]][["scatLK"]] <<- set_ggChainScatter(gg_DFscat = ggdataSampScat, meanL = groPars[[sexDrop]]$LHat, meanK = groPars[[sexDrop]]$kHat)
                               ### MCMC chain Boxplot Tau
                               sprePlot[[sexDrop]][["cohoPreciGG"]] <<- set_ggTausBox(df_taus = groPars[[sexDrop]]$taus[, 1:(max(groMixout[[sexDrop]]$Age) + 1)], tauPalette = outPalette, numCoho = nCoho)
                               ### MCMC Boxplot Sigma
                               sprePlot[[sexDrop]][["cohoVariGG"]] <<- set_ggSigmaBox(df_sigma = groPars[[sexDrop]]$sigma2s[, 1:(max(groMixout[[sexDrop]]$Age) + 1)], sigPalette = outPalette, numCoho = nCoho)
                               cat("Done!", sep = "")
                               
                               coho_AL <- ddply(groMixout[[sexDrop]], .(Age), summarise,
                                                coh.mean = mean(Length), coh.var = var(Length), coh.num = length(Length)
                               )
                               
                               cat("\n\tSetting Age-Length plots... ", sep = "")
                               ### MCMC Plot Age-Length
                               sprePlot[[sexDrop]][["ageLength"]] <<- set_ggAgeLength(df_mix = groMixout[[sexDrop]], mixPalette = outPalette)
                               ### MCMC Age-Length Key
                               sprePlot[[sexDrop]][["ageLengthTbl"]] <<- set_tblAgeLength(df_mix = groMixout[[sexDrop]])
                               ### MCMC output cohort stats
                               sprePlot[[sexDrop]][["cohoStatTbl"]] <<- set_tblCohoStat(df_coho = coho_AL)
                               cat("Done!", sep = "")
                               
                               growPath <- data.frame(
                                 Birth = rep(min(groMixout[[sexDrop]]$Birth):(min(groMixout[[sexDrop]]$Birth) + 11), each = length(levels(groMixout[[sexDrop]]$CatcDate))),
                                 Date = rep(levels(groMixout[[sexDrop]]$CatcDate), times = length(min(groMixout[[sexDrop]]$Birth):(min(groMixout[[sexDrop]]$Birth) + 11))),
                                 Length = NA
                               )
                               growPath$Age <- as.numeric(strtrim(growPath$Date, 4)) - growPath$Birth + as.numeric(substr(growPath$Date, 6, 7)) / 12
                               
                               if (selCurve == "von Bertalanffy") {
                                 growPath$Length <- calcVonBert(groPars[[sexDrop]]$LHat, groPars[[sexDrop]]$kHat, growPath$Age)
                               } else {
                                 growPath$Length <- calcGomp(groPars[[sexDrop]]$LHat, groPars[[sexDrop]]$kHat, growPath$Age)
                               }
                               
                               growPath$Date <- factor(growPath$Date, levels = levels(groMixout[[sexDrop]]$CatcDate))
                               growPath <- growPath[growPath$Length > floor(min(groMixout[[sexDrop]]$Length)), ]
                               
                               cat("\n\tSetting Population plots... ", sep = "")
                               ### MCMC quarter vertical hist
                               sprePlot[[sexDrop]][["histBirth"]] <<- set_ggHistBirth(df_mix = groMixout[[sexDrop]], df_grow = growPath)
                               
                               ### MCMC calc birth
                               out_birth <- table(paste(groMixout[[sexDrop]]$Year, groMixout[[sexDrop]]$Quarter, sep = "_"), groMixout[[sexDrop]]$Birth)
                               birth_melt <- melt(out_birth)
                               names(birth_melt) <- c("Catch", "Birth", "Qty")
                               birth_melt$Catch <- factor(birth_melt$Catch, levels = paste(rep(levels(groMixout[[sexDrop]]$Year), each = 4),
                                                                                           rep(1:4, times = length(levels(groMixout[[sexDrop]]$Year))),
                                                                                           sep = "_"
                               ))
                               birth_melt$Birth <- as.factor(birth_melt$Birth)
                               birth_melt <- birth_melt[birth_melt$Qty != 0, ]
                               
                               ### MCMC Catch * Quarters
                               sprePlot[[sexDrop]][["lineCatch"]] <<- set_ggCatchLine(df_birth = birth_melt)
                               
                               ### MCMC Calc Survivors
                               tot_count <- apply(out_birth, 2, sum)
                               surv_tbl <- out_birth
                               for (i in 1:nrow(out_birth)) {
                                 surv_tbl[i, ] <- tot_count
                                 tot_count <- tot_count - out_birth[i, ]
                               }
                               
                               surv_melt <- melt(surv_tbl)
                               names(surv_melt) <- c("Catch", "Birth", "Qty")
                               surv_melt$Catch <- factor(surv_melt$Catch, levels = paste(rep(levels(groMixout[[sexDrop]]$Year), each = 4),
                                                                                         rep(1:4, times = length(levels(groMixout[[sexDrop]]$Year))),
                                                                                         sep = "_"
                               ))
                               surv_melt <- surv_melt[!duplicated(surv_melt[, 2:3], fromLast = TRUE), ]
                               surv_melt <- surv_melt[surv_melt$Qty != 0, ]
                               surv_melt$Age <- as.numeric(strtrim(surv_melt$Catch, 4)) - surv_melt$Birth + as.numeric(substr(surv_melt$Catch, 6, 7)) / 4
                               surv_melt$Birth <- as.factor(surv_melt$Birth)
                               surv_melt$QtyNorm <- 100 * round(as.numeric(surv_melt$Qty / apply(surv_tbl, 2, max)[surv_melt$Birth]), 2)
                               # surv_melt$QtyNorm <- 100*round(as.numeric(surv_melt$Qty/max(surv_tbl)), 1)
                               
                               surv_melt$Zeta <- 0
                               for (i in unique(surv_melt$Birth)) {
                                 tmp_surv_i <- surv_melt[surv_melt$Birth == i, ]
                                 surv_melt$Zeta[surv_melt$Birth == i] <- c(0, -diff(tmp_surv_i$Qty) / diff(tmp_surv_i$Age) / tmp_surv_i$Qty[1])
                                 # surv_melt$Zeta[surv_melt$Birth == i] <- c(0,1/diff(tmp_surv_i$Age)*log(tmp_surv_i$Qty[-nrow(tmp_surv_i)]/tmp_surv_i$Qty[-1]))
                                 # surv_melt$Zeta[surv_melt$Birth == i] <- c(-diff(tmp_surv_i$Qty)/tmp_surv_i$Qty[-nrow(tmp_surv_i)], 0)
                               }
                               # surv_melt$Zeta <- 0.2*(surv_melt$Zeta)/(1/surv_melt$Zeta)
                               
                               ### MCMC Survivors * quarter
                               sprePlot[[sexDrop]][["lineSurv"]] <<- set_ggSurvLine(df_surv = surv_melt)
                               cat("Done!\n", sep = "")
                             },
                             calcMixDate = function(nAdap = 100, nSamp = 2000, nIter = 500, sexDrop = "Female", curveSel = "von Bertalanffy") {
                               cat("\n\tGetting mcmc samples... ", sep = "")
                               getMCsamps(numAdap = nAdap, numSamp = nSamp, numIter = nIter, sexDrop = sexDrop, curveSel = curveSel)
                               cat("Done!", sep = "")
                               cat("\n\tGetting growth parameters... ", sep = "")
                               getGrowPar(sexDrop = sexDrop)
                               cat("Done!", sep = "")
                               cat("\n\tGetting age estimates... ", sep = "")
                               getMCage(sexDrop = sexDrop)
                               cat("Done!", sep = "")
                               setMCplot(sexDrop = sexDrop, selCurve = curveSel)
                             },
                             ggplotMcmcOut = function(selCompo = "MCMC", selSex = "Female") {
                               switch(selCompo,
                                      MCMC = suppressWarnings(grid.arrange(sprePlot[[selSex]][["traceChain"]],
                                                                           sprePlot[[selSex]][["scatLK"]],
                                                                           sprePlot[[selSex]][["cohoPreciGG"]],
                                                                           sprePlot[[selSex]][["cohoVariGG"]],
                                                                           layout_matrix = rbind(
                                                                             c(1, 1, 1, 2),
                                                                             c(1, 1, 1, 2),
                                                                             c(4, 4, 5, 5)
                                                                           )
                                      )),
                                      Key = suppressWarnings(grid.arrange(sprePlot[[selSex]][["ageLength"]],
                                                                          sprePlot[[selSex]][["ageLengthTbl"]],
                                                                          sprePlot[[selSex]][["cohoStatTbl"]],
                                                                          layout_matrix = rbind(
                                                                            c(1, 1, 2),
                                                                            c(1, 1, 2),
                                                                            c(1, 1, 3)
                                                                          )
                                      )),
                                      Birth = suppressWarnings(grid.arrange(sprePlot[[selSex]][["histBirth"]],
                                                                            sprePlot[[selSex]][["lineCatch"]],
                                                                            sprePlot[[selSex]][["lineSurv"]],
                                                                            layout_matrix = rbind(
                                                                              c(1, 1),
                                                                              c(1, 1),
                                                                              c(2, 3)
                                                                            )
                                      ))
                               )
                             }
                           )
)



#### FishFleet################################################
#' FishFleet
#'
#' The \code{FishFleet} class implements the class of SMART
#' to manage fleet data.
#'
#' @docType class
#' @usage NULL
#' @keywords data
#' @return Object of \code{\link{R6Class}} with attributes and methods for the fishery data.
#'
#' @format \code{\link{R6Class}} object.
#'
#' @field rawRegister data.frame, raw fleet register data.
#' @field vmsRegister data.frame, raw fleet register data for vms vessels only.
#' @field rawEffort list of DF, raw effort data.
#' @field dayEffoMatr list of DF, daily aggregated effort data.
#' @field prodMatr list of DF, production data.
#' @field effoProd list of DF, merged effort and production data.
#' @field effoProdMont list of DF, monthly aggregated effort and production data.
#' @field effoMont list of DF, monthly aggregated effort data.
#' @field effoProdAll data.frame, monthly aggregated effort and production data.
#' @field effoAll data.frame, monthly aggregated effort data.
#' @field regHarbsUni data.frame, Harbours name, longitude, latitude and
#' distance from the environment grid.
#' @field regHarbsBox data.frame, Harbours name, longitude, latitude, number of
#' registered vessels and distance from the environment grid within the grid
#' box.
#' @field rawProduction list of DF, raw production data.
#' @field rawEconomy data.frame, raw economic data.
#' @field registerIds character, vessel identification from fleet register.
#' @field predProd list of matrix, simulated production.
#' @field productionIds list of int, vessel ids with production data available.
#' @field prodSpec list of character, species with production data.
#' @field specSett list of DF, logit parameter settings by species.
#' @field specLogit list, logit results by species.
#' @field effortIds list of int, vessel ids with effort data available.
#' @field idsEffoProd list of int, merged vessel ids with both effort and
#' production data available.
#' @field effoProdAllLoa data.frame, monthly aggregated effort, production and
#' loa data.
#' @field effoAllLoa  data.frame, monthly aggregated effort and loa data.
#' @field effortIndex data.frame, effort index by vessel, year and month with
#' loa data.
#' @field daysAtSea data.frame, days at sea index by vessel, year, month with
#' loa and Kw data.
#' @field prodIndex data.frame, production index by vessel, year and month.
#' @field resNNLS list, lander results by species.
#' @field betaMeltYear list of DF, melted yearly productivity by species.
#' @field prodMeltYear list of DF, melted yearly production by species.
#' @field fishPoinPara data.frame, fishing point parameters.
#' @field ecoPrice list of DF, price/size by species.
#' @field inSpatialReg data.frame, input for spatial index regression.
#' @field inEffortReg data.frame, input for effort index regression.
#' @field inProductionReg data.frame, input for production index regression.
#' @field outSpatialReg list, output for spatial index regression.
#' @field outEffortReg list, output for effort index regression.
#' @field outProductionReg list, output for production index regression.
#' @field plotSpatialReg ggplot, spatial index regression results.
#' @field plotEffortReg ggplot, effort index regression results.
#' @field plotProductionReg ggplot, production index regression results.
#'
#' @section Methods:
#' \describe{
#'   \item{\code{setVmsRegister()}}{This method is used to exclude the fleet
#'   register records of vessels without vms}
#'   \item{\code{setRegHarbs()}}{This method is used to fetch the harbours
#'   coordinates}
#'   \item{\code{setEcoPrice(sel_specie, price_df)}}{This method is used to set
#'   the price/size attribute by species}
#'   \item{\code{saveFleetHarb(harb_path)}}{This method is used to export the
#'   rds with the harbours' coordinates}
#'   \item{\code{loadFleetHarb(harb_path)}}{This method is used to import the
#'   rds with the harbours' coordinates}
#'   \item{\code{loadFleetRegis(register_path)}}{This method is used to load the
#'   raw fleet register}
#'   \item{\code{loadMatEffort(effort_path)}}{This method is used to import the
#'   raw effort matrix}
#'   \item{\code{loadRawEconomy(economic_path)}}{This method is used to load the
#'   raw csv file with economic data}
#'   \item{\code{setInSpatial()}}{This method is used to setup the input for the
#'   spatial regression}
#'   \item{\code{setInEffort()}}{This method is used to setup the input for the
#'   effprt regression}
#'   \item{\code{getRegSpatial()}}{This method is used to compute the spatial
#'   regression}
#'   \item{\code{getRegEffort()}}{This method is used to compute the effort
#'   regression}
#'   \item{\code{getRegProduction()}}{This method is used to compute the
#'   production regression}
#'   \item{\code{getCostOutput()}}{This method is a wrapper function to get the
#'   economic regressions}
#'   \item{\code{setCostPlot()}}{This method is used to setup the plot of the
#'   economic regression}
#'   \item{\code{loadProduction(production_path)}}{This method is used to load
#'   the raw csv of the production data}
#'   \item{\code{setFishPoinPara(speed_range, depth_range)}}{This method is
#'   used to setup the fishign points parameters}
#'   \item{\code{setWeekMonthNum()}}{This method is used to assign the week and
#'   month num to the raw effort data}
#'   \item{\code{setFishPoin()}}{This method is used to filter the
#'   fishing points}
#'   \item{\code{plotFishPoinStat()}}{This method is used to show the basic
#'   statistics for the fishing points}
#'   \item{\code{plotSpeedDepth(which_year, speed_range, depth_range)}}{
#'   This method is used to show the speed/depth profile}
#'   \item{\code{setEffortIds()}}{This method is used to set the distinct
#'   vessel' ids in the effort dataset}
#'   \item{\code{setProdSpec()}}{This method is used to set the distinct species
#'   in the production dataset}
#'   \item{\code{setBetaMeltYear(species)}}{This method is used to set the melted
#'   yearly productivity by species}
#'   \item{\code{setProdMeltYear(species)}}{This method is used to set the melted
#'   yearly production by species}
#'   \item{\code{plotTotProd(species)}}{This method is used to plot the total
#'   production by species}
#'   \item{\code{plotNNLS(species, thresR2)}}{This method is used to show the
#'   NNLS results}
#'   \item{\code{setSpecSettItm(species, thresh, brea, max_xlim)}}{
#'   This method is used to set the logit parameters by species}
#'   \item{\code{plotLogitROC(selSpecie)}}{This method is used to show the
#'   ROC of the logit results}
#'   \item{\code{setSpecLogitConf(selSpecie, cutoff)}}{This method is used to
#'   set the confusion matrix of the logit results by species}
#'   \item{\code{setLogitTrain(selSpecie, train, cp_val, cv_val)}}{
#'   This method is used to setup the train dataset for the logit model}
#'   \item{\code{setLogitTest(selSpecie, test)}}{This method is used to setup
#'   the test dataset for the logit model}
#'   \item{\code{setLogitPred(selSpecie, test)}}{This method is used to compute
#'   the prediction for the logit model}
#'   \item{\code{setLogitCut(selSpecie)}}{This method is used to tune the
#'   cutoff of the logit model}
#'   \item{\code{setLogitRoc(selSpecie)}}{This method is used to set the ROC of
#'   the logit model}
#'   \item{\code{setLogitConf(selSpecie, test)}}{This method is used to
#'   set the confusion matrix of the logit results}
#'   \item{\code{setSpecLogit(selSpecie, selModel, cp, cv)}}{This method is
#'    a wrapper function to get the logit model}
#'   \item{\code{getMatSpeLand(species)}}{This method is used to get the input
#'   data for the logit model}
#'   \item{\code{setEffoProdAll()}}{This method is used to combine the
#'   effort/production data from the yearly list into a single data.frame}
#'   \item{\code{setEffoAll()}}{This method is used to combine the
#'   effort data from the yearly list into a single data.frame}
#'   \item{\code{setEffoProdAllLoa()}}{This method is used to add the LOA data
#'   to the effort/production data}
#'   \item{\code{setEffoAllLoa()}}{This method is used to add the LOA data
#'   to the effort data}
#'   \item{\code{setProdIds()}}{This method is used to get the vessel ids with
#'   production data available}
#'   \item{\code{setIdsEffoProd()}}{This method is used to get the vessel ids
#'    with both effort and production data available}
#'   \item{\code{plotCountIDsEffoProd()}}{This method is used to set the plot of
#'   the basic statistics of the effort/production data}
#'   \item{\code{plotCountIDsEffo()}}{This method is used to set the plot of
#'   the basic statistics of the effort data}
#'   \item{\code{plotCountIDsProd()}}{This method is used to set the plot of
#'   the basic statistics of the production data}
#'   \item{\code{setEffoProdMatr()}}{This method is used to merge the effort
#'   and production data}
#'   \item{\code{setEffoProdMont()}}{This method is used to aggregate the
#'   effort/production data by month}
#'   \item{\code{setEffoMont()}}{This method is used to aggregate the
#'   effort data by month}
#'   \item{\code{setProdMatr()}}{This method is used to create the production
#'   matrix from the raw production data}
#'   \item{\code{setDayEffoMatrGround(maxFG)}}{This method is used to assign
#'   the fishing grounds to the raw effort data}
#'   \item{\code{readRegisterEU(reg_path)}}{This method is used to load the
#'   raw european fleet register}
#'   \item{\code{cleanRegister()}}{This method is used to clean the raw data in
#'   the fleet register}
#'   \item{\code{plotRegSum()}}{This method is used to plot the basic statistics
#'   for the fleet register data}
#'   \item{\code{setRegIds()}}{This method is used to get the distinct vessels
#'   ids in the fleet register}
#'   }

FishFleet <- R6Class("fishFleet",
                     portable = FALSE,
                     class = TRUE,
                     public = list(
                       rawRegister = NULL,
                       vmsRegister = NULL,
                       rawEffort = NULL,
                       weekEffoMatr = NULL,
                       dayEffoMatr = NULL,
                       prodMatr = NULL,
                       effoProd = NULL,
                       effoProdMont = NULL,
                       effoMont = NULL,
                       effoProdAll = NULL,
                       effoAll = NULL,
                       trackHarbs = NULL,
                       regHarbsUni = NULL,
                       regHarbsBox = NULL,
                       rawSelectivity = NULL,
                       rawProduction = NULL,
                       rawEconomy = NULL,
                       registerIds = NULL,
                       predProd = NULL,
                       productionIds = NULL,
                       prodIdsLoa = NULL,
                       prodSpec = NULL,
                       specSett = NULL,
                       specLogit = NULL,
                       effortIds = NULL,
                       idsEffoProd = NULL,
                       effoProdAllLoa = NULL,
                       effoAllLoa = NULL,
                       effortIndex = NULL,
                       daysAtSea = NULL,
                       prodIndex = NULL,
                       resNNLS = NULL,
                       betaAvg = NULL,
                       effortAvg = NULL,
                       betaMeltYear = NULL,
                       prodMeltYear = NULL,
                       fishPoinPara = NULL,
                       ecoPrice = NULL,
                       inSpatialReg = NULL,
                       inEffortReg = NULL,
                       inProductionReg = NULL,
                       outSpatialReg = NULL,
                       outEffortReg = NULL,
                       outProductionReg = NULL,
                       plotSpatialReg = NULL,
                       plotEffortReg = NULL,
                       plotProductionReg = NULL,
                       setVmsRegister = function() {
                         vmsRegister <<- suppressWarnings(rawRegister[as.numeric(substr(rawRegister$CFR, 4, nchar(rawRegister$CFR[1]))) %in% effortIds$All, ])
                       },
                       setRegHarbs = function() {
                         cat("\n\tGetting Harbours coordinates...\n", cat = "")
                         regHarbsUni <<- nominatim_osm(address = sort(unique(vmsRegister$Port.Name)))
                         cat("\n\t\tHarbours geocoding completed!", cat = "")
                       },
                       setEcoPrice = function(sel_specie, price_df) {
                         if (is.null(ecoPrice)) ecoPrice <<- list()
                         ecoPrice[[sel_specie]] <<- price_df
                       },
                       saveFleetHarb = function(harb_path) {
                         saveRDS(object = regHarbsUni, file = harb_path)
                       },
                       loadFleetHarb = function(harb_path) {
                         regHarbsUni <<- readRDS(file = harb_path)
                       },
                       setBetaAvg = function(sel_specie) {
                         tmp_df <- data.frame(
                           Month = resNNLS[[sel_specie]]$SceMat$MONTH,
                           resNNLS[[sel_specie]]$bmat
                         )
                         betaAvg <<- aggregate(. ~ Month, tmp_df, mean)
                       },
                       setEffortAvg = function() {
                         tmp_effo <- aggregate(. ~ Year + MonthNum, effoAllLoa[, -c(2, ncol(effoAllLoa))], sum)
                         effortAvg <<- aggregate(. ~ MonthNum, tmp_effo[, -1], mean)
                       },
                       loadFleetRegis = function(register_path) {
                         cat("\nLoading raw Fleet Register data... ", sep = "")
                         rawRegister <<- readRegisterEU(register_path)
                         cat("\nFleet Register loading completed!", sep = "")
                       },
                       loadMatEffort = function(effort_path) {
                         cat("\nLoading Effort data... ", sep = "")
                         rawEffort <<- readRDS(effort_path)
                         cat("Done!", sep = "")
                       },
                       loadRawEconomy = function(economic_path) {
                         cat("\nLoading Economic data... ", sep = "")
                         rawEconomy <<- read.csv(file = economic_path, stringsAsFactors = FALSE)
                         cat("Done!", sep = "")
                       },
                       setYearEconomy = function() {
                         rawEconomy$Year <<- as.numeric(substr(rawEconomy$DateStart, 7, nchar(rawEconomy$DateStart)))
                       },
                       setInSpatial = function() {
                         agg_EffInd <- aggregate(EffInd ~ I_NCEE + Year + Loa, effortIndex, sum)
                         tmp_spatCost <- rawEconomy[, c("VessID", "Year", "SpatialCost")]
                         inSpatialReg <<- merge(
                           x = tmp_spatCost, y = agg_EffInd,
                           by.x = c("VessID", "Year"), by.y = c("I_NCEE", "Year")
                         )
                       },
                       setInEffort = function() {
                         agg_DaysAtSea <- aggregate(Freq ~ I_NCEE + Year + Loa + Kw, daysAtSea, sum)
                         tmp_effoCost <- rawEconomy[, c("VessID", "Year", "EffortCost")]
                         inEffortReg <<- merge(
                           x = tmp_effoCost, y = agg_DaysAtSea,
                           by.x = c("VessID", "Year"), by.y = c("I_NCEE", "Year")
                         )
                       },
                       getRegSpatial = function() {
                         outSpatialReg <<- lm(formula = SpatialCost ~ EffInd + Loa - 1, data = inSpatialReg)
                       },
                       getRegEffort = function() {
                         outEffortReg <<- lm(formula = EffortCost ~ Freq + Loa + Kw - 1, data = inEffortReg)
                       },
                       getRegProduction = function() {
                         outProductionReg <<- lm(formula = ProductionCost ~ Production - 1, data = inProductionReg)
                       },
                       getCostOutput = function() {
                         getRegSpatial()
                         getRegEffort()
                         getRegProduction()
                       },
                       setCostPlot = function() {
                         plotSpatialReg <<- ggplot_spatialRegression(df_spatialIn = inSpatialReg, reg_spatialOut = outSpatialReg)
                         plotEffortReg <<- ggplot_effortRegression(df_effortIn = inEffortReg, reg_effortOut = outEffortReg)
                         plotProductionReg <<- ggplot_productionRegression(df_productionIn = inProductionReg, reg_productionOut = outProductionReg)
                       },
                       loadProduction = function(production_path) {
                         cat("\nLoading Production data... ", sep = "")
                         sort_files <- sort(production_path)
                         rawProduction <<- list()
                         for (i in 1:length(sort_files)) {
                           tmp_mat <- read.csv(sort_files[i], stringsAsFactors = FALSE)
                           numYea <- as.character(sort(unique(years(tmp_mat$UTC_S))))
                           for (yea in 1:length(numYea)) {
                             if (is.null(rawProduction[[numYea[yea]]])) {
                               rawProduction[[numYea[yea]]] <<- tmp_mat[years(tmp_mat$UTC_S) == numYea[yea], ]
                             } else {
                               rawProduction[[numYea[yea]]] <<- rbind(rawProduction[[numYea[yea]]], tmp_mat[years(tmp_mat$UTC_S) == numYea[yea], ])
                             }
                           }
                         }
                         cat("Done!", sep = "")
                       },
                       setFishPoinPara = function(speed_range, depth_range) {
                         fishPoinPara <<- data.frame(
                           min_spe = speed_range[1],
                           max_spe = speed_range[2],
                           min_dep = depth_range[1],
                           max_dep = depth_range[2]
                         )
                       },
                       setWeekMonthNum = function() {
                         cat("\nAdding week and month number to year ", sep = "")
                         for (j in names(rawEffort)) {
                           cat(j, "... ", sep = "")
                           tmp_dat <- rawEffort[[j]][, c("DATE")]
                           tmp_date <- as.Date(chron(tmp_dat))
                           rawEffort[[j]]$WeekNum <<- as.numeric(format(tmp_date, "%V"))
                           rawEffort[[j]]$MonthNum <<- as.numeric(format(tmp_date, "%m"))
                         }
                         cat("Done!", sep = "")
                       },
                       setFishPoin = function() {
                         cat("\nComputing fishing points year ", sep = "")
                         for (j in names(rawEffort)) {
                           cat(j, "... ", sep = "")
                           tmp_dat <- rawEffort[[j]][, c("SPE", "DEPTH")]
                           tmp_dat$FishSpeed <- tmp_dat$SPE >= as.numeric(fishPoinPara[1]) & tmp_dat$SPE <= as.numeric(fishPoinPara[2])
                           tmp_dat$FishDepth <- tmp_dat$DEPTH <= as.numeric(fishPoinPara[3]) & tmp_dat$DEPTH >= as.numeric(fishPoinPara[4])
                           rawEffort[[j]]$FishPoint <<- tmp_dat$FishSpeed & tmp_dat$FishDepth
                         }
                         cat("Done!", sep = "")
                       },
                       plotFishPoinStat = function() {
                         tmp_stat <- data.frame()
                         for (j in names(rawEffort)) {
                           tmp_sum <- sum(rawEffort[[j]]$FishPoint)
                           tmp_stat <- rbind(tmp_stat, cbind(j, c("Fishing", "Not Fishing"), rbind(tmp_sum, length(rawEffort[[j]]$FishPoint) - tmp_sum)))
                         }
                         colnames(tmp_stat) <- c("Year", "Status", "Value")
                         rownames(tmp_stat) <- NULL
                         tmp_stat$Value <- as.numeric(as.character(tmp_stat$Value))
                         print(tmp_stat)
                         fishStatPlot <- ggplot(data = tmp_stat, aes(x = Year, y = Value, fill = Status)) +
                           geom_bar(stat = "identity", position = position_dodge(), colour = "black") +
                           ggtitle("Number of Fishing Points each Year") +
                           scale_fill_manual(values = c("gainsboro", "grey50")) +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "right",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10),
                             legend.text = element_text(size = 8),
                             legend.title = element_text(size = 10)
                           )
                         
                         print(fishStatPlot)
                       },
                       plotSpeedDepth = function(which_year, speed_range, depth_range) {
                         tmp_dat <- rawEffort[[which_year]][, c("SPE", "DEPTH")]
                         op <- par(no.readonly = TRUE)
                         par(mfrow = c(2, 1), mar = c(3, 0, 2, 0))
                         speed_hist <- hist(tmp_dat$SPE[which(tmp_dat$SPE <= quantile(tmp_dat$SPE, 0.99) & tmp_dat$SPE > 0)], 100, plot = FALSE)
                         plot(speed_hist, xlab = "Speed", main = "")
                         abline(v = speed_range[1], col = "red", lty = 2)
                         abline(v = speed_range[2], col = "red", lty = 2)
                         text(x = speed_range[1] + ((speed_range[2] - speed_range[1]) / 2), y = max(speed_hist$counts) / 2, labels = "FISHING", col = 2)
                         title(main = paste("Speed/Depth profile of ", which_year, sep = ""))
                         depth_hist <- hist(tmp_dat$DEPTH[which(tmp_dat$DEPTH >= quantile(tmp_dat$DEPTH, 0.01) & tmp_dat$DEPTH <= 0)], 100, plot = FALSE)
                         plot(depth_hist, xlab = "Depth", main = "")
                         abline(v = depth_range[1], col = "red", lty = 2)
                         abline(v = depth_range[2], col = "red", lty = 2)
                         text(x = depth_range[1] + ((depth_range[2] - depth_range[1]) / 2), y = max(depth_hist$counts) / 2, labels = "FISHING", col = 2)
                         par(op)
                       },
                       setEffortIds = function() {
                         cat("\nSetting Effort IDs x year\n", sep = "")
                         effortIds <<- list()
                         for (i in names(rawEffort)) {
                           cat(i, "... ", sep = "")
                           tmp_ids <- unique(rawEffort[[i]][, 1])
                           tmp_key <- i
                           effortIds[[tmp_key]] <<- tmp_ids
                         }
                         effortIds[["All"]] <<- unique(unlist(effortIds))
                         cat("Done!\n", sep = "")
                       },
                       setProdSpec = function() {
                         prodSpec <<- list()
                         cat("\nSetting species list x year\n", sep = "")
                         for (i in names(effoProdMont)) {
                           cat(i, "... ", sep = "")
                           prodSpec[[i]] <<- colnames(effoProdMont[[i]])[ncol(dayEffoMatr[[i]]):ncol(effoProdMont[[i]])]
                           # if(i == names(prodSpec)[1]){
                           #   prodSpec[["Cross"]] <<- prodSpec[[i]]
                           # }else{
                           #   prodSpec[["Cross"]] <<- intersect(prodSpec[["Cross"]], prodSpec[[i]])
                           # }
                         }
                         prodSpec[["Cross"]] <<- sort(unique(unlist(prodSpec)))
                         setSpecSett()
                         setNNLS()
                         cat("Done!\n", sep = "")
                       },
                       setSpecSett = function() {
                         specSett <<- vector(mode = "list", length = length(prodSpec[["Cross"]]))
                         names(specSett) <<- sort(prodSpec[["Cross"]])
                       },
                       setNNLS = function() {
                         resNNLS <<- vector(mode = "list", length = length(prodSpec[["Cross"]]))
                         names(resNNLS) <<- sort(prodSpec[["Cross"]])
                       },
                       setBetaMeltYear = function(species) {
                         tmp_df <- data.frame(
                           Year = names(effoProd)[resNNLS[[species]]$SceMat$YEAR],
                           resNNLS[[species]]$bmat
                         )
                         tmp_df_agg <- aggregate(. ~ Year, tmp_df, sum)
                         betaMeltYear[[species]] <<- melt(
                           data = tmp_df_agg, id.vars = "Year",
                           measure.vars = c(2:ncol(tmp_df_agg)),
                           variable.name = "FishGround", value.name = "Productivity"
                         )
                       },
                       setProdMeltYear = function(species) {
                         tmp_df <- data.frame(
                           Year = as.character(effoAllLoa[, 1]),
                           predProd[[species]]
                         )
                         tmp_df_agg <- aggregate(. ~ Year, tmp_df, sum)
                         prodMeltYear[[species]] <<- melt(
                           data = tmp_df_agg, id.vars = "Year",
                           measure.vars = c(2:ncol(tmp_df_agg)),
                           variable.name = "FishGround", value.name = "Production"
                         )
                       },
                       plotTotProd = function(species) {
                         year_Prod <- aggregate(. ~ Year, prodMeltYear[[species]][, c(1, 3)], sum)
                         year_Prod[, 1] <- as.numeric(as.character(year_Prod[, 1]))
                         tot_prod <- ggplot_TotalProduction(year_Prod)
                         fg_prod <- ggplot_FGProduction(prodMeltYear[[species]])
                         grid.arrange(tot_prod,
                                      fg_prod,
                                      layout_matrix = rbind(c(1, 1, 2), c(1, 1, 2)),
                                      top = "Production"
                         )
                       },
                       plotNNLS = function(species, thresR2) {
                         tmp_df <- data.frame(
                           R2 = "R2",
                           Values = as.numeric(resNNLS[[species]][["nnls_r2"]])
                         )
                         bp <- suppressMessages(
                           ggplot(tmp_df, aes(x = R2, y = Values)) +
                             geom_violin(fill = "grey30", colour = "grey90", alpha = 0.05) +
                             geom_boxplot(fill = "grey90", width = 0.5) +
                             stat_boxplot(geom = "errorbar", width = 0.25) +
                             theme(axis.text.x = element_blank()) +
                             ylim(0, 1) +
                             labs(title = "R^2 values") +
                             geom_hline(aes(yintercept = thresR2),
                                        linetype = "dashed", size = 0.5, colour = "red"
                             ) +
                             theme_tufte(base_size = 14, ticks = F) +
                             theme(
                               legend.position = "none",
                               plot.title = element_text(size = 14),
                               axis.text.x = element_text(size = 8),
                               axis.title = element_blank(),
                               panel.grid = element_line(size = 0.1, linetype = 2, colour = "grey20"),
                               axis.text.y = element_text(size = 10),
                               axis.ticks.y = element_blank()
                             )
                         )
                         tmp_reg <- data.frame(
                           Observed = resNNLS[[species]]$obsY,
                           Fitted = resNNLS[[species]]$fittedY
                         )
                         
                         reg_p <- suppressMessages(
                           ggplot(tmp_reg, aes(y = Fitted, x = Observed)) +
                             geom_point(alpha = 0.25, size = 0.3) + stat_smooth(method = "lm") +
                             labs(title = "Observed VS Fitted") +
                             scale_x_log10() +
                             scale_y_log10() +
                             annotation_logticks() +
                             theme_tufte(base_size = 14) +
                             theme(
                               legend.position = "none",
                               plot.title = element_text(size = 14),
                               axis.text.x = element_text(size = 8),
                               axis.title = element_blank(),
                               panel.grid = element_line(size = 0.1, linetype = 2, colour = "grey20"),
                               axis.text.y = element_text(size = 10),
                               axis.ticks.y = element_blank()
                             )
                         )
                         suppressWarnings(
                           grid.arrange(reg_p, bp, layout_matrix = rbind(c(1, 1, 2), c(1, 1, 2)))
                         )
                       },
                       setSpecSettItm = function(species, thresh, brea, max_xlim) {
                         specSett[[species]] <<- data.frame(
                           threshold = thresh,
                           breaks = brea,
                           max_x = max_xlim
                         )
                       },
                       plotLogitROC = function(selSpecie) {
                         plot(specLogit[[selSpecie]]$logit$Roc,
                              print.cutoffs.at = seq(0, 1, 0.1),
                              text.adj = c(-0.2, 1.7), main = "Logit ROC results"
                         )
                       },
                       setSpecLogitConf = function(selSpecie, cutoff = specLogit[[selSpecie]]$logit$Cut) {
                         predict <- factor(as.numeric(specLogit[[selSpecie]]$logit$Predict > cutoff))
                         truth <- factor(1 * (specLogit[[selSpecie]]$Landings[-specLogit[[selSpecie]]$logit$Split] > specSett[[selSpecie]]$threshold))
                         tmp_Tbl <- table(predict, truth)
                         specLogit[[selSpecie]]$logit$Confusion <<- caret::confusionMatrix(tmp_Tbl)
                       },
                       setLogitTrain = function(selSpecie, train, cp_val = 0.01, cv_val = 2) {
                         specLogit[[selSpecie]]$logit$Model <<- switch(specLogit[[selSpecie]]$logit$Name,
                                                                       GLM = {
                                                                         glm(Target ~ ., family = binomial(logit), data = train)
                                                                       },
                                                                       CART = {
                                                                         rpart::rpart(Target ~ .,
                                                                                      data = train, method = "class",
                                                                                      control = rpart.control(cp = cp_val)
                                                                         )
                                                                       },
                                                                       RF = {
                                                                         caret::train(Target ~ .,
                                                                                      data = train, method = "rf",
                                                                                      trControl = trainControl(method = "cv", number = cv_val),
                                                                                      prox = TRUE, allowParallel = TRUE, metric = "Kappa",
                                                                                      maximize = TRUE
                                                                         )
                                                                       },
                                                                       NN = {   }
                         )
                       },
                       setLogitTest = function(selSpecie, test) {
                         specLogit[[selSpecie]]$logit$Predict <<- switch(specLogit[[selSpecie]]$logit$Name,
                                                                         GLM = {
                                                                           predict(specLogit[[selSpecie]]$logit$Model,
                                                                                   newdata = test, type = "response"
                                                                           )
                                                                         },
                                                                         CART = {
                                                                           predict(specLogit[[selSpecie]]$logit$Model,
                                                                                   newdata = test, type = "prob"
                                                                           )[, 2]
                                                                         },
                                                                         RF = {
                                                                           predict(specLogit[[selSpecie]]$logit$Model,
                                                                                   newdata = test, type = "prob"
                                                                           )[, 2]
                                                                         },
                                                                         NN = {   }
                         )
                       },
                       setLogitPred = function(selSpecie, test) {
                         specLogit[[selSpecie]]$logit$Prediction <<- ROCR::prediction(specLogit[[selSpecie]]$logit$Predict, test$Target)
                       },
                       setLogitCut = function(selSpecie) {
                         perf <- ROCR::performance(specLogit[[selSpecie]]$logit$Prediction, "acc")
                         specLogit[[selSpecie]]$logit$Cut <<- perf@x.values[[1]][which.max(perf@y.values[[1]])]
                       },
                       setLogitRoc = function(selSpecie) {
                         specLogit[[selSpecie]]$logit$Roc <<- ROCR::performance(specLogit[[selSpecie]]$logit$Prediction, "tpr", "fpr")
                       },
                       setLogitConf = function(selSpecie, test) {
                         tryCatch(
                           expr = {
                             specLogit[[selSpecie]]$logit$Confusion <<- caret::confusionMatrix(
                               factor(specLogit[[selSpecie]]$logit$Predict > specLogit[[selSpecie]]$logit$Cut, levels = c(FALSE, TRUE)),
                               test$Target
                             )
                           },
                           error = function(error_message) {
                             message("An error has occurred, different categories to compute the confusion matrix")
                             message(error_message)
                           }
                         )
                       },
                       setSpecLogit = function(selSpecie, selModel = c("GLM", "CART", "RF")[1],
                                               cp = 0.01, cv = 2) {
                         if (is.null(specLogit)) specLogit <<- list()
                         if (is.null(specLogit[[selSpecie]])) specLogit[[selSpecie]] <<- list()
                         
                         tmp_mat <- getMatSpeLand(selSpecie)
                         
                         colnames(tmp_mat)[ncol(tmp_mat)] <- "Target"
                         specLogit[[selSpecie]]$Landings <<- tmp_mat$Target
                         tmp_mat$Target <- factor(tmp_mat$Target > specSett[[selSpecie]]$threshold, levels = c(FALSE, TRUE))
                         
                         split <- caret::createDataPartition(y = tmp_mat$Target, p = 0.7, list = FALSE)[, 1]
                         train <- tmp_mat[split, ]
                         test <- tmp_mat[-split, ]
                         specLogit[[selSpecie]]$logit$Split <<- split
                         specLogit[[selSpecie]]$logit$Name <<- selModel
                         # Train
                         setLogitTrain(selSpecie, train, cp_val = cp, cv_val = cv)
                         # Test
                         setLogitTest(selSpecie, test)
                         # Prediction
                         setLogitPred(selSpecie, test)
                         # Cutoff
                         setLogitCut(selSpecie)
                         # ROC
                         setLogitRoc(selSpecie)
                         # Confusion
                         setLogitConf(selSpecie, test)
                       },
                       getMatSpeLand = function(species) {
                         tmp_mat <- effoProdAllLoa[, c(
                           3:(ncol(dayEffoMatr[[1]])),
                           which(colnames(effoProdAllLoa) == "Loa"),
                           which(colnames(effoProdAllLoa) == species)
                         )]
                         tmp_mat$MonthNum <- as.factor(tmp_mat$MonthNum)
                         return(tmp_mat)
                       },
                       setEffoProdAll = function() {
                         cat("\nSetting Effort x Production year\n", sep = "")
                         tmp_spe <- sort(prodSpec[["Cross"]])
                         for (i in names(effoProdMont)) {
                           cat(i, "... ", sep = "")
                           tmp_nam <- colnames(effoProdMont[[i]])
                           tmp_cols <- which(tmp_nam %in% tmp_spe)
                           if (i == names(effoProdMont)[1]) {
                             effoProdAll <<- cbind(Year = i, effoProdMont[[i]][, c(1:(ncol(dayEffoMatr[[i]]) - 1), tmp_cols[order(tmp_nam[tmp_cols])])])
                           } else {
                             effoProdAll <<- rbind(effoProdAll, cbind(Year = i, effoProdMont[[i]][, c(1:(ncol(dayEffoMatr[[i]]) - 1), tmp_cols[order(tmp_nam[tmp_cols])])]))
                           }
                         }
                         cat("Done!\n", sep = "")
                       },
                       setEffoAll = function() {
                         cat("\nSetting Effort x Year\n", sep = "")
                         for (i in names(effoMont)) {
                           cat(i, "... ", sep = "")
                           if (i == names(effoMont)[1]) {
                             effoAll <<- cbind(Year = i, effoMont[[i]])
                           } else {
                             effoAll <<- rbind(effoAll, cbind(Year = i, effoMont[[i]]))
                           }
                         }
                         cat("Done!", sep = "")
                       },
                       setEffoProdAllLoa = function() {
                         tmp_effoProd <- effoProdAll
                         tmp_loa <- rawRegister[, c("CFR", "Loa")]
                         tmp_loa$CFR <- substr(tmp_loa$CFR, 4, nchar(tmp_loa$CFR[1]))
                         names(tmp_loa) <- c("I_NCEE", "Loa")
                         tmp_allLoa <- sqldf("select * from tmp_effoProd left join (select * from tmp_loa) using (I_NCEE)")
                         naFind <- which(is.na(tmp_allLoa), arr.ind = TRUE)
                         if (length(naFind) > 0) {
                           effoProdAllLoa <<- tmp_allLoa[-naFind[, 1], ]
                         } else {
                           effoProdAllLoa <<- tmp_allLoa
                         }
                       },
                       setEffoAllLoa = function() {
                         cat("\nSetting Effort LOA... ", sep = "")
                         tmp_effo <- effoAll
                         tmp_loa <- rawRegister[, c("CFR", "Loa")]
                         tmp_loa$CFR <- substr(tmp_loa$CFR, 4, nchar(tmp_loa$CFR[1]))
                         names(tmp_loa) <- c("I_NCEE", "Loa")
                         tmp_Loa <- sqldf("select * from tmp_effo left join (select * from tmp_loa) using (I_NCEE)")
                         naFind <- which(is.na(tmp_Loa), arr.ind = TRUE)
                         if (length(naFind) > 0) {
                           effoAllLoa <<- tmp_Loa[-naFind[, 1], ]
                         } else {
                           effoAllLoa <<- tmp_Loa
                         }
                         cat("Done!\n", sep = "")
                       },
                       setProdIds = function() {
                         cat("\nSetting Production IDs x year\n", sep = "")
                         productionIds <<- list()
                         for (i in names(rawProduction)) {
                           cat(i, "... ", sep = "")
                           tmp_ids <- unique(rawProduction[[i]][, 1])
                           tmp_key <- i
                           productionIds[[tmp_key]] <<- tmp_ids
                         }
                         productionIds[["All"]] <<- unique(unlist(productionIds))
                         cat("Done!\n", sep = "")
                       },
                       setIdsEffoProd = function() {
                         ###   Set IDs cross match effort/production
                         to_match <- names(effortIds)[names(effortIds) %in% names(productionIds)]
                         idsEffoProd <<- list()
                         for (i in to_match) {
                           idsEffoProd[[i]] <<- effortIds[[i]][effortIds[[i]] %in% productionIds[[i]]]
                         }
                       },
                       plotCountIDsEffoProd = function() {
                         tmp_effo <- data.frame(
                           "Year" = names(effortIds),
                           "Ids" = unlist(lapply(sapply(effortIds, unique, simplify = FALSE), length)),
                           "Dataset" = "Effort"
                         )
                         tmp_prod <- data.frame(
                           "Year" = names(productionIds),
                           "Ids" = unlist(lapply(sapply(productionIds, unique, simplify = FALSE), length)),
                           "Dataset" = "Production"
                         )
                         tmp_comb <- data.frame(
                           "Year" = names(idsEffoProd),
                           "Ids" = unlist(lapply(sapply(idsEffoProd, unique, simplify = FALSE), length)),
                           "Dataset" = "Overlap"
                         )
                         tmp_df <- rbind(tmp_effo, tmp_prod, tmp_comb)
                         rownames(tmp_df) <- NULL
                         tmp_plot <- ggplot(tmp_df, aes(x = Year, y = Ids, fill = Dataset)) +
                           geom_bar(position = position_dodge(), stat = "identity", alpha = 0.75) +
                           geom_text(aes(y = Ids, label = Ids),
                                     position = position_dodge(width = 1),
                                     vjust = 2.5, color = "grey20"
                           ) +
                           scale_fill_brewer(palette = "Set2") +
                           ggtitle("Count of Distinct Vessels") +
                           ylab("N. of IDs") +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "right",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.1, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10),
                             legend.text = element_text(size = 8),
                             legend.title = element_text(size = 10)
                           )
                         print(tmp_plot)
                       },
                       plotCountIDsEffo = function() {
                         tmp_df <- data.frame(
                           "Year" = names(effortIds),
                           "Ids" = unlist(lapply(sapply(effortIds, unique), length))
                         )
                         tmp_plot <- ggplot(tmp_df, aes(x = Year, y = Ids)) + geom_bar(stat = "identity") +
                           geom_text(aes(y = Ids, label = Ids),
                                     position = position_dodge(width = 1),
                                     vjust = 2.5, color = "white"
                           ) +
                           ggtitle("Count of Distinct Vessels - Effort Dataset") +
                           ylab("N. of IDs") +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "none",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10),
                             legend.text = element_text(size = 8),
                             legend.title = element_text(size = 10)
                           )
                         print(tmp_plot)
                       },
                       plotCountIDsProd = function() {
                         tmp_df <- data.frame(
                           "Year" = names(productionIds),
                           "Ids" = unlist(lapply(unique(productionIds), length))
                         )
                         tmp_plot <- ggplot(tmp_df, aes(x = Year, y = Ids)) + geom_bar(stat = "identity") +
                           geom_text(aes(y = Ids, label = Ids),
                                     position = position_dodge(width = 1),
                                     vjust = 2.5, color = "white"
                           ) +
                           ggtitle("Count of Distinct Vessels - Production Dataset") +
                           ylab("N. of IDs") +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "none",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10),
                             legend.text = element_text(size = 8),
                             legend.title = element_text(size = 10)
                           )
                         print(tmp_plot)
                       },
                       setEffoProdMatr = function() {
                         effoProd <<- list()
                         cat("\nCreating Effort x Production matrix\n", sep = "")
                         for (i in intersect(names(dayEffoMatr), names(prodMatr))) {
                           cat(i, "... ", sep = "")
                           tmp_effo <- dayEffoMatr[[i]]
                           tmp_prod <- prodMatr[[i]]
                           tmp_effoProd <- merge(x = cbind(tmp_effo, NUMUE = tmp_effo$I_NCEE), y = tmp_prod, by.x = "I_NCEE", by.y = "NUMUE")
                           tmp_ids <- which(tmp_effoProd$DATE >= tmp_effoProd$UTC_S & tmp_effoProd$DATE <= tmp_effoProd$UTC_E)
                           effoProd[[i]] <<- tmp_effoProd[tmp_ids,]
                           rm(tmp_effo, tmp_prod, tmp_effoProd)
                           gc()
                         }
                         cat("Done!\n")
                       },
                       setEffoProdMont = function() {
                         effoProdMont <<- list()
                         cat("\nMatching Effort x FG with Production\n", sep = "")
                         for (i in names(effoProd)) {
                           cat(i, "... ", sep = "")
                           dis_vesmon <- unique(effoProd[[i]][, c("I_NCEE", "MonthNum")])
                           effoProdMont[[i]] <<- data.frame(matrix(data = 0, nrow = nrow(dis_vesmon), ncol = ncol(effoProd[[i]]) - 5))
                           colnames(effoProdMont[[i]]) <<- c(colnames(dayEffoMatr[[i]])[-2], colnames(prodMatr[[i]])[-c(1:4)])
                           effoProdMont[[i]][, 1:2] <<- dis_vesmon
                           for (j in 1:nrow(dis_vesmon)) {
                             tmp_itm <- effoProd[[i]][which(effoProd[[i]]$I_NCEE == dis_vesmon[j, 1] & effoProd[[i]]$MonthNum == dis_vesmon[j, 2]), ]
                             effoProdMont[[i]][j, 3:(ncol(dayEffoMatr[[i]]) - 1)] <<- apply(unique(tmp_itm[, 4:ncol(dayEffoMatr[[i]])]), 2, sum)
                             tmp_prod_itm <- unique(tmp_itm[, c(ncol(dayEffoMatr[[i]]) + 1, (ncol(dayEffoMatr[[i]]) + 5):ncol(tmp_itm))])
                             effoProdMont[[i]][j, (ncol(dayEffoMatr[[i]])):ncol(effoProdMont[[i]])] <<- apply(tmp_prod_itm[, 2:ncol(tmp_prod_itm)], 2, sum)
                           }
                         }
                         cat("Done!\n")
                       },
                       setEffoMont = function() {
                         effoMont <<- list()
                         cat("\nAggregating year by month\n", sep = "")
                         for (i in names(dayEffoMatr)) {
                           cat(i, "... ", sep = "")
                           dis_vesmon <- unique(dayEffoMatr[[i]][, c("I_NCEE", "MonthNum")])
                           effoMont[[i]] <<- data.frame(matrix(data = 0, nrow = nrow(dis_vesmon), ncol = ncol(dayEffoMatr[[i]]) - 1))
                           colnames(effoMont[[i]]) <<- c(colnames(dayEffoMatr[[i]])[-2])
                           effoMont[[i]][, 1:2] <<- dis_vesmon
                           for (j in 1:nrow(dis_vesmon)) {
                             tmp_itm <- dayEffoMatr[[i]][which(dayEffoMatr[[i]]$I_NCEE == dis_vesmon[j, 1] & dayEffoMatr[[i]]$MonthNum == dis_vesmon[j, 2]), ]
                             effoMont[[i]][j, 3:(ncol(dayEffoMatr[[i]]) - 1)] <<- apply(unique(tmp_itm[, 4:ncol(dayEffoMatr[[i]])]), 2, sum)
                           }
                         }
                         cat("Done!")
                       },
                       setProdMatr = function() {
                         prodMatr <<- list()
                         for (i in names(rawProduction)) {
                           tmp_prod <- rawProduction[[i]]
                           tmp_prod <- tmp_prod[tmp_prod$NUMUE %in% idsEffoProd[[i]], ]
                           tmp_matrix <- dcast(tmp_prod,
                                               NUMUE + UTC_S + UTC_E ~ SPECIES,
                                               fun.aggregate = sum,
                                               na.rm = TRUE, value.var = "KGS"
                           )
                           tmp_matrix <- cbind("prodID" = 1:nrow(tmp_matrix), tmp_matrix)
                           prodMatr[[i]] <<- tmp_matrix
                         }
                       },
                       setDayEffoMatrGround = function(maxFG = max(rawEffort[[1]]$FishGround[!is.na(rawEffort[[1]]$FishGround)])) {
                         dayEffoMatr <<- list()
                         cat("\nCreating daily effort x Fishing Ground matrix\n", sep = "")
                         for (j in names(rawEffort)) {
                           cat(j, "... ", sep = "")
                           tmp_dat <- rawEffort[[j]][rawEffort[[j]]$FishPoint == TRUE & rawEffort[[j]]$P_INT == 1 & !is.na(rawEffort[[j]]$Cell), c("I_NCEE", "DATE", "MonthNum", "FishGround", "FishPoint")]
                           tmp_dat$DATE <- ceiling(tmp_dat$DATE)
                           tmp_matrix <- dcast(tmp_dat,
                                               formula = I_NCEE + DATE + MonthNum ~ FishGround,
                                               fun.aggregate = sum, na.rm = TRUE, value.var = "FishPoint"
                           )
                           ## points to hours: interpolation interval 10 min
                           tmp_matrix[, 4:ncol(tmp_matrix)] <- tmp_matrix[, 4:ncol(tmp_matrix)] / 6
                           # miss_cols <- setdiff(as.character(unique(rawEffort[[j]]$FishGround[!is.na(rawEffort[[j]]$FishGround)])),
                           #                      names(tmp_matrix)[4:ncol(tmp_matrix)])
                           miss_cols <- setdiff(
                             1:maxFG,
                             names(tmp_matrix)[4:ncol(tmp_matrix)]
                           )
                           if (length(miss_cols) > 0) {
                             # tmp_matrix[,miss_cols] <- 0
                             tmp_matrix[, paste(miss_cols)] <- 0
                             # tmp_matrix <- tmp_matrix[,c(1:4, 4+order(as.numeric(names(tmp_matrix)[5:ncol(tmp_matrix)])))]
                             tmp_matrix <- tmp_matrix[, c(1:3, 3 + order(as.numeric(names(tmp_matrix)[4:ncol(tmp_matrix)])))]
                           }
                           dayEffoMatr[[j]] <<- tmp_matrix
                         }
                         cat("Done!\n", sep = "")
                       },
                       getLoa4Prod = function() {
                         if (!is.null(productionIds) & !is.null(registerIds)) {
                           tmp_reg <- rawRegister[, c("CFR", "Loa")]
                           tmp_reg[, 1] <- substr(tmp_reg[, 1], 4, nchar(tmp_reg[1, 1]))
                           
                           tmp_pro <- data.frame("CFR" = productionIds[["All"]])
                           tmp_pro[, 1] <- paste(unlist(lapply(mapply(rep, times = 9 - nchar(tmp_pro[, 1]), x = 0), paste, collapse = "")),
                                                 tmp_pro[, 1],
                                                 sep = ""
                           )
                           
                           prodIdsLoa <<- merge(x = tmp_pro, y = tmp_reg, by = "CFR")
                         }
                       },
                       plotLoaProd = function() {
                         tmp_tab <- table(round(prodIdsLoa[, 2]))
                         tmp_df <- data.frame(
                           "Length" = names(tmp_tab),
                           "Count" = as.numeric(tmp_tab)
                         )
                         ggplot(tmp_df, aes(x = Length, y = Count)) + geom_bar(stat = "identity") +
                           geom_text(aes(y = Count, label = Count),
                                     position = position_dodge(width = 1),
                                     vjust = -.5, color = "black"
                           ) +
                           ggtitle("Count of Distinct Vessels") +
                           ylab("N. of IDs")
                       },
                       readRegisterEU = function(reg_path) {
                         cat("\n\tChecking EU Fleet Register format...", sep = "")
                         two_rows <- readLines(con = reg_path, n = 2)
                         last_char <- substr(two_rows[2], nchar(two_rows[2]), nchar(two_rows[2]))
                         raw_fleet <- readLines(con = reg_path, n = -1)
                         if (last_char == ";") {
                           cat("\n\tTrailing character found! Cleaning...", sep = "")
                           tmp_flee <- paste(unlist(lapply(strsplit(raw_fleet, split = ";"), paste, collapse = ";")), collapse = "\n")
                         } else {
                           tmp_flee <- paste(raw_fleet, collapse = "\n")
                         }
                         if (substr(reg_path, nchar(reg_path) - 12, nchar(reg_path)) != "_smart-ed.csv") {
                           tmp_flee <- gsub("\\;", ",", tmp_flee)
                           new_path <- paste(substr(reg_path, 1, nchar(reg_path) - 4), "_smart-ed",
                                             substr(reg_path, nchar(reg_path) - 3, nchar(reg_path)),
                                             sep = ""
                           )
                           cat("\nWriting edited Fleet register in:\n", new_path, "\n", sep = "")
                           write(tmp_flee, file = new_path)
                           reg_path <- new_path
                         }
                         cat("\nLoading file... ", sep = "")
                         re_fleet <- read.csv(reg_path, stringsAsFactors = FALSE)
                         cat("OK", sep = "")
                         return(re_fleet)
                       },
                       cleanRegister = function() {
                         cat("\n\tOrdering Fleet Register by CFR... ", sep = "")
                         rawRegister$CFR <<- as.character(rawRegister$CFR)
                         rawRegister <<- rawRegister[order(rawRegister$CFR), ]
                         rawRegister$Country.Code <<- as.character(rawRegister$Country.Code)
                         setRegIds()
                         rawRegister$Loa <<- as.numeric(as.character(rawRegister$Loa))
                         rawRegister$Power.Main <<- as.numeric(as.character(rawRegister$Power.Main))
                         cat("OK\n", sep = "")
                       },
                       plotRegSum = function() {
                         cat("Plotting Fleet register summary statistics... ", sep = "")
                         def.par <- par(no.readonly = TRUE)
                         layout(matrix(c(1, 2, 3, 4, 5, 6), 2, 3, byrow = TRUE))
                         plotBarReg(regVar = "Gear.Main.Code", title = "Main Gear")
                         plotBarReg(regVar = "Gear.Sec.Code", title = "Secondary Gear")
                         plotBarReg(regVar = "Hull.Material", title = "Hull Material")
                         plotBoxReg(regVar = "Construction.Year", title = "Year of Construction")
                         plotBoxReg(regVar = "Loa", title = "Length Overall")
                         plotBoxReg(regVar = "Power.Main", title = "Power")
                         par(def.par)
                         cat("Completed\n", sep = "")
                       },
                       plotBarReg = function(regVar, p_las = 2, title = regVar) {
                         barplot(table(rawRegister[, regVar]), las = p_las, main = title)
                       },
                       plotBoxReg = function(regVar, title = regVar) {
                         boxplot(rawRegister[, regVar], main = title)
                       },
                       setRegIds = function() {
                         registerIds <<- rawRegister$CFR
                       }
                     )
)



#### SampleMap################################################
#' SampleMap
#'
#' The \code{SampleMap} class implements the class of SMART
#' to control geographical data.
#'
#' @docType class
#' @usage NULL
#' @keywords data
#' @return Object of \code{\link{R6Class}} with attributes and methods for the Environmental data.
#'
#' @format \code{\link{R6Class}} object.
#'
#' @field gridPath Stores the file path to the selected Environment grid shapefile.
#' @field gridName Stores the file name of the Environment grid shapefile.
#' @field gridShp Stores the SpatialPoligon object of the Environment grid shapefile.
#' @field gridBbox Stores the bounding box coordinates of the Environment grid shapefile.
#' @field gridBboxExt Stores the extended bounding box coordinates of the Environment grid shapefile.
#' @field gridBboxSP Stores the bounding box of the Environment grid as a SpatialPoligon.
#' @field areaGrid Stores the total area covered by the Environment grid.
#' @field areaStrata Stores the area covered by depth strata.
#' @field weightStrata Stores the area covered by depth strata relative to the total area.
#' @field harbDbf Stores coordinates and names of the harbours.
#' @field bioPath Stores the file path of the Substrate map.
#' @field bioName Stores the file name of the Substrate map.
#' @field bioShp Stores the SpatialPoligon object of the Substrate map.
#' @field bioDF Stores the data.frame representation of the Substrate map.
#' @field gridPolySet Stores the PolySet object of the Environment grid.
#' @field gridFortify Stores the fortified SpatialPoligon object of the Environment grid.
#' @field nCells Stores the number of cells in the Environment grid.
#' @field griCent Stores the coordinates of the cells' centroids.
#' @field gridBathy Stores the bathymetric matrix.
#' @field centDept Stores the depth of the cells' centroids.
#' @field clusInpu Stores the input data for the spatial clustering.
#' @field clusMat Stores the results of the spatial clustering.
#' @field indSil Stores the silhouette index of the spatial clustering result.
#' @field cutFG Stores the number of cuts for the spatial clustering.
#' @field availData Stores the names of the variables for the spatial clustering.
#' @field rawInpu Stores the raw input for the spatial clustering.
#' @field cutResult Stores the summary data of the spatial clustering result.
#' @field cutResEffo Stores the average effort data from the spatial clustering result.
#' @field cutResShp Stores the SpatialPoligon object of the spatial clustering result.
#' @field cutResShpCent Stores the coordinates of the clusters' centroids.
#' @field cutResShpFort Stores the fortified SpatialPoligon object of the spatial clustering result.
#' @field fgWeigDist Stores the weighted distance between harbours and fishing grounds.
#' @field ggBioDF Stores the plot of the Substrate map.
#' @field ggDepth Stores the plot of the Bathymetric map.
#' @field ggDepthFGbox Stores the boxplot of the depth values of each fishing ground.
#' @field ggEffoFGbox Stores the boxplot of the effort values of each fishing ground.
#' @field ggEffoFGmap Stores the plot of the Effort map.
#' @field ggBioFGmat Stores the tilemap of the substrate values of each fishing ground.
#' @field ggCutFGmap Stores the plot of the Fishing ground configuration.
#' @field ggSilFGlin Stores the plot of the silhouette index.
#' @field ggBetaFGmap Stores the plot of the Productivity map.
#' @field ggBetaFGbox Stores the boxplot of the Productivity values of each fishing ground.
#' @field ggProdFGmap Stores the plot of the Production map.
#' @field ggProdFGbox Stores the boxplot of the Production values of each fishing ground.
#' @field ggMapFgFishery Stores the plot of the Fishery data coordinates.
#' @field ggMapFgSurvey Stores the plot of the Survey data coordinates.
#' @field gooMap Stores the satellite view of the area of study.
#' @field gooMapPlot Stores the satellite plot of the area of study.
#' @field gooGrid Stores the plot of the Environment Grid.
#' @field gooBbox Stores the plot of the Bounding Box of the Environment Grid.
#' @field sampColScale Stores the color scale for the species plots.
#' @field plotRange Stores the plot ranges for the Environmental Grid.
#'
#' @section Methods:
#' \describe{
#'   \item{\code{setAreaGrid()}}{This method is used to compute the total area covered by the environmental grid.}
#'   \item{\code{setAreaStrata(vectorStrata)}}{This method is used to compute the area covered by each depth strata.}
#'   \item{\code{setWeightStrata()}}{This method is used to compute the area covered by each depth strata relative to the total area of the grid.}
#'   \item{\code{loadHarbDbf(dbf_path)}}{This method is used to load a dbf file of coordinates and harbours names.}
#'   \item{\code{set_ggMapFgSurvey(rawSampCoo)}}{This method is used to setup the plot of the spatial distribution of survey data.}
#'   \item{\code{set_ggMapFgFishery(rawSampCoo)}}{This method is used to setup the plot of the spatial distribution of fishery data.}
#'   \item{\code{createGridBbox()}}{This method is used to setup the bounding box of the environment grid.}
#'   \item{\code{getGooMap()}}{This method is used to retrieve the satellite view of the area of study.}
#'   \item{\code{setGooPlot()}}{This method is used to setup the base plot of the area of study.}
#'   \item{\code{setPlotRange()}}{This method is used to setup the ranges of the base plot.}
#'   \item{\code{setGooGrid()}}{This method is used to setup the plot of the environment grid.}
#'   \item{\code{plotGooGrid()}}{This method is used to plot the environment grid.}
#'   \item{\code{plotGooGridData(grid_data)}}{This method is used to plot the environment grid.}
#'   \item{\code{setSampColScale(fac_col)}}{This method is used to setup the color scale for the species' plots.}
#'   \item{\code{plotGooSpeSur(poi_data)}}{This method is used to plot the spatial distribution of the survey data}
#'   \item{\code{plotGooSpeFis(poi_data)}}{This method is used to plot the spatial distribution of the fishery data}
#'   \item{\code{setGooBbox()}}{This method is used to setup the bounding box of the environment grid.}
#'   \item{\code{plotGooBbox()}}{This method is used to plot the bounding box of the environment grid.}
#'   \item{\code{setGridPath(path2grid)}}{This method is used to store the path to the grid file.}
#'   \item{\code{setGridName()}}{This method is used to store the name of the grid file.}
#'   \item{\code{loadGridShp()}}{This method is used to load the grid file.}
#'   \item{\code{setBioPath()}}{This method is used to store the path of the seabed substrates file.}
#'   \item{\code{setBioName()}}{This method is used to store the name of the seabed substrates file.}
#'   \item{\code{loadBioShp()}}{This method is used to load the seabed substrates file.}
#'   \item{\code{addBioShp(bio_path)}}{This method is used store the path and name of the seabed file and then load the SpatialPoligon object.}
#'   \item{\code{loadBioDF(bio_path)}}{This method is used to load a Data.Frame of seabed substrates.}
#'   \item{\code{plotBioDF()}}{This method is used to plot the map of substrates.}
#'   \item{\code{setGgBioDF()}}{This method is used to setup the plot of substrates.}
#'   \item{\code{ggplotBioDF()}}{This method is used to plot the map of substrates.}
#'   \item{\code{createPolySet()}}{This method is used to store the PolySet object of the Environment grid.}
#'   \item{\code{fortifyGridShp()}}{This method is used to fortify the SpatialPolygon od the Environment grid.}
#'   \item{\code{setNumCell()}}{This method is used to setup the number of cells in the grid.}
#'   \item{\code{setGridCenter()}}{This method is used to store the coordinates of cells centroids.}
#'   \item{\code{getGridBath()}}{This method is used to retrieve the bathymetric matrix.}
#'   \item{\code{saveGridBath(bathy_path)}}{This method is used to save the bathymetric matrix to file.}
#'   \item{\code{loadGridBath(bathy_path)}}{This method is used to load the bathymetric matrix from file.}
#'   \item{\code{getCentDept()}}{This method is used to assign the depth to each cell centroids.}
#'   \item{\code{setGgDepth(isoLine)}}{This method is used to setup the bathymetric plot.}
#'   \item{\code{ggplotGridBathy()}}{This method is used to plot the bathymetric map.}
#'   \item{\code{plotSamMap(title, celCol)}}{This method is used to plot the Environment map.}
#'   \item{\code{plotCoho(abbs)}}{This method is used to plot the spatial distribution of the species.}
#'   \item{\code{setClusInpu(whiData, howData)}}{This method is used to setup the input for the spatial clustering}
#'   \item{\code{calcFishGrou(numCuts, minsize, maxsize, modeska, skater_method, nei_queen)}}{This method is used to run the spatial clustering routine}
#'   \item{\code{plotFishGrou(ind_clu)}}{This method is used to plot the fishing ground configuration}
#'   \item{\code{setCutResult(ind_clu)}}{This method is used to choose a fishing ground configuration}
#'   \item{\code{setDepthFGbox()}}{This method is used to setup the boxplot of depth by fishing ground}
#'   \item{\code{setEffoFGbox()}}{This method is used to setup the boxplot of effort by fishing ground}
#'   \item{\code{setEffoFGmap()}}{This method is used to setup the map of effort by fishing ground}
#'   \item{\code{setBioFGmat()}}{This method is used to setup the tileplot of substrate by fishing ground}
#'   \item{\code{setCutFGmap()}}{This method is used to plot the fishing ground map}
#'   \item{\code{setSilFGlin(numCut)}}{This method is used to setup the plot of the Silhouette index}
#'   }

SampleMap <- R6Class("sampleMap",
                     portable = FALSE,
                     class = TRUE,
                     public = list(
                       gridPath = NULL,
                       gridName = NULL,
                       gridShp = NULL,
                       gridBbox = NULL,
                       gridBboxExt = NULL,
                       gridBboxSP = NULL,
                       areaGrid = NULL,
                       areaStrata = NULL,
                       weightStrata = NULL,
                       harbDbf = NULL,
                       bioPath = NULL,
                       bioName = NULL,
                       bioShp = NULL,
                       bioDF = NULL,
                       gridPolySet = NULL,
                       gridFortify = NULL,
                       nCells = NULL,
                       griCent = NULL, # gCenter
                       gridBathy = NULL,
                       centDept = NULL,
                       clusInpu = NULL, # calcfish input
                       clusMat = NULL, # matrix output calcfish
                       indSil = NULL, # vect clusters silhouette output calcfish
                       # indCH = NULL, # vect index CH output calcfish
                       cutFG = NULL,
                       availData = NULL,
                       rawInpu = NULL,
                       cutResult = NULL,
                       cutResEffo = NULL,
                       cutResShp = NULL,
                       cutResShpCent = NULL,
                       cutResShpFort = NULL,
                       fgWeigDist = NULL,
                       gooEnv = NULL,
                       ggBioDF = NULL,
                       ggDepth = NULL,
                       ggDepthFGbox = NULL,
                       ggEffoFGbox = NULL,
                       ggEffoFGmap = NULL,
                       ggBioFGmat = NULL,
                       ggCutFGmap = NULL,
                       # ggIchFGlin = NULL,
                       ggSilFGlin = NULL,
                       ggBetaFGmap = NULL,
                       ggBetaFGbox = NULL,
                       ggProdFGmap = NULL,
                       ggProdFGbox = NULL,
                       ggMapFgFishery = NULL,
                       ggMapFgSurvey = NULL,
                       gooMap = NULL,
                       gooMapPlot = NULL,
                       gooGrid = NULL,
                       gooBbox = NULL,
                       sampColScale = NULL,
                       plotRange = NULL,
                       initialize = function(grid_path = "") {
                         if (grid_path != "") {
                           setGridPath(grid_path)
                           setGridName()
                           loadGridShp()
                           setNumCell()
                           createPolySet()
                           fortifyGridShp()
                           setGridCenter()
                           createGridBbox()
                         }
                       },
                       exportEnv = function() {
                         envOut <- list()
                         envOut$gridName <- gridName
                         envOut$gridShp <- gridShp
                         envOut$nCells <- nCells
                         envOut$gridPolySet <- gridPolySet
                         envOut$gridFortify <- gridFortify
                         envOut$griCent <- griCent
                         envOut$gridBbox <- gridBbox
                         envOut$gridBboxExt <- gridBboxExt
                         envOut$gridBboxSP <- gridBboxSP
                         envOut$gooMap <- gooMap
                         envOut$gooMapPlot <- gooMapPlot
                         envOut$plotRange <- plotRange
                         envOut$gooGrid <- gooGrid
                         envOut$gooBbox <- gooBbox
                         envOut$gridBathy <- gridBathy
                         envOut$centDept <- centDept
                         envOut$ggDepth <- ggDepth
                         envOut$bioDF <- bioDF
                         envOut$ggBioDF <- ggBioDF
                         return(envOut)
                       },
                       exportFG = function() {
                         fgOut <- list()
                         fgOut$cutFG <- cutFG
                         fgOut$rawInpu <- rawInpu
                         fgOut$clusMat <- clusMat
                         fgOut$indSil <- indSil
                         # fgOut$indCH <- indCH
                         return(fgOut)
                       },
                       importFG = function(fgList) {
                         cutFG <<- fgList$cutFG
                         rawInpu <<- fgList$rawInpu
                         clusMat <<- fgList$clusMat
                         indSil <<- fgList$indSil
                         # indCH <<- fgList$indCH
                       },
                       setAreaGrid = function() {
                         cat("\n\nComputing Total Area... ", sep = "")
                         clipDept <- as.SpatialGridDataFrame(gridBathy)
                         tmp_grid <- gridShp
                         proj4string(tmp_grid) <- proj4string(clipDept)
                         naFind <- which(is.na(over(clipDept, tmp_grid))[, 1])
                         if (length(naFind) > 0) clipDept[naFind] <- NA
                         clipDept <- as.bathy(clipDept)
                         areaGrid <<- get.area(clipDept, level.inf = -Inf, level.sup = 0)[[1]]
                         cat("Completed!", sep = "")
                       },
                       setAreaStrata = function(vectorStrata = c(0, 10, 100, 1000, Inf)) {
                         cat("\n\nComputing Area of Strata: ", sep = "")
                         clipDept <- as.SpatialGridDataFrame(gridBathy)
                         tmp_grid <- gridShp
                         proj4string(tmp_grid) <- proj4string(clipDept)
                         naFind <- which(is.na(over(clipDept, tmp_grid))[, 1])
                         if (length(naFind) > 0) clipDept[naFind] <- NA
                         clipDept <- as.bathy(clipDept)
                         strataList <- list()
                         for (stratum in 1:(length(vectorStrata) - 1)) {
                           stratum_ith <- paste(vectorStrata[stratum],
                                                "-",
                                                vectorStrata[stratum + 1],
                                                sep = ""
                           )
                           cat("\n\t", stratum_ith, "... ", sep = "")
                           strataList[[stratum_ith]] <- get.area(clipDept, level.inf = -vectorStrata[stratum + 1], level.sup = -vectorStrata[stratum])
                           cat("Done!", sep = "")
                         }
                         areaStrata <<- strataList
                         cat("\nCompleted!", sep = "")
                       },
                       setWeightStrata = function() {
                         cat("\n\nComputing Strata Weighted Area... ", sep = "")
                         weightStrata <<- unlist(lapply(areaStrata, "[[", 1)) / areaGrid
                         cat("Completed!")
                       },
                       loadHarbDbf = function(dbf_path) {
                         tmp_dbf <- read.dbf(file = dbf_path)
                         colnames(tmp_dbf) <- c("XCOORD", "YCOORD", "Name")
                         harbDbf <<- tmp_dbf
                       },
                       set_ggMapFgSurvey = function(rawSampCoo) {
                         if (length(cutFG) > 0) {
                           if (length(gridShp@polygons) == (cutFG + 1)) {
                             tmp_coo <- data.frame(coordinates(gridShp), cell_id = 1:length(gridShp))
                             colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                           } else {
                             tmp_coo <- cutResShpCent
                           }
                         } else {
                           tmp_coo <- cutResShpCent
                         }
                         
                         ggMapFgSurvey <<- suppressMessages(
                           gooMapPlot +
                             geom_polygon(
                               data = cutResShpFort,
                               aes(x = long, y = lat, group = group),
                               colour = "grey10", size = 0.1, alpha = 0.8
                             ) +
                             theme_tufte(base_size = 14, ticks = F) +
                             theme(
                               legend.position = "right",
                               panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                               axis.text.x = element_text(size = 5),
                               axis.title.x = element_text(size = 7),
                               axis.text.y = element_text(size = 5),
                               legend.key.size = unit(0.5, "cm"),
                               legend.text = element_text(size = 4),
                               legend.title = element_text(size = 5),
                               axis.title.y = element_text(size = 7)
                             ) +
                             geom_jitter(
                               data = rawSampCoo[, c("Lon", "Lat", "Species")],
                               aes(x = Lon, y = Lat, shape = Species, color = Species),
                               size = 1, width = 0.1, height = 0.1, alpha = 0.35
                             ) +
                             geom_point(
                               data = unique(rawSampCoo[, c("Lon", "Lat")]),
                               aes(x = Lon, y = Lat),
                               color = "darkolivegreen1", shape = 4, size = 0.5, alpha = 0.6
                             ) +
                             geom_text(
                               data = tmp_coo,
                               aes(label = FG, x = Lon, y = Lat),
                               size = 2, color = "grey72"
                             )
                         )
                       },
                       set_ggMapFgFishery = function(rawSampCoo) {
                         if (length(cutFG) > 0) {
                           if (length(gridShp@polygons) == (cutFG + 1)) {
                             tmp_coo <- data.frame(coordinates(gridShp), cell_id = 1:length(gridShp))
                             colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                           } else {
                             tmp_coo <- cutResShpCent
                           }
                         } else {
                           tmp_coo <- cutResShpCent
                         }
                         
                         ggMapFgFishery <<- suppressMessages(
                           gooMapPlot +
                             geom_polygon(
                               data = cutResShpFort,
                               aes(x = long, y = lat, group = group),
                               colour = "grey10", size = 0.1, alpha = 0.8
                             ) +
                             theme(
                               legend.position = "right",
                               axis.text.x = element_text(size = 5),
                               axis.title.x = element_text(size = 7),
                               axis.text.y = element_text(size = 5),
                               legend.key.size = unit(0.5, "cm"),
                               legend.text = element_text(size = 4),
                               legend.title = element_text(size = 5),
                               axis.title.y = element_text(size = 7)
                             ) +
                             geom_jitter(
                               data = rawSampCoo[, c("Lon", "Lat", "Species")],
                               aes(x = Lon, y = Lat, shape = Species, color = Species),
                               size = 1, width = 0.1, height = 0.1, alpha = 0.35
                             ) +
                             geom_point(
                               data = unique(rawSampCoo[, c("Lon", "Lat")]),
                               aes(x = Lon, y = Lat),
                               color = "darkolivegreen1", shape = 4, size = 0.5, alpha = 0.6
                             ) +
                             geom_text(
                               data = tmp_coo,
                               aes(label = FG, x = Lon, y = Lat),
                               size = 2, color = "grey72"
                             )
                         )
                       },
                       createGridBbox = function() {
                         gridBbox <<- bbox(gridShp)
                         lon_range <- extendrange(range(gridBbox[1, ], na.rm = TRUE), f = 0.1)
                         lat_range <- extendrange(range(gridBbox[2, ], na.rm = TRUE), f = 0.1)
                         gridBboxExt <<- c(
                           left = lon_range[1],
                           bottom = lat_range[1],
                           right = lon_range[2],
                           top = lat_range[2]
                         )
                         polyext <- Polygon(cbind(
                           c(gridBboxExt[["left"]], gridBboxExt[["left"]], gridBboxExt[["right"]], gridBboxExt[["right"]], gridBboxExt[["left"]]),
                           c(gridBboxExt[["bottom"]], gridBboxExt[["top"]], gridBboxExt[["top"]], gridBboxExt[["bottom"]], gridBboxExt[["bottom"]])
                         ))
                         
                         polypolyext <- Polygons(list(polyext), "s1")
                         
                         gridBboxSP <<- SpatialPolygons(list(polypolyext))
                       },
                       getGooMap = function() {
                         gooMap <<- ggplot() +
                           borders(
                             fill = "black", colour = "black",
                             xlim = gridBboxExt[c(1, 3)],
                             ylim = gridBboxExt[c(2, 4)]
                           ) +
                           coord_map("bonne",
                                     xlim = gridBboxExt[c(1, 3)],
                                     ylim = gridBboxExt[c(2, 4)],
                                     lat0 = mean(gridBboxExt[c(2, 4)])
                           )
                         setGooPlot()
                         setPlotRange()
                       },
                       setGooPlot = function() {
                         gooMapPlot <<- gooMap + xlab("Longitude") + ylab("Latitude")
                       },
                       setPlotRange = function() {
                         plotRange <<- data.frame(
                           xmin = gridBboxExt[1],
                           xmax = gridBboxExt[3],
                           ymin = gridBboxExt[2],
                           ymax = gridBboxExt[4]
                         )
                       },
                       setGooGrid = function() {
                         gooGrid <<- suppressMessages(gooMapPlot + geom_polygon(aes(x = long, y = lat, group = group, fill = ""),
                                                                                size = 0.1,
                                                                                color = "gainsboro", data = gridFortify, alpha = 0.5
                         ) +
                           scale_fill_manual("Case Study Cells", values = "grey50") +
                           xlab("Longitude") + ylab("Latitude") +
                           ggtitle("Grid") +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "bottom",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10)
                           ))
                       },
                       setGooEnv = function() {
                         if (is.null(gooGrid)) {
                           tmp_grid <- ggplot() + geom_blank() + ggtitle("No grid Loaded")
                         } else {
                           tmp_grid <- gooGrid
                         }
                         if (is.null(ggDepth)) {
                           tmp_dept <- ggplot() + geom_blank() + ggtitle("No depth Loaded")
                         } else {
                           tmp_dept <- ggDepth
                         }
                         if (is.null(ggBioDF)) {
                           tmp_bioc <- ggplot() + geom_blank() + ggtitle("No seabed Loaded")
                         } else {
                           tmp_bioc <- ggBioDF
                         }
                         
                         gooEnv <<- suppressWarnings(grid.arrange(tmp_grid,
                                                                  tmp_dept,
                                                                  tmp_bioc,
                                                                  layout_matrix = matrix(1:3, 1, 3)
                         ))
                       },
                       plotGooEnv = function() {
                         suppressWarnings(grid.draw(gooEnv))
                       },
                       plotGooGrid = function() {
                         suppressWarnings(print(gooGrid))
                       },
                       plotGooGridData = function(grid_data) {
                         gooMapPlot + geom_polygon(aes(x = X, y = Y, group = PID),
                                                   fill = "grey", size = 0.2,
                                                   color = "gainsboro", data = grid_data, alpha = 0.5
                         )
                       },
                       setSampColScale = function(fac_col) {
                         myColors <- brewer.pal(length(fac_col), "Set1")
                         names(myColors) <- fac_col
                         sampColScale <<- scale_colour_manual(name = "Species", values = myColors)
                       },
                       plotGooSpeSur = function(poi_data) {
                         temp_pos <- suppressMessages(gooGrid + geom_jitter(
                           data = poi_data,
                           aes(x = Lon, y = Lat, shape = Species, color = Species),
                           width = 0.05, height = 0.05, alpha = 0.95
                         ) + sampColScale)
                         suppressWarnings(print(temp_pos))
                       },
                       plotGooSpeFis = function(poi_data) {
                         temp_pos <- suppressMessages(gooGrid + geom_jitter(
                           data = poi_data,
                           aes(x = Lon, y = Lat, shape = Species, color = Species),
                           width = 0.05, height = 0.05, alpha = 0.95
                         ))
                         # temp_pos <- gooGrid + geom_jitter(data = poi_data,
                         #                             aes(x = Lon, y = Lat, shape = Species, color = Species),
                         #                             width = 0.05, height = 0.05, alpha = 0.95)
                         # print(temp_pos)
                         suppressWarnings(print(temp_pos))
                       },
                       setGooBbox = function() {
                         text_x <- mean(gridBboxExt[c(1, 3)])
                         text_y <- mean(gridBboxExt[c(2, 4)])
                         gooBbox <<- suppressMessages(gooGrid + geom_rect(
                           data = plotRange, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
                           color = "firebrick",
                           fill = alpha("red", 0.2),
                           inherit.aes = FALSE
                         ) +
                           annotate("label",
                                    x = text_x, y = text_y,
                                    label = "Bounding\nBox", family = "serif", fontface = "italic",
                                    colour = "firebrick", size = 6, fill = "grey80"
                           ) +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "none",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10)
                           ))
                       },
                       plotGooBbox = function() {
                         suppressWarnings(print(gooBbox))
                       },
                       setGridPath = function(path2grid) {
                         gridPath <<- path2grid
                       },
                       setGridName = function() {
                         tmp_name <- unlist(strsplit(gridPath, "/"))
                         tmp_name <- tmp_name[length(tmp_name)]
                         gridName <<- substr(tmp_name, 1, nchar(tmp_name) - 4)
                       },
                       loadGridShp = function() {
                         gridShp <<- readOGR(gridPath)
                       },
                       setBioPath = function(path2bio) {
                         bioPath <<- path2bio
                       },
                       setBioName = function() {
                         tmp_name <- unlist(strsplit(bioPath, "/"))
                         tmp_name <- tmp_name[length(tmp_name)]
                         bioName <<- substr(tmp_name, 1, nchar(tmp_name) - 4)
                       },
                       loadBioShp = function() {
                         bioShp <<- readOGR(bioPath)
                       },
                       addBioShp = function(bio_path) {
                         setBioPath(bio_path)
                         setBioName()
                         loadBioShp()
                       },
                       loadBioDF = function(bio_path) {
                         bioDF <<- readRDS(bio_path)
                         setGgBioDF()
                       },
                       plotBioDF = function() {
                         def.par <- par(no.readonly = TRUE)
                         par(mar = c(2.5, 2.5, 3, 1))
                         layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(6, 1))
                         
                         vec_bio <- apply(bioDF, 1, function(x) which(x == 1))
                         color_clas <- rainbow(max(vec_bio))
                         plotSamMap(title = "Biocenosis", celCol = color_clas[vec_bio])
                         
                         par(mar = c(5, 1, 1, 1))
                         par(mar = c(1, 0, 1, 0))
                         
                         plot(NULL, xlim = c(0, 1), ylim = c(0, 1), bty = "n", axes = FALSE, ann = FALSE)
                         legend(x = 0, y = 0.5, legend = colnames(bioDF), fill = color_clas, bty = "n")
                         par(def.par)
                       },
                       setGgBioDF = function() {
                         cell_bed <- apply(bioDF, 1, which.max)
                         
                         zero_cell <- which(apply(bioDF, 1, sum) == 0)
                         
                         tmp_dat <- colnames(bioDF)[cell_bed]
                         
                         if (length(zero_cell) > 0) tmp_dat[zero_cell] <- "No Data"
                         
                         color_clas <- rainbow(max(cell_bed))
                         names(tmp_dat) <- 1:length(tmp_dat)
                         all_cell <- merge(
                           x = gridPolySet$PID,
                           data.frame(
                             x = as.numeric(names(tmp_dat)), y = tmp_dat,
                             stringsAsFactors = FALSE
                           ), all = TRUE
                         )
                         all_cell[is.na(all_cell)] <- 0
                         grid_data <- cbind(gridPolySet, Seabed = all_cell[, 2])
                         ggBioDF <<- suppressMessages(gooMapPlot +
                                                        geom_polygon(aes(x = X, y = Y, group = PID, fill = Seabed),
                                                                     size = 0.2,
                                                                     data = grid_data, alpha = 0.8
                                                        ) +
                                                        xlab("Longitude") + ylab("Latitude") +
                                                        ggtitle("Seabed") +
                                                        theme_tufte(base_size = 14, ticks = T) +
                                                        theme(
                                                          legend.position = "bottom",
                                                          axis.text.x = element_text(size = 8),
                                                          axis.title.x = element_text(size = 10),
                                                          panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                                          axis.text.y = element_text(size = 8),
                                                          axis.title.y = element_text(size = 10),
                                                          legend.text = element_text(size = 8),
                                                          legend.title = element_text(size = 10)
                                                        ))
                       },
                       ggplotBioDF = function() {
                         suppressWarnings(print(ggBioDF))
                       },
                       createPolySet = function() {
                         gridPolySet <<- as.data.frame(SpatialPolygons2PolySet(gridShp))
                       },
                       fortifyGridShp = function() {
                         gridFortify <<- fortify(gridShp)
                       },
                       setNumCell = function() {
                         nCells <<- length(gridShp@polygons)
                       },
                       setGridCenter = function() {
                         griCent <<- coordinates(gridShp) # or gCentroid from rgeos
                       },
                       getGridBath = function() {
                         lon_ran <- extendrange(griCent[, 1], f = 0.05)
                         lat_ran <- extendrange(griCent[, 2], f = 0.05)
                         gridBathy <<- getNOAA.bathy(
                           lon1 = lon_ran[1],
                           lon2 = lon_ran[2],
                           lat1 = lat_ran[1],
                           lat2 = lat_ran[2], resolution = 1
                         )
                         getCentDept()
                         setGgDepth()
                       },
                       saveGridBath = function(bathy_path) {
                         saveRDS(object = gridBathy, file = bathy_path)
                       },
                       loadGridBath = function(bathy_path) {
                         gridBathy <<- readRDS(bathy_path)
                         getCentDept()
                         setGgDepth()
                       },
                       getCentDept = function() {
                         centDept <<- get.depth(gridBathy, x = griCent[, 1], y = griCent[, 2], locator = FALSE)
                       },
                       setGgDepth = function(isoLine = c(-200, -1000)) {
                         f_bathy <- fortify.bathy(gridBathy)
                         f_bathy$z[f_bathy$z > 0] <- 0
                         colnames(f_bathy) <- c("lon", "lat", "Depth")
                         ggDepth <<- suppressMessages(gooMapPlot +
                                                        geom_contour(aes(x = lon, y = lat, z = Depth, colour = factor(..level..)),
                                                                     data = f_bathy,
                                                                     linetype = "solid", size = 0.35,
                                                                     breaks = isoLine, alpha = 1
                                                        ) +
                                                        scale_colour_brewer(palette = "Accent", name = "Isobath") +
                                                        guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1))) +
                                                        xlab("Longitude") + ylab("Latitude") +
                                                        ggtitle("Depth") +
                                                        theme_tufte(base_size = 14, ticks = T) +
                                                        theme(
                                                          legend.position = "bottom",
                                                          axis.text.x = element_text(size = 8),
                                                          axis.title.x = element_text(size = 10),
                                                          panel.grid = element_line(size = 0.5, linetype = 2, colour = "grey20"),
                                                          axis.text.y = element_text(size = 8),
                                                          axis.title.y = element_text(size = 10),
                                                          legend.text = element_text(size = 8),
                                                          legend.title = element_text(size = 10)
                                                        ))
                       },
                       ggplotGridBathy = function() {
                         suppressWarnings(print(ggDepth))
                       },
                       plotSamMap = function(title = "", celCol = NULL) {
                         par(mar = c(3, 3, 1.5, 0.5))
                         plotPolys(gridPolySet,
                                   main = title, ylab = "", xlab = "", col = celCol, axes = TRUE,
                                   xlim = extendrange(c(gridShp@bbox[1, 1], gridShp@bbox[1, 2]), f = 0.05),
                                   ylim = extendrange(c(gridShp@bbox[2, 1], gridShp@bbox[2, 2]), f = 0.05)
                         )
                         mtext("Longitude", side = 1, line = 2.1, cex = 1.5)
                         mtext("Latitude", side = 2, line = 2, cex = 1.5)
                         map("worldHires", fill = T, col = "gainsboro", add = TRUE)
                         # map.axes(cex.axis=0.8)
                         map.scale(cex = 0.75, ratio = FALSE)
                       },
                       plotCoho = function(abbs) {
                         x_rang <- extendrange(c(gridShp@bbox[1, 1], gridShp@bbox[1, 2]), f = 0.05)
                         y_rang <- extendrange(c(gridShp@bbox[2, 1], gridShp@bbox[2, 2]), f = 0.07)
                         plotPolys(gridPolySet,
                                   ylab = "", xlab = "",
                                   xlim = x_rang,
                                   ylim = y_rang
                         )
                         mtext("Longitude", side = 1, line = 2.1, cex = 1.1)
                         mtext("Latitude", side = 2, line = 2, cex = 1.1)
                         smoothScatter(griCent[rep.int(1:nrow(griCent), abbs + 1), ],
                                       colramp = colorRampPalette(c("white", "white", "blue", "purple"), bias = 0.77, alpha = 0.5),
                                       add = TRUE, nbin = 250, xlab = NULL, ylab = NULL, nrpoints = 0
                         )
                         plot(gridShp,
                              ylab = "", xlab = "",
                              add = TRUE, lwd = 0.25
                         )
                         map("worldHires",
                             fill = T, col = "gainsboro",
                             xlim = x_rang,
                             ylim = y_rang, add = TRUE
                         )
                         map.scale(x = min(abs(x_rang)) + (diff(x_rang) / 10), y = min(abs(y_rang)) + (diff(y_rang) / 7), cex = 0.65, ratio = FALSE)
                       },
                       setClusInpu = function(howData = rep(1, 3)) {
                         tmp_lst <- list()
                         cat("\n   -   Seabed")
                         tmp_lst <- c(tmp_lst, list(rawInpu[[1]] * howData[1]))
                         cat("\t\t-   Set!")
                         cat("\n   -   Effort")
                         tmp_lst <- c(tmp_lst, list(vegan::decostand(log10(1 + rawInpu[[2]]), method = "range", MARGIN = 2) * howData[2]))
                         cat("\t\t-   Set!")
                         cat("\n   -   Depth")
                         tmp_inpu <- -rawInpu[[3]]
                         tmp_inpu[tmp_inpu < 0] <- 0
                         tmp_lst <- c(tmp_lst, Depth = list(vegan::decostand(log10(1 + tmp_inpu), method = "range", MARGIN = 2) * howData[3]))
                         cat("\t\t-   Set!\n")
                         clusInpu <<- do.call(cbind, tmp_lst)
                       },
                       calcFishGrou = function(numCuts = 50,
                                               minsize = 10,
                                               maxsize = 100,
                                               modeska = "S",
                                               skater_method,
                                               nei_queen = TRUE) {
                         set.seed(123)
                         # Build the neighboorhod list
                         Grid.bh <- gridShp[1]
                         ##  Construct neighbours list from polygon list
                         bh.nb <- poly2nb(Grid.bh, queen = nei_queen)
                         bh.mat <- cbind(rep(1:length(bh.nb), lapply(bh.nb, length)), unlist(bh.nb))
                         # Compute the basic objects for clustering
                         ##  Cost of each edge as the distance between nodes
                         lcosts <- nbcosts(bh.nb, clusInpu)
                         ##  Spatial weights for neighbours lists
                         nb.w <- nb2listw(bh.nb, lcosts, style = modeska)
                         ##  Find the minimal spanning tree
                         mst.bh <- mstree(nb.w, ini = 1)
                         clusMat <<- matrix(NA, nCells, numCuts)
                         cat("\nClustering", sep = "")
                         res1 <- skater(
                           edges = mst.bh[, 1:2],
                           data = clusInpu,
                           ncuts = 1,
                           crit = minsize,
                           method = skater_method
                         )
                         clusMat[, 1] <<- res1$groups
                         
                         # Perform the first CC (without removing spurious clusters)
                         for (nCuts in 2:numCuts) {
                           cat(".", sep = "")
                           ##  Spatial 'K'luster Analysis by Tree Edge Removal
                           #                            res1 <- skater(mst.bh[,1:2], cells_data, ncuts = nCuts, minsize,
                           #                                           method = skater_method)
                           
                           if (nCuts > ceiling(nrow(clusInpu) / maxsize)) {
                             res1 <- skater(res1, clusInpu,
                                            ncuts = 1,
                                            crit = c(minsize, maxsize),
                                            method = skater_method
                             )
                           } else {
                             res1 <- skater(res1, clusInpu,
                                            ncuts = 1,
                                            crit = minsize,
                                            method = skater_method
                             )
                           }
                           clusMat[, nCuts] <<- res1$groups
                         }
                         cat(" Done!", sep = "")
                         
                         indSil <<- numeric(ncol(clusMat))
                         # indCH <<- numeric(ncol(clusMat))
                         for (i in 1:ncol(clusMat)) {
                           ##  Compute silhouette information according to a given clustering in k clusters
                           indSil[i] <<- summary(silhouette(clusMat[, i], dist(clusInpu,
                                                                               method = skater_method
                           )))$avg.width
                           ##  Compute CH index for a given partition of a data set
                           # indCH[i] <<- get_CH(as.matrix(clusInpu), clusMat[, i])
                         }
                       },
                       setCutResult = function(ind_clu) {
                         # cutResult <<- data.frame(clusInpu, FG = as.factor(clusMat[,ind_clu]))
                         cutResult <<- data.frame(do.call(cbind, rawInpu), FG = as.factor(clusMat[, ind_clu]))
                         cutResEffo <<- data.frame(
                           Effort = apply(cutResult[, grep("Year", colnames(cutResult))], 1, mean),
                           Cluster = cutResult[, ncol(cutResult)]
                         )
                         cutResShp <<- unionSpatialPolygons(gridShp, IDs = clusMat[, ind_clu])
                         cutResShp@plotOrder <<- 1:ind_clu
                         plot.new()
                         cutResShpCent <<- as.data.frame(polygonsLabel(cutResShp, method = "inpolygon", doPlot = FALSE))
                         cutResShpCent$id <<- names(cutResShp)
                         names(cutResShpCent) <<- c("Lon", "Lat", "FG")
                         cutResShpFort <<- fortify(cutResShp)
                         cutResShpFort$FG <<- as.factor(cutResShpFort$id)
                         setDepthFGbox()
                         setEffoFGbox()
                         setEffoFGmap()
                         setBioFGmat()
                         setCutFGmap()
                         # setIchFGlin(numCut = ind_clu)
                         setSilFGlin(numCut = ind_clu)
                       },
                       setDepthFGbox = function() {
                         ggDepthFGbox <<- suppressMessages(ggplot(cutResult, aes(x = FG, y = Depth, group = FG)) +
                                                             geom_boxplot(color = "grey23") +
                                                             coord_flip() +
                                                             xlab("Fishing Ground") +
                                                             theme_tufte(base_size = 14, ticks = T) +
                                                             ylim(NA, 0) +
                                                             theme(
                                                               legend.position = "none",
                                                               axis.text.x = element_text(size = 8),
                                                               axis.title.x = element_text(size = 10),
                                                               panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                               axis.text.y = element_text(size = 8),
                                                               axis.title.y = element_text(size = 10),
                                                               legend.text = element_text(size = 8),
                                                               legend.title = element_text(size = 10)
                                                             ))
                       },
                       setEffoFGbox = function() {
                         ggEffoFGbox <<- suppressMessages(ggplot(cutResEffo, aes(x = Cluster, y = Effort, group = Cluster)) +
                                                            geom_boxplot(color = "grey23") +
                                                            coord_flip() +
                                                            xlab("Fishing Ground") +
                                                            theme_tufte(base_size = 14, ticks = T) +
                                                            theme(
                                                              legend.position = "none",
                                                              axis.text.x = element_text(size = 8),
                                                              axis.title.x = element_text(size = 10),
                                                              panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                              axis.text.y = element_text(size = 8),
                                                              axis.title.y = element_text(size = 10),
                                                              legend.text = element_text(size = 8),
                                                              legend.title = element_text(size = 10)
                                                            ))
                       },
                       setEffoFGmap = function() {
                         agg_eff <- aggregate(formula = Effort ~ Cluster, data = cutResEffo, FUN = mean)
                         all_cell <- merge(
                           x = cutResShpFort$id,
                           data.frame(x = agg_eff$Cluster, y = agg_eff$Effort), all = TRUE
                         )
                         all_cell[is.na(all_cell)] <- 0
                         grid_data <- cbind(cutResShpFort, Hours = all_cell[, 2])
                         
                         if (length(cutFG) > 0) {
                           if (length(gridShp@polygons) == (cutFG + 1)) {
                             tmp_coo <- data.frame(coordinates(gridShp), cell_id = 1:length(gridShp))
                             colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                           } else {
                             tmp_coo <- cutResShpCent
                           }
                         } else {
                           tmp_coo <- cutResShpCent
                         }
                         
                         ggEffoFGmap <<- suppressMessages(gooMapPlot + geom_polygon(aes(x = long, y = lat, group = group, fill = Hours),
                                                                                    colour = "black", size = 0.1,
                                                                                    data = grid_data, alpha = 0.8
                         ) +
                           scale_fill_gradient(low = "Yellow", high = "coral", trans = "sqrt") +
                           geom_text(aes(label = FG, x = Lon, y = Lat),
                                     data = tmp_coo, size = 2
                           ) +
                           ggtitle("Average Effort Intensity") +
                           theme_tufte(base_size = 14, ticks = T) +
                           theme(
                             legend.position = "none",
                             axis.text.x = element_text(size = 8),
                             axis.title.x = element_text(size = 10),
                             panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                             axis.text.y = element_text(size = 8),
                             axis.title.y = element_text(size = 10)
                           ))
                       },
                       setBioFGmat = function() {
                         tmp_bio <- data.frame(FG = cutResult$FG, bioDF)
                         bio2plot <- melt(tmp_bio, id.vars = "FG", variable.name = "Substrate")
                         bio2plot <- bio2plot[bio2plot$value == 1, 1:2]
                         ggBioFGmat <<- suppressMessages(ggplot(bio2plot, aes(x = FG, y = Substrate, fill = Substrate)) +
                                                           geom_tile() +
                                                           coord_flip() +
                                                           annotate("text",
                                                                    colour = "grey30", y = 1:length(levels(bio2plot$Substrate)),
                                                                    x = rep(4, length(levels(bio2plot$Substrate))),
                                                                    label = levels(bio2plot$Substrate),
                                                                    angle = rep(90, length(levels(bio2plot$Substrate)))
                                                           ) +
                                                           theme_tufte(base_size = 14, ticks = T) +
                                                           xlab("Fishing Ground") +
                                                           theme(
                                                             legend.position = "none",
                                                             axis.text.x = element_text(size = 8),
                                                             axis.title.x = element_text(size = 10, colour = "white"),
                                                             panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                             axis.text.y = element_text(size = 8),
                                                             axis.title.y = element_text(size = 10),
                                                             legend.text = element_text(size = 8),
                                                             legend.title = element_text(size = 10)
                                                           ))
                         # }
                       },
                       setCutFGmap = function() {
                         if (length(cutFG) > 0) {
                           if (length(gridShp@polygons) == (cutFG + 1)) {
                             tmp_coo <- data.frame(coordinates(gridShp), cell_id = 1:length(gridShp))
                             colnames(tmp_coo) <- c("Lon", "Lat", "FG")
                           } else {
                             tmp_coo <- cutResShpCent
                           }
                         } else {
                           tmp_coo <- cutResShpCent
                         }
                         
                         ggCutFGmap <<- suppressMessages(gooMapPlot +
                                                           geom_polygon(aes(x = long, y = lat, group = group, fill = FG),
                                                                        colour = "black", size = 0.1,
                                                                        data = cutResShpFort, alpha = 0.8
                                                           ) +
                                                           geom_text(aes(label = FG, x = Lon, y = Lat),
                                                                     data = tmp_coo, size = 2
                                                           ) +
                                                           scale_fill_manual(values = colorRampPalette(brewer.pal(8, "Accent"))(length(unique(cutResShpFort$FG)))) +
                                                           ggtitle("Regions") +
                                                           theme_tufte(base_size = 14, ticks = T) +
                                                           theme(
                                                             legend.position = "none",
                                                             axis.text.x = element_text(size = 8),
                                                             axis.title.x = element_text(size = 10),
                                                             panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                             axis.text.y = element_text(size = 8),
                                                             axis.title.y = element_text(size = 10)
                                                           ))
                       },
                       # setIchFGlin = function(numCut) {
                       #   ch_df <- data.frame(
                       #     cut = 1:length(indCH),
                       #     CH_index = indCH
                       #   )
                       #   ggIchFGlin <<- suppressMessages(ggplot(ch_df, aes(x = cut, y = CH_index)) +
                       #     geom_line() +
                       #     geom_vline(aes(xintercept = numCut),
                       #       linetype = "dashed", size = 0.5, colour = "red"
                       #     ) +
                       #     ggtitle("Calinski-Harabasz") +
                       #     ylab("Index") +
                       #     theme_tufte(base_size = 14, ticks = T) +
                       #     theme(
                       #       legend.position = "bottom",
                       #       axis.text.x = element_text(size = 8),
                       #       axis.title.x = element_blank(),
                       #       panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                       #       axis.text.y = element_text(size = 8),
                       #       axis.title.y = element_text(size = 10)
                       #     ))
                       # },
                       setSilFGlin = function(numCut) {
                         sil_df <- data.frame(
                           cut = 1:length(indSil),
                           Sil_index = indSil
                         )
                         ggSilFGlin <<- suppressMessages(ggplot(sil_df, aes(x = cut, y = Sil_index)) +
                                                           geom_line() +
                                                           geom_vline(aes(xintercept = numCut),
                                                                      linetype = "dashed", size = 0.5, colour = "red"
                                                           ) +
                                                           ggtitle("Silhouette") +
                                                           ylab("Index") +
                                                           theme_tufte(base_size = 14, ticks = T) +
                                                           theme(
                                                             legend.position = "bottom",
                                                             axis.text.x = element_text(size = 8),
                                                             axis.title.x = element_blank(),
                                                             panel.grid = element_line(size = 0.05, linetype = 2, colour = "grey20"),
                                                             axis.text.y = element_text(size = 8),
                                                             axis.title.y = element_text(size = 10)
                                                           ))
                       }
                     )
)

Try the smartR package in your browser

Any scripts or data that you put into this service are public.

smartR documentation built on July 2, 2020, 2:25 a.m.