R/iSCAM2OM.R

Defines functions iSCAM2Data LinInt iSCAMcomps iSCAMinds strip.static.params strip.areas.groups calc.mcmc extract.area.sex.matrices extract.group.matrices read.mcmc read.par.file read.projection.file read.control.file read.data.file read.report.file fetch.file.names load.iscam.files iSCAM2OM

Documented in fetch.file.names iSCAM2Data iSCAM2OM iSCAMcomps iSCAMinds load.iscam.files read.control.file read.data.file read.mcmc read.par.file read.projection.file read.report.file

# === OM specification using iSCAM stock assessment ====================


#' Reads MLE estimates from iSCAM file structure into an operating model
#'
#' @description A function that uses the file location of a fitted iSCAM
#' model including input files to population the various slots of an
#' operating model parameter estimates. iSCAM2OM relies on several
#' functions written by Chris Grandin (DFO PBS).
#' @param iSCAMdir A folder with iSCAM input and output files in it
#' @param nsim The number of simulations to take for parameters with
#' uncertainty (for OM@cpars custom parameters)
#' @param proyears The number of MSE projection years
#' @param mcmc Whether to use mcmc samples to create custom parameters cpars
#' @param Name The name of the operating model
#' @param Source Reference to assessment documentation e.g. a url
#' @param length_timestep How long is a model time step in years
#' (e.g. a quarterly model is 0.25, a monthly model 1/12)
#' @param nyr_par_mu integer, the number of recent years to estimate vulnerability over for future projections
#' @param LowerTri Integer. The number of recent years for which model estimates of recruitment are ignored (not reliably estimated by the assessment)
#' @param Author Who did the assessment
#' @param report logical should a numbers at age reconstruction plot be produced?
#' @param silent logical should progress reporting be printed to the console?
#' @author T. Carruthers
#' @importFrom grDevices dev.off gray jpeg png
#' @importFrom coda mcmc
#' @importFrom graphics arrows contour
#' @importFrom stats acf aggregate qnorm window
#' @export iSCAM2OM
iSCAM2OM<-function(iSCAMdir, nsim=48, proyears=50, mcmc=F, Name=NULL, Source="No source provided",
                 length_timestep=1, nyr_par_mu=2,  Author="No author provided", report=F, silent=F){

  nseas<-1/length_timestep # defaults to an annual model
  message("-- Using function of Chris Grandin (DFO PBS) to extract data from iSCAM file structure --")

  replist<-load.iscam.files(iSCAMdir)

  message("-- End of iSCAM extraction operations --")

  # get dimensions
  nyears<-(replist$dat$end.yr-replist$dat$start.yr+1)
  maxage<-replist$dat$end.age
  sage<-replist$dat$start.age

  # filling dimensions
  aind<-sage:maxage # for filling all quantities (that do not include age zero)
  nafill<-length(aind)

  ageArray<-array(rep(1:maxage,each=nsim),c(nsim,maxage,nyears))

  # make matrices
  naa<-Mataa<-laa<-array(NA,c(nsim,maxage,nyears))
  faa<-Maa<-waa<-array(0,c(nsim,maxage,nyears))

  # growth parameters
  Linf=replist$dat$linf[1]
  K=replist$dat$k[1]
  t0<-replist$dat$to[1]

  #ageM<-replist$dat$age.at.50.mat
  #ageMsd<-replist$dat$sd.at.50.mat

  # rip F, N and wt matrices
  faat<-t(replist$mpd$F) # maxage-1 x nyears-1
  naat<-t(replist$mpd$N) # maxage-1 x nyears
  naat<-rbind(c(naat[1,2:(nyears+1)],NA)*exp(replist$mpd$m),naat)
  naat<-naat[,1:nyears]
  waat<-t(replist$mpd$d3_wt_mat)[,1:nyears] # maxage-1 x nyears

  # Numbers
  naa<-array(rep(naat,each=nsim),c(nsim,maxage,nyears)) # !

  # Mat at age
  Maa[,aind,]<-rep(t(replist$mpd$M),each=nsim)   # !

  # Weight at age
  waa[,aind,]<-array(rep(waat,each=nsim),c(nsim,nafill,nyears)) # !

  # Fishing mortality rate  for F = 0 years
  Vtemp<-apply(faat,1,mean)
  Vtemp<-Vtemp/max(Vtemp)*1E-5 #fill with a very low typical vulnerability
  tofill<-apply(faat,2,function(x)all(x==0))
  faat[,tofill]<-rep(Vtemp,sum(tofill))

  faa[,aind,]<-array(rep(faat,each=nsim),c(nsim,nafill,nyears))
  #faa[,aind,1]<-rep(Vtemp,each=nsim)

  #Mataa<- 1/(1 + exp((ageM - ageArray)/ageMsd))
  Mataat<-c(0,replist$mpd$ma)
  Mataa[]<-rep(Mataat,each=nsim)
  laa<-Linf*(1-exp(-K*(ageArray-t0)))
  h<-rep(replist$mpd$steepness,nsim)

  # make the OM
  OM<-VPA2OM(Name="A fishery made by VPA2OM",
             proyears=50, interval=2, CurrentYr=2019,
             h=h,
             Obs = DLMtool::Imprecise_Unbiased, Imp=DLMtool::Perfect_Imp,
             naa, faa, waa, Mataa, Maa, laa,
             nyr_par_mu = nyr_par_mu, LowerTri=1,
             recind=2, plusgroup=TRUE, altinit=2, fixq1=TRUE,
             report=report, silent=FALSE)

  # Observation model parameters ==============================================================================


  dat<-iSCAM2Data

  # --- mcmc functionality ------------------------------------

  if(mcmc){

    message("Attempting to read mcmc file to assign posterior samples to custom parameters")

    model.dir=iSCAMdir

    if(!file.exists(model.dir))stop(paste("Could not find the mcmc subfolder:",model.dir))

    tmp<-read.mcmc(model.dir)
    nmcmc<-nrow(tmp$params)

    if(nsim<nmcmc){
      samp<-sample(1:nmcmc,size=nsim)
    }else{
      message("You requested a greater number of simulations than the number of mcmc samples that are available - sampling with replacement")
      samp<-sample(1:nmcmc,size=nsim,replace=T)
    }

    #@nyears<-nyears<-ncol(tmp$sbt[[1]])
    #M<-tmp$params$m_gs1[samp]
    #OM@M<-quantile(M,c(0.05,0.95))
    hs<-tmp$params$h_gr1[samp]
    #OM@h<-quantile(hs,c(0.05,0.95))
    R0<-(1E6)*tmp$params$ro_gr1[samp]
    #OM@R0<-quantile(R0, c(0.05,0.95))

    #ssb_r <-tmp$params$sbo[samp]/R0
    D<-tmp$sbt[[1]][samp,nyears]/tmp$params$sbo[samp]#*ssb_r
    #OM@D<-quantile(D,c(0.05,0.95))

    #recdevs<-tmp$rdev[[1]][samp,]

    #procsd<-apply(recdevs,1,sd)
    #OM@Perr<-quantile(procsd,c(0.05,0.95))

    #nrecs<-ncol(recdevs)
    #AC<-apply(recdevs,1,function(x)acf(x)$acf[2,1,1])
    #OM@AC<-quantile(AC,c(0.05,0.95))

    #procmu <- -0.5 * (procsd)^2  # adjusted log normal mean

    Perr<-matrix(rnorm(nsim*(maxage+nyears+proyears-1),rep(procmu,maxage+nyears+proyears-1),rep(procsd,maxage+nyears+proyears-1)),nrow=nsim)
    Perr[,(maxage-1)+(1:nrecs)]<-as.matrix(recdevs) # generate a bunch of simulations with uncertainty

    for (y in c(2:(maxage-1),(-1:(proyears-1))+(maxage+nyears))) Perr[, y] <- AC * Perr[, y - 1] +   Perr[, y] * (1 - AC * AC)^0.5
    Perr<-exp(Perr)

    agesind<-ncol(replist$mpd$d3_wt_mat)
    surv<-exp(cumsum(c(rep(-replist$mpd$m,maxage))))[maxage-((agesind-1):0)]

    agearr<-maxage-((agesind-1):0)
    survarr<-exp(-array(rep(M,agesind)*rep(agearr,each=nsim),c(nsim,agesind)))
    N0vec<-array(rep(R0,agesind)*rep(surv,each=nsim),c(nsim,agesind))
    N0pred<-array(rep(replist$mpd$N[1,],each=nsim)*1E6,dim(N0vec))

    RecDev0<-log(N0pred/N0vec)
    Perr[,ncol(N0vec):1]<-RecDev0

    nfleet<-length(tmp$ft[[1]])
    FM<-tmp$ft[[1]][[1]][samp,1:nyears]
    for(ff in 2:nfleet)FM<-FM+tmp$ft[[1]][[ff]][samp,1:nyears]
    Find<-as.matrix(FM/apply(FM,1,mean))

    OM@cpars<-list(V=V,Perr=Perr,Wt_age=Wt_age2,K=K,Linf=Linf,hs=hs,Find=Find,D=D,M=M,R0=R0,AC=AC)

  }

  OM

}


#' Reads iSCAM files into a hierarchical R list object
#'
#' @description A function for reading iSCAM input and output files
#' into R
#' @param model.dir An iSCAM directory
#' @param burnin The initial mcmc samples to be discarded
#' @param thin The degree of chain thinning 1 in every thin
#' iterations is kept
#' @param verbose Should detailed outputs be provided.
#' @author Chris Grandin (DFO PBS)
#' @export load.iscam.files
load.iscam.files <- function(model.dir,
                             burnin = 1000,
                             thin = 1,
                             verbose = FALSE){
  ## Load all the iscam files for output and input, and return the model object.
  ## If MCMC directory is present, load that and perform calculations for mcmc
  ##  parameters.

  starter.file.name <- "iscam.dat"
  par.file <- "iscam.par"
  rep.file <- "iscam.rep"
  mcmc.file <- "iscam_mcmc.csv"
  mcmc.biomass.file <- "iscam_sbt_mcmc.csv"
  mcmc.recr.file <- "iscam_rt_mcmc.csv"
  mcmc.recr.devs.file <- "iscam_rdev_mcmc.csv"
  mcmc.fishing.mort.file <- "iscam_ft_mcmc.csv"
  mcmc.fishing.mort.u.file <- "iscam_ut_mcmc.csv"
  mcmc.vuln.biomass.file <- "iscam_vbt_mcmc.csv"
  mcmc.proj.file <- "iscammcmc_proj_Gear1.csv"
  mpd.proj.file <- "iscammpd_proj_Gear1.csv"

  model <- list()
  model$path <- model.dir
  ## Get the names of the input files
  inp.files <- fetch.file.names(model.dir, starter.file.name)
  model$dat.file <- inp.files[[1]]
  model$ctl.file <- inp.files[[2]]
  #model$proj.file <- inp.files[[3]]

  ## Load the input files
  model$dat <- read.data.file(model$dat.file)
  model$ctl <- read.control.file(model$ctl.file,
                                 model$dat$num.gears,
                                 model$dat$num.age.gears)
  #model$proj <- read.projection.file(model$proj.file)

  model$par <- read.par.file(file=file.path(model.dir, par.file))
  ## Load MPD results
  model$mpd <- read.report.file(file.path(model.dir, rep.file))
  model.dir.listing <- dir(model.dir)
  ## Set default mcmc members to NA. Later code depends on this.
  model$mcmc <- NA
  ## Set the mcmc path. This doesn't mean it exists.
  model$mcmcpath <- file.path(model.dir, "mcmc")

  ## If it has an 'mcmc' sub-directory, load it
  if(dir.exists(model$mcmcpath)){

   # model$mcmc <- read.mcmc(model$mcmcpath)
    model$mcmc <- read.mcmc(model.dir)
    ## Do the mcmc quantile calculations
    model$mcmccalcs <- calc.mcmc(model,
                                 burnin,
                                 thin,
                                 lower = 0.025,
                                 upper = 0.975)
  }
  model
}

#' Reads iSCAM Data, Control and Projection files
#'
#' @description A function for returning the three types of
#' iSCAM input and output files
#' @param path File path
#' @param filename The filename
#' @author Chris Grandin (DFO PBS)
#' @export fetch.file.names
fetch.file.names <- function(path, ## Full path to the file
                             filename){
  ## Read the starter file and return a list with 3 elements:
  ## 1. Data file name
  ## 2. Control file name
  ## 3. Projection file name

  ## Get the path the file is in
  d <- readLines(file.path(path, filename), warn = FALSE)
  ## Remove comments
  d <- gsub("#.*", "", d)
  ## Remove trailing whitespace
  d <- gsub(" +$", "", d)
  list(file.path(path, d[1]),
       file.path(path, d[2]),
       file.path(path, d[3]))
}


#' Reads iSCAM Rep file
#'
#' @description A function for returning the results of the
#' .rep iscam file
#' @param fn File location
#' @author Chris Grandin (DFO PBS)
#' @export read.report.file
read.report.file <- function(fn){
  # Read in the data from the REP file given as 'fn'.
  # File structure:
  # It is assumed that each text label will be on its own line,
  # followed by one or more lines of data.
  # If the label is followed by a single value or line of data,
  #  a vector will be created to hold the data.
  # If the label is followed by multiple lines of data,
  #  a matrix will be created to hold the data. The matrix might be
  #  ragged so a check is done ahead of time to ensure correct
  #  matrix dimensions.
  #
  # If a label has another label following it but no data,
  #  that label is thrown away and not included in the returned list.
  #
  # A label must start with an alphabetic character followed by
  # any number of alphanumeric characters (includes underscore and .)

  dat <- readLines(fn, warn = FALSE)
  # Remove preceeding and trailing whitespace on all elements,
  #  but not 'between' whitespace.
  dat <- gsub("^[[:blank:]]+", "", dat)
  dat <- gsub("[[:blank:]]+$", "", dat)

  # Find the line indices of the labels
  # Labels start with an alphabetic character followed by
  # zero or more alphanumeric characters
  idx  <- grep("^[[:alpha:]]+[[:alnum:]]*", dat)
  objs <- dat[idx]     # A vector of the object names
  nobj <- length(objs) # Number of objects
  ret  <- list()
  indname <- 0

  for(obj in 1:nobj){
    indname <- match(objs[obj], dat)
    if(obj != nobj){ # If this is the last object
      inddata <- match(objs[obj + 1], dat)
    }else{
      inddata <- length(dat) + 1 # Next row
    }
    # 'inddiff' is the difference between the end of data
    # and the start of data for this object. If it is zero,
    # throw away the label as there is no data associated with it.
    inddiff <- inddata - indname
    tmp <- NA
    if(inddiff > 1){
      if(inddiff == 2){
        # Create and populate a vector
        vecdat <- dat[(indname + 1) : (inddata - 1)]
        vecdat <- strsplit(vecdat,"[[:blank:]]+")[[1]]
        vecnum <- as.numeric(vecdat)
        ret[[objs[obj]]] <- vecnum
      }else if(inddiff > 2){
        # Create and populate a (possible ragged) matrix
        matdat <- dat[(indname + 1) : (inddata - 1)]
        matdat <- strsplit(c(matdat), "[[:blank:]]+")
        # Now we have a vector of strings, each representing a row
        # of the matrix, and not all may be the same length
        rowlengths <- unlist(lapply(matdat, "length"))
        nrow <- max(rowlengths)
        ncol <- length(rowlengths)
        # Create a new list with elements padded out by NAs
        matdat <- lapply(matdat, function(.ele){c(.ele, rep(NA, nrow))[1:nrow]})
        matnum <- do.call(rbind, matdat)
        mode(matnum) <- "numeric"
        ret[[objs[obj]]] <- matnum
      }
    }else{
      # Throw away this label since it has no associated data.
    }
  }
  return(ret)
}



#' Reads iSCAM dat file
#'
#' @description A function for returning the results of the
#' .dat iscam file
#' @param file File location
#' @param verbose should detailed results be printed to console
#' @author Chris Grandin (DFO PBS)
#' @export read.data.file
read.data.file <- function(file = NULL,
                           verbose = FALSE){
  ## Read in the iscam datafile given by 'file'
  ## Parses the file into its constituent parts
  ## And returns a list of the contents

  data <- readLines(file, warn=FALSE)
  tmp <- list()
  ind <- 0

  # Remove any empty lines
  data <- data[data != ""]

  # remove preceeding whitespace if it exists
  data <- gsub("^[[:blank:]]+", "", data)

  # Get the element number for the "Gears" names if present
  dat <- grep("^#.*Gears:.+", data)
  tmp$has.gear.names <- FALSE
  if(length(dat > 0)){
    gear.names.str <- gsub("^#.*Gears:(.+)", "\\1", data[dat])
    gear.names <- strsplit(gear.names.str, ",")[[1]]
    tmp$gear.names <- gsub("^[[:blank:]]+", "", gear.names)
    tmp$has.gear.names <- TRUE
  }


  ## Get the element number for the "CatchUnits" if present
  dat <- grep("^#.*CatchUnits:.+", data)
  if(length(dat > 0)){
    catch.units.str <- gsub("^#.*CatchUnits:(.+)", "\\1", data[dat])
    tmp$catch.units <- gsub("^[[:blank:]]+", "", catch.units.str)
  }

  ## Get the element number for the "IndexUnits" if present
  dat <- grep("^#.*IndexUnits:.+", data)
  if(length(dat > 0)){
    index.units.str <- gsub("^#.*IndexUnits:(.+)", "\\1", data[dat])
    tmp$index.units <- gsub("^[[:blank:]]+", "", index.units.str)
  }

  ## Save the number of specimens per year (comment at end of each age comp
  ##  line), eg. #135 means 135 specimens contributed to the age proportions for
  ##  that year
  age.n <- vector()
  ## Match age comp lines which have N's as comments
  tmp$has.age.comp.n <- FALSE
  pattern <- "^[[:digit:]]{4}[[:space:]]+[[:digit:]][[:space:]]+[[:digit:]][[:space:]]+[[:digit:]][[:space:]]+[[:digit:]].*#([[:digit:]]+).*"
  dat <- data[grep(pattern, data)]
  if(length(dat) > 0){
    for(n in 1:length(dat)){
      age.n[n] <- sub(pattern, "\\1", dat[n])
    }
  }
  ## age.n is now a vector of values of N for the age comp data.
  ## The individual gears have not yet been parsed out, this will
  ##  happen later when the age comps are read in.

  ## Get the element numbers which start with #.
  dat <- grep("^#.*", data)
  ## Remove the lines that start with #.
  dat <- data[-dat]

  ## Remove comments which come at the end of a line
  dat <- gsub("#.*", "", dat)

  ## Remove preceeding and trailing whitespace
  dat <- gsub("^[[:blank:]]+", "", dat)
  dat <- gsub("[[:blank:]]+$", "", dat)

  ## Now we have a nice bunch of string elements which are the inputs for iscam.
  ## Here we parse them into a list structure
  ## This is dependent on the current format of the DAT file and needs to
  ##  be updated whenever the DAT file changes format
  tmp$num.areas  <- as.numeric(dat[ind <- ind + 1])
  tmp$num.groups <- as.numeric(dat[ind <- ind + 1])
  tmp$num.sex    <- as.numeric(dat[ind <- ind + 1])
  tmp$start.yr   <- as.numeric(dat[ind <- ind + 1])
  tmp$end.yr     <- as.numeric(dat[ind <- ind + 1])
  tmp$start.age  <- as.numeric(dat[ind <- ind + 1])
  tmp$end.age    <- as.numeric(dat[ind <- ind + 1])
  tmp$num.gears  <- as.numeric(dat[ind <- ind + 1])

  ## Gear allocation
  tmp$gear.alloc  <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  if(!tmp$has.gear.names){
    tmp$gear.names <- 1:length(tmp$gear.alloc)
  }

  ## Age-schedule and population parameters
  tmp$linf      <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$k         <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$to        <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$lw.alpha  <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$lw.beta   <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$age.at.50.mat <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$sd.at.50.mat  <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$use.mat   <- as.numeric(dat[ind <- ind + 1])
  tmp$mat.vec   <- as.numeric(strsplit(dat[ind <- ind + 1],", ")[[1]])

  ## Delay-difference options
  tmp$dd.k.age   <- as.numeric(dat[ind <- ind + 1])
  tmp$dd.alpha.g <- as.numeric(dat[ind <- ind + 1])
  tmp$dd.rho.g   <- as.numeric(dat[ind <- ind + 1])
  tmp$dd.wk      <- as.numeric(dat[ind <- ind + 1])

  ## Catch data
  tmp$num.catch.obs <- as.numeric(dat[ind <- ind + 1])
  tmp$catch         <- matrix(NA, nrow = tmp$num.catch.obs, ncol = 7)

  for(row in 1:tmp$num.catch.obs){
    tmp$catch[row,] <- as.numeric(strsplit(dat[ind <- ind + 1], "[[:blank:]]+")[[1]])
  }
  colnames(tmp$catch) <- c("year", "gear", "area", "group", "sex", "type", "value")
  ## Abundance indices are a ragged object and are stored as a list of matrices
  tmp$num.indices     <- as.numeric(dat[ind <- ind + 1])
  tmp$num.index.obs   <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$survey.type <- as.numeric(strsplit(dat[ind <- ind + 1], "[[:blank:]]+")[[1]])
  ##nrows <- sum(tmp$nitnobs)
  tmp$indices <- list()
  for(index in 1:tmp$num.indices){
    nrows <- tmp$num.index.obs[index]
    ncols <- 8
    tmp$indices[[index]] <- matrix(NA, nrow = nrows, ncol = ncols)
    for(row in 1:nrows){
      tmp$indices[[index]][row,] <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
    }
    colnames(tmp$indices[[index]]) <- c("iyr","it","gear","area","group","sex","wt","timing")
  }
  ## Age composition data are a ragged object and are stored as a list of matrices
  tmp$num.age.gears <- as.numeric(dat[ind <- ind + 1])
  ##if(!tmp$hasAgeGearNames){
  ##  tmp$ageGearNames <- 1:length(tmp$nagears)
  ##}

  tmp$num.age.gears.vec       <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$num.age.gears.start.age <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$num.age.gears.end.age   <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$eff                     <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$age.comp.flag           <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  tmp$age.comps <- NULL
  ## One list element for each gear (tmp$nagears)
  ## Check to see if there are age comp data
  if(tmp$num.age.gears.vec[1] > 0){
    tmp$age.comps <- list()
    for(gear in 1:tmp$num.age.gears){
      nrows <- tmp$num.age.gears.vec[gear]
      ## 5 of the 6 here is for the header columns
      ncols <- tmp$num.age.gears.end.age[gear] - tmp$num.age.gears.start.age[gear] + 6
      tmp$age.comps[[gear]] <- matrix(NA, nrow = nrows, ncol = ncols)
      for(row in 1:nrows){
        tmp$age.comps[[gear]][row,] <- as.numeric(strsplit(dat[ind <- ind + 1], "[[:blank:]]+")[[1]])
      }
      colnames(tmp$age.comps[[gear]]) <- c("year",
                                           "gear",
                                           "area",
                                           "group",
                                           "sex",
                                           tmp$num.age.gears.start.age[gear]:tmp$num.age.gears.end.age[gear])
    }
  }
  ## Build a list of age comp gear N's
  tmp$age.gears.n <- list()
  start <- 1
  for(ng in 1:length(tmp$num.age.gears.vec)){
    end <- start + tmp$num.age.gears.vec[ng] - 1
    tmp$age.gears.n[[ng]] <- age.n[start:end]
    start <- end + 1
  }
  ## Empirical weight-at-age data
  tmp$num.weight.tab <- as.numeric(dat[ind <- ind + 1])
  tmp$num.weight.obs <- as.numeric(dat[ind <- ind + 1])
  tmp$waa <- NULL

  if(tmp$num.weight.obs > 0){
    ## Parse the weight-at-age data
    nrows       <- tmp$num.weight.obs
    ncols       <- tmp$end.age - tmp$start.age + 6
    tmp$weight.at.age <- matrix(NA, nrow = nrows, ncol = ncols)
    for(row in 1:nrows){
      tmp$weight.at.age[row,] <-
        as.numeric(strsplit(dat[ind <- ind + 1], "[[:blank:]]+")[[1]])
    }
    colnames(tmp$weight.at.age) <- c("year",
                                     "gear",
                                     "area",
                                     "group",
                                     "sex",
                                     tmp$start.age:tmp$end.age)
  }

  ## Annual Mean Weight data
  ## Catch data
  tmp$num.mean.weight <- as.numeric(dat[ind <- ind + 1])
  tmp$num.mean.weight.obs <- as.numeric(dat[ind <- ind + 1])
  if(tmp$num.mean.weight.obs >0){
    tmp$mean.weight.data  <- matrix(NA, nrow = sum(tmp$num.mean.weight.obs), ncol = 7)
    for(row in 1:sum(tmp$num.mean.weight.obs)){
      tmp$mean.weight.data[row,] <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
    }
    colnames(tmp$mean.weight.data) <- c("year",
                                        "meanwt",
                                        "gear",
                                        "area",
                                        "group",
                                        "sex",
                                        "timing")
  }
  tmp$eof <- as.numeric(dat[ind <- ind + 1])
  tmp
}


#' Reads iSCAM control file
#'
#' @description A function for returning the results of the
#' iscam control file
#' @param file File location
#' @param num.gears The number of gears
#' @param num.age.gears The number age-gears
#' @param verbose should detailed results be printed to console
#' @author Chris Grandin (DFO PBS)
#' @export read.control.file
read.control.file <- function(file = NULL,
                              num.gears = NULL,
                              num.age.gears = NULL,
                              verbose = FALSE){
  ## Read in the iscam control file given by 'file'
  ## Parses the file into its constituent parts and returns a list of the
  ##  contents.
  ## num.gears is the total number of gears in the datafile
  ## num.age.gears in the number of gears with age composition information in the
  ##  datafile

  if(is.null(num.gears)){
    cat("You must supply the total number of gears (num.gears). ",
         "Returning NULL.")
    return(NULL)
  }
  if(is.null(num.age.gears)){
    cat("You must supply the number of gears with age composition ",
         "(num.age.gears). Returning NULL.")
    return(NULL)
  }

  data <- readLines(file, warn = FALSE)

  ## Remove any empty lines
  data <- data[data != ""]

  ## Remove preceeding whitespace if it exists
  data <- gsub("^[[:blank:]]+", "", data)

  ## Get the element numbers which start with #.
  dat <- grep("^#.*", data)
  ## Remove the lines that start with #.
  dat <- data[-dat]

  ## Save the parameter names, since they are comments and will be deleted in
  ##  subsequent steps.
  ## To get the npar, remove any comments and preceeding and trailing
  ##  whitespace for it.
  dat1 <- gsub("#.*", "", dat[1])
  dat1 <- gsub("^[[:blank:]]+", "", dat1)
  dat1 <- gsub("[[:blank:]]+$", "", dat1)
  n.par <- as.numeric(dat1)
  param.names <- vector()
  ## Lazy matching with # so that the first instance matches, not any other
  pattern <- "^.*?#([[:alnum:]]+_*[[:alnum:]]*).*"
  for(param.name in 1:n.par){
    ## Each parameter line in dat which starts at index 2,
    ##  retrieve the parameter name for that line
    param.names[param.name] <- sub(pattern, "\\1", dat[param.name + 1])
  }
  ## Now that parameter names are stored, parse the file.
  ##  remove comments which come at the end of a line
  dat <- gsub("#.*", "", dat)

  ## Remove preceeding and trailing whitespace
  dat <- gsub("^[[:blank:]]+", "", dat)
  dat <- gsub("[[:blank:]]+$", "", dat)

  ## Now we have a nice bunch of string elements which are the inputs for iscam.
  ## Here we parse them into a list structure.
  ## This is dependent on the current format of the CTL file and needs to
  ## be updated whenever the CTL file changes format.
  tmp <- list()
  ind <- 0
  tmp$num.params <- as.numeric(dat[ind <- ind + 1])
  tmp$params <- matrix(NA, nrow = tmp$num.params, ncol = 7)
  for(param in 1:tmp$num.params){
    tmp$params[param,] <-
      as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  }
  colnames(tmp$params) <- c("ival","lb","ub","phz","prior","p1","p2")
  ## param.names is retreived at the beginning of this function
  rownames(tmp$params) <- param.names

  ## Age and size composition control parameters and likelihood types
  nrows <- 8
  ncols <- num.age.gears
  tmp$age.size <- matrix(NA, nrow = nrows, ncol = ncols)
  for(row in 1:nrows){
    tmp$age.size[row,] <-
      as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  }
  ## Rownames here are hardwired, so if you add a new row you must add a name
  ##  for it here
  rownames(tmp$age.size) <- c("gearind",
                              "likelihoodtype",
                              "minprop",
                              "comprenorm",
                              "logagetau2phase",
                              "phi1phase",
                              "phi2phase",
                              "degfreephase")
  ## Ignore the int check value
  ind <- ind + 1

  ## Selectivity parameters for all gears
  nrows <- 10
  ncols <- num.gears
  tmp$sel <- matrix(NA, nrow = nrows, ncol = ncols)
  for(row in 1:nrows){
    tmp$sel[row,] <-
      as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  }
  ## Rownames here are hardwired, so if you add a new row you must add a name
  ##  for it here
  rownames(tmp$sel) <- c("iseltype",
                         "agelen50log",
                         "std50log",
                         "nagenodes",
                         "nyearnodes",
                         "estphase",
                         "penwt2nddiff",
                         "penwtdome",
                         "penwttvs",
                         "nselblocks")

  ## Start year for time blocks, one for each gear
  max.block <- max(tmp$sel[10,])
  tmp$start.yr.time.block <- matrix(nrow = num.gears, ncol = max.block)
  for(ng in 1:num.gears){
    ## Pad the vector with NA's to make it the right size if it isn't
    ##  maxblocks size.
    tmp.vec <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
    if(length(tmp.vec) < max.block){
      for(i in (length(tmp.vec) + 1):max.block){
        tmp.vec[i] <- NA
      }
    }
    tmp$start.yr.time.block[ng,] <- tmp.vec
  }

  ## Priors for survey Q, one column for each survey
  tmp$num.indices <- as.numeric(dat[ind <- ind + 1])
  nrows <- 3
  ncols <- tmp$num.indices
  tmp$surv.q <- matrix(NA, nrow = nrows, ncol = ncols)
  for(row in 1:nrows){
    tmp$surv.q[row,] <-
      as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])
  }
  ## Rownames here are hardwired, so if you add a new row you must add a name
  ##  for it here.
  rownames(tmp$surv.q) <- c("priortype",
                            "priormeanlog",
                            "priorsd")

  ## Controls for fitting to mean weight data
  tmp$fit.mean.weight <- as.numeric(dat[ind <- ind + 1])
  tmp$num.mean.weight.cv <- as.numeric(dat[ind <- ind + 1])
  n.vals <- tmp$num.mean.weight.cv
  tmp$weight.sig <-  vector(length = n.vals)
  for(val in 1:n.vals){
    tmp$weight.sig[val] <- as.numeric(dat[ind <- ind + 1])
  }

  ## Miscellaneous controls
  n.rows <- 16
  tmp$misc <- matrix(NA, nrow = n.rows, ncol = 1)
  for(row in 1:n.rows){
    tmp$misc[row, 1] <- as.numeric(dat[ind <- ind + 1])
  }
  ## Rownames here are hardwired, so if you add a new row you must add a name
  ##  for it here.
  rownames(tmp$misc) <- c("verbose",
                          "rectype",
                          "sdobscatchfirstphase",
                          "sdobscatchlastphase",
                          "unfishedfirstyear",
                          "maternaleffects",
                          "meanF",
                          "sdmeanFfirstphase",
                          "sdmeanFlastphase",
                          "mdevphase",
                          "sdmdev",
                          "mnumestnodes",
                          "fracZpriorspawn",
                          "agecompliketype",
                          "IFDdist",
                          "fitToMeanWeight")
  tmp$eof <- as.numeric(dat[ind <- ind + 1])
  tmp
}


#' Reads iSCAM projection file
#'
#' @description A function for returning the results of the
#' iscam projection file
#' @param file File location
#' @param verbose should detailed results be printed to console
#' @author Chris Grandin (DFO PBS)
#' @export read.projection.file
read.projection.file <- function(file = NULL,
                                 verbose = FALSE){
  ## Read in the projection file given by 'file'
  ## Parses the file into its constituent parts
  ##  and returns a list of the contents

  data <- readLines(file, warn = FALSE)

  ## Remove any empty lines
  data <- data[data != ""]

  ## remove preceeding whitespace if it exists
  data <- gsub("^[[:blank:]]+", "", data)

  ## Get the element numbers which start with #.
  dat <- grep("^#.*", data)
  ## remove the lines that start with #.
  dat <- data[-dat]

  ## remove comments which come at the end of a line
  dat <- gsub("#.*", "", dat)

  ## remove preceeding and trailing whitespace
  dat <- gsub("^[[:blank:]]+", "", dat)
  dat <- gsub("[[:blank:]]+$", "", dat)

  ## Now we have a nice bunch of string elements which are the inputs for iscam.
  ## Here we parse them into a list structure.
  ## This is dependent on the current format of the DAT file and needs to
  ##  be updated whenever the proj file changes format.
  tmp <- list()
  ind <- 0

  ## Get the TAC values
  tmp$num.tac  <- as.numeric(dat[ind <- ind + 1])
  for(tac in 1:tmp$num.tac){
    ## Read in the tacs, one, per line
    tmp$tac.vec[tac] <- as.numeric(dat[ind <- ind + 1])
  }

  ## If the tac vector is on one line
  ##tmp$tac.vec <- as.numeric(strsplit(dat[ind <- ind + 1],"[[:blank:]]+")[[1]])

  ## Get the control options vector
  tmp$num.ctl.options <- as.numeric(dat[ind <- ind + 1])
  n.rows <- tmp$num.ctl.options
  n.cols <- 1
  tmp$ctl.options  <- matrix (NA, nrow = n.rows, ncol = n.cols)
  for(row in 1:n.rows){
    tmp$ctl.options[row, 1] <- as.numeric(dat[ind <- ind + 1])
  }
  ## Rownames here are hardwired, so if you add a new row you must add a name
  ##  or it here.
  option.names <- c("syrmeanm",
                    "nyrmeanm",
                    "syrmeanfecwtageproj",
                    "nyrmeanfecwtageproj",
                    "syrmeanrecproj",
                    "nyrmeanrecproj",
                    "shortcntrlpts",
                    "longcntrlpts",
                    "bmin")
  rownames(tmp$ctl.options) <- option.names[1:tmp$num.ctl.options]
  tmp$eof <- as.numeric(dat[ind <- ind + 1])
  tmp
}

#' Reads iSCAM parameter file
#'
#' @description A function for returning the results of the
#' iscam .par file
#' @param file File location
#' @param verbose should detailed results be printed to console
#' @author Chris Grandin (DFO PBS)
#' @export read.par.file
read.par.file <- function(file = NULL,
                          verbose = FALSE){
  ## Read in the parameter estimates file given by 'file'
  ## Parses the file into its constituent parts
  ## And returns a list of the contents

  data <- readLines(file, warn = FALSE)
  tmp <- list()
  ind <- 0

  ## Remove preceeding #
  conv.check <- gsub("^#[[:blank:]]*", "", data[1])
  ## Remove all letters, except 'e'
  ##convCheck <- gsub("[[:alpha:]]+","",convCheck)
  convCheck <- gsub("[abcdfghijklmnopqrstuvwxyz]",
                    "",
                    conv.check,
                    ignore.case = TRUE)
  ## Remove the equals signs
  conv.check <- gsub("=", "", conv.check)
  ## Remove all preceeding and trailing whitespace
  conv.check <- gsub("^[[:blank:]]+", "", conv.check)
  conv.check <- gsub("[[:blank:]]+$", "", conv.check)
  ## Remove the non-numeric parts
  conv.check <- strsplit(conv.check, " +")[[1]]
  conv.check <- conv.check[grep("^[[:digit:]]", conv.check)]
  ## The following values are saved for appending to the tmp list later

  num.params   <- conv.check[1]
  obj.fun.val <-  format(conv.check[2], digits = 6, scientific = FALSE)
  max.gradient <-  format(conv.check[3], digits = 8, scientific = FALSE)

  ##Remove the first line from the par data since we already parsed it and saved the values
  data <- data[-1]

  ## At this point, every odd line is a comment and every even line is the value.
  ## Parse the names from the odd lines (oddData) and parse the
  ## values from the even lines (evenData)
  odd.elem <- seq(1, length(data), 2)
  even.elem <- seq(2, length(data), 2)
  odd.data <- data[odd.elem]
  even.data <- data[even.elem]

  ## Remove preceeding and trailing whitespace if it exists from both
  ##  names and values.
  names <- gsub("^[[:blank:]]+", "", odd.data)
  names <- gsub("[[:blank:]]+$", "", names)
  values <- gsub("^[[:blank:]]+", "", even.data)
  values <- gsub("[[:blank:]]+$", "", values)

  ## Remove the preceeding # and whitespace and the trailing : from the names
  pattern <- "^#[[:blank:]]*(.*)[[:blank:]]*:"
  names <- sub(pattern, "\\1", names)

  ## Remove any square brackets from the names
  names <- gsub("\\[|\\]", "", names)

  data.length <- length(names)
  for(item in 1:(data.length)){
    tmp[[item]] <-
      as.numeric(strsplit(values[ind <- ind + 1], "[[:blank:]]+")[[1]])
  }

  names(tmp) <- names
  tmp$num.params <- num.params
  tmp$obj.fun.val <- as.numeric(obj.fun.val)
  #tmp$max.gradient <- as.numeric(max.gradient)
  tmp
}

#' Reads iSCAM mcmc output files
#'
#' @description A function for returning the results of the
#' iscam mcmc files
#' @param model.dir Folder name
#' @param verbose should detailed results be printed to console
#' @author Chris Grandin (DFO PBS)
#' @export read.mcmc
read.mcmc <- function(model.dir = NULL,
                      verbose = TRUE){
  ## Read in the MCMC results from an iscam model run found in the directory
  ##  model.dir.
  ## Returns a list of the mcmc outputs, or NULL if there was a problem or
  ##  there are no MCMC outputs.


  mcmc.file <- "iscam_mcmc.csv"
  mcmc.biomass.file <- "iscam_sbt_mcmc.csv"
  mcmc.recr.file <- "iscam_rt_mcmc.csv"
  mcmc.recr.devs.file <- "iscam_rdev_mcmc.csv"
  mcmc.fishing.mort.file <- "iscam_ft_mcmc.csv"
  mcmc.fishing.mort.u.file <- "iscam_ut_mcmc.csv"
  mcmc.vuln.biomass.file <- "iscam_vbt_mcmc.csv"
  mcmc.proj.file <- "iscammcmc_proj_Gear1.csv"
  mpd.proj.file <- "iscammpd_proj_Gear1.csv"


  if(is.null(model.dir)){
    cat("You must supply a directory name (model.dir). Returning NULL.")
    return(NULL)
  }
  mcmcfn     <- file.path(model.dir, mcmc.file)
  mcmcsbtfn  <- file.path(model.dir, mcmc.biomass.file)
  mcmcrtfn   <- file.path(model.dir, mcmc.recr.file)
  mcmcrdevfn <- file.path(model.dir, mcmc.recr.devs.file)
  mcmcftfn   <- file.path(model.dir, mcmc.fishing.mort.file)
  mcmcutfn   <- file.path(model.dir, mcmc.fishing.mort.u.file)
  mcmcvbtfn  <- file.path(model.dir, mcmc.vuln.biomass.file)
  mcmcprojfn <- file.path(model.dir, mcmc.proj.file)

  tmp        <- list()
  if(file.exists(mcmcfn)){
    tmp$params <- read.csv(mcmcfn)
  }
  if(file.exists(mcmcsbtfn)){
    sbt        <- read.csv(mcmcsbtfn)
    tmp$sbt    <- extract.group.matrices(sbt, prefix = "sbt")
  }
  if(file.exists(mcmcrtfn)){
    rt         <- read.csv(mcmcrtfn)
    tmp$rt     <- extract.group.matrices(rt, prefix = "rt")
  }
  if(file.exists(mcmcftfn)){
    ft         <- read.csv(mcmcftfn)
    tmp$ft     <- extract.area.sex.matrices(ft, prefix = "ft")
  }
  if(file.exists(mcmcutfn)){
    ut         <- read.csv(mcmcutfn)
    tmp$ut     <- extract.area.sex.matrices(ut, prefix = "ut")
  }
  if(file.exists(mcmcrdevfn)){
    rdev       <- read.csv(mcmcrdevfn)
    tmp$rdev   <- extract.group.matrices(rdev, prefix = "rdev")
  }
  if(file.exists(mcmcvbtfn)){
    vbt        <- read.csv(mcmcvbtfn)
    tmp$vbt    <- extract.area.sex.matrices(vbt, prefix = "vbt")
  }
  tmp$proj <- NULL
  if(file.exists(mcmcprojfn)){
    tmp$proj   <- read.csv(mcmcprojfn)
  }
  tmp
}

extract.group.matrices <- function(data = NULL,
                                   prefix = NULL){
  ## Extract the data frame given (data) by unflattening into a list of matrices
  ## by group. The group number is located in the names of the columns of the
  ## data frame in this format: "prefix[groupnum]_year" where [groupnum] is one
  ## or more digits representing the group number and prefix is the string
  ## given as an argument to the function.
  ## Returns a list of matrices, one element per group.

  if(is.null(data) || is.null(prefix)){
    cat("You must give two arguments (data & prefix). Returning NULL.")
    return(NULL)
  }
  tmp <- list()

  names <- names(data)
  pattern <- paste0(prefix, "([[:digit:]]+)_[[:digit:]]+")
  groups  <- sub(pattern, "\\1", names)
  unique.groups <- unique(as.numeric(groups))
  tmp <- vector("list", length = length(unique.groups))
  ## This code assumes that the groups are numbered sequentially from 1,2,3...N
  for(group in 1:length(unique.groups)){
    ## Get all the column names (group.names) for this group by making a specific
    ##  pattern for it
    group.pattern <- paste0(prefix, group, "_[[:digit:]]+")
    group.names   <- names[grep(group.pattern, names)]
    ## Remove the group number in the name, as it is not needed anymore
    pattern      <- paste0(prefix, "[[:digit:]]+_([[:digit:]]+)")
    group.names   <- sub(pattern, "\\1", group.names)

    # Now, the data must be extracted
    # Get the column numbers that this group are included in
    dat <- data[,grep(group.pattern, names)]
    colnames(dat) <- group.names
    tmp[[group]]  <- dat
  }
  tmp
}

extract.area.sex.matrices <- function(data = NULL,
                                      prefix = NULL){
  ## Extract the data frame given (data) by unflattening into a list of matrices
  ##  by area-sex and gear. The area-sex number is located in the names of the
  ##  columns of the data frame in this format:
  ##  "prefix[areasexnum]_gear[gearnum]_year" where [areasexnum] and [gearnum]
  ##  are one or more digits and prefix is the string given as an argument
  ##  to the function.
  ## Returns a list (area-sex) of lists (gears) of matrices, one element
  ##  per group.

  if(is.null(data) || is.null(prefix)){
    cat("You must give two arguments (data & prefix). Returning NULL.")
    return(NULL)
  }

  names <- names(data)
  pattern <- paste0(prefix, "([[:digit:]]+)_gear[[:digit:]]+_[[:digit:]]+")
  groups  <- sub(pattern, "\\1", names)
  unique.groups <- unique(as.numeric(groups))
  tmp <- vector("list", length = length(unique.groups))
  ## This code assumes that the groups are numbered sequentially from 1,2,3...N
  for(group in 1:length(unique.groups)){
    ## Get all the column names (group.names) for this group by making a
    ##  specific pattern for it
    group.pattern <- paste0(prefix, group, "_gear[[:digit:]]+_[[:digit:]]+")
    group.names <- names[grep(group.pattern, names)]
    ## Remove the group number in the name, as it is not needed anymore
    pattern <- paste0(prefix, "[[:digit:]]+_gear([[:digit:]]+_[[:digit:]]+)")
    group.names <- sub(pattern, "\\1", group.names)
    ## At this point, group.names' elements look like this: 1_1963
    ## The first value is the gear, and the second, the year.
    ## Get the unique gears for this area-sex group
    pattern <- "([[:digit:]]+)_[[:digit:]]+"
    gears <- sub(pattern, "\\1", group.names)
    unique.gears <- unique(as.numeric(gears))
    tmp2 <- vector("list", length = length(unique.gears))
    for(gear in 1:length(unique.gears)){
      gear.pattern <- paste0(prefix, group,"_gear", gear, "_[[:digit:]]+")
      ## Now, the data must be extracted
      ## Get the column numbers that this group are included in
      dat <- data[,grep(gear.pattern, names)]
      ##colnames(dat) <- groupNames
      tmp2[[gear]] <- dat
    }
    tmp[[group]] <- tmp2
  }
  tmp
}

calc.mcmc <- function(model,
                      burnin = 1000,
                      thin = 1,
                      lower = 0.025,
                      upper = 0.975){
  ## Do the mcmc calculations, e.g. quantiles for sbt, recr, recdevs, F, U, vbt
  ## Returns a list of them all
  ##
  ## mcmc - output of the read.mcmc function
  ## burnin - the number of posteriors to remove from the data
  ## thin - the thinning to apply to the posterior samples
  ## lower - lower quantile for confidence interval calcs
  ## upper - upper quantile for confidence interval calcs

  mcmc.thin <- function(mcmc.dat){
    ## apply burnin and thinning to the data

    nm <- names(mcmc.dat)
    mcmc.obj <- mcmc.dat#apply(mcmc.dat, 2, mcmc)
    mcmc.window <- NULL
    for(col in 1:ncol(mcmc.obj)){
      tmp <- window(mcmc.obj[,col],
                    start = burnin + 1,
                    thin = thin)
      mcmc.window <- cbind(mcmc.window, tmp)
    }
    mcmc.window <- as.data.frame(mcmc.window)
    names(mcmc.window) <- nm
    mcmc.window
  }

  probs <- c(lower, 0.5, upper)

  ## Parameters
  mc <- model$mcmc
  params.dat <- mc$params
  params.dat <- strip.areas.groups(params.dat)
  params.dat <- strip.static.params(model, params.dat)
  nm <- names(params.dat)

  p.dat <- params.dat[ , -which(nm %in% c("msy",
                                          "fmsy",
                                          "bmsy",
                                          "umsy",
                                          "ssb",
                                          "bo"))]
  p.names <- names(p.dat)
  p.dat <- mcmc.thin(p.dat)
  p.quants <- apply(p.dat, 2, quantile, prob = probs)

  ## Reference points
  r.dat <- params.dat[ , which(nm %in% c("msy",
                                         "fmsy",
                                         "bmsy",
                                         "umsy",
                                         "bo"))]
  r.names <- names(r.dat)
  r.dat <- mcmc.thin(r.dat)
  r.quants <- apply(r.dat, 2, quantile, prob = probs)

  ## Spawning biomass
  sbt.dat <- mcmc.thin(mc$sbt[[1]])
  sbt.quants <- apply(sbt.dat,
                      2,
                      quantile,
                      prob = probs)
  ## Depletion
  depl.dat <- apply(sbt.dat,
                    2,
                    function(x){x / r.dat$bo})
  depl.quants <- apply(sbt.dat / r.dat$bo,
                       2,
                       quantile,
                       prob = probs)
  ## Recruitment
  recr.dat <- mcmc.thin(mc$rt[[1]])
  recr.mean <- apply(recr.dat,
                     2,
                     mean)
  recr.quants <- apply(recr.dat,
                       2,
                       quantile,
                       prob = probs)
  ## Recruitment deviations
  recr.devs.dat <- mcmc.thin(mc$rdev[[1]])
  recr.devs.quants <- apply(recr.devs.dat,
                            2,
                            quantile,
                            prob = probs)
  ## Vulnerable biomass by gear (list of data frames)
  vuln.dat <- lapply(mc$vbt[[1]], mcmc.thin)
  vuln.quants <- lapply(vuln.dat,
                        function(x){
                          apply(x,
                                2,
                                quantile,
                                prob = lower,
                                na.rm = TRUE)})
  ## Fishing mortalities by gear (list of data frames)
  f.mort.dat <- lapply(mc$ft[[1]], mcmc.thin)
  f.mort.quants <- lapply(f.mort.dat,
                          function(x){
                            apply(x,
                                  2,
                                  quantile,
                                  prob = lower,
                                  na.rm = TRUE)})

  u.mort.dat <- lapply(mc$ut[[1]], mcmc.thin)
  u.mort.quants <- lapply(u.mort.dat,
                          function(x){
                            apply(x,
                                  2,
                                  quantile,
                                  prob = lower,
                                  na.rm = TRUE)})

  sapply(c("p.dat",
           "p.quants",
           "r.dat",
           "r.quants",
           "sbt.dat",
           "sbt.quants",
           "depl.dat",
           "depl.quants",
           "recr.dat",
           "recr.quants",
           "recr.devs.dat",
           "recr.devs.quants",
           "vuln.dat",
           "vuln.quants",
           "f.mort.dat",
           "f.mort.quants",
           "u.mort.dat",
           "u.mort.quants"),
         function(x){get(x)})
}

strip.areas.groups <- function(dat){
  ## This is a hack function to remove the area and group prefixes for the
  ##  mcmc data 'dat'. The reason is that for now we are just working with a
  ##  single group and area, and the extra text in the parameter names are
  ##  confusing, eg. 'ro_gr1' will become just 'ro'. If you make a model with
  ##  more than one group or area this will need to be revisited. Also, this
  ##  removes 'f' which is assumed to be the objective function value. Note
  ##  that q1, q2, q3... will stay the same and m1 and m2 will remain if the
  ##  model was two-sex.

  pnames <- names(dat)
  ## M will only ever be 1 or 2, for each sex
  pnames <- gsub("m_gs1", "m1", pnames)
  pnames <- gsub("m_gs2", "m2", pnames)

  pnames <- gsub("msy1", "msy", pnames)
  pnames <- gsub("fmsy1", "fmsy", pnames)
  pnames <- gsub("SSB1", "ssb", pnames)
  pnames <- gsub("sel_sd([0-9]+)", "selsd\\1", pnames)
  pnames <- gsub("sel_g([0-9]+)", "sel\\1", pnames)
  ## Remove underscores
  names(dat) <- gsub("_+.*", "", pnames)
  ## Remove objective function value
  dat[,names(dat) != "f"]
}

strip.static.params <- function(model, dat){
  ## Strip out the static (non-estimated) parameters from the mcmc output data
  ##  for the given scenario. We only need to see estimated parameters on the
  ##  diagnostic plots. If there are no static parameters, NULL will be returned

  # Check the control file to see which parameters were static
  inp <- as.data.frame(model$ctl$param)
  static <- inp[inp$phz <= 0,]
  snames <- rownames(static)

  ## Now remove those from the mcmc data
  pnames <- names(dat)
  ## remove the log_ stuff from the input parameter names
  snames <- gsub("log_", "", snames)
  ## There will be either one "m" or "m1" and "m2" in pnames.
  ## If "m" is in snames, remove the "m1", and "m2" from pnames as well if they
  ##  exist
  if("m" %in% snames){
    ms <- c("m1", "m2")
    snames <- c(snames, "m1", "m2")
  }
  ## The following also removes "m" in a combined sex model
  dat <- dat[,!(pnames %in% snames)]

  ## Remove static selectivity params
  sel.params <- as.data.frame(model$ctl$sel)
  est.phase <- sel.params["estphase",]
  static.sel <- est.phase < 1
  sel.post.names <- names(dat)[grep("sel[0-9]+",
                                    names(dat))]
  sel.post.names <- sel.post.names[static.sel]
  sel.sd.post.names <- names(dat)[grep("selsd[0-9]+",
                                       names(dat))]
  sel.sd.post.names <- sel.sd.post.names[static.sel]
  dat.names <- names(dat)
  static.sel.inds <- NULL
  if(length(sel.post.names) > 0){
    ## If there are static parameters, remove them
    for(static.sel in 1:length(sel.post.names)){
      static.sel.inds <- c(static.sel.inds,
                           grep(sel.post.names[static.sel],
                                dat.names))
      static.sel.inds <- c(static.sel.inds,
                           grep(sel.sd.post.names[static.sel],
                                dat.names))
    }
    dat <- dat[,-static.sel.inds]
  }
  dat
}


#' Combines indices into a single index using linear modelling
#'
#' @description iSCAM assessments often make use of multiple indices of abundance.
#' The DLMtool data object and MPs currently only make use of a single index.
#' combiSCAMinds is a function that creates a single index from many using
#' linear modelling. It is a simple way of providing initial calculations of
#' management recommendations and it should be noted that this process
#' is important and in a real application would require due diligence (ie
#' peer reviewed data workshop).
#' @param idata List: the indices recorded in a read from an iSCAM data folder, e.g. replist$data$indices
#' @param Year Integer vector: the years of the DLMtool data object ie Data@Year
#' @param fleeteffect Logical: should a fleet effect be added to the linear model?
#' @author T. Carruthers
#' @export iSCAMinds
iSCAMinds<-function(idata,Year,fleeteffect=T){

  ind<-NULL
  for(i in 1:length(idata)){

    edat<-as.data.frame(idata[[i]])
    index<-edat$it/mean(edat$it)
    ind<-rbind(ind,cbind(edat$iyr,rep(i,nrow(edat)),index))

  }

  ind<-as.data.frame(ind)
  names(ind)<-c("Y","FF","I")
  ind$Y<-as.factor(ind$Y)
  ind$FF<-as.factor(ind$FF)

  if(fleeteffect)lm<-lm(log(I)~Y+FF,dat=ind)
  if(!fleeteffect)lm<-lm(log(I)~Y,dat=ind)
  Years<-Year[Year%in%ind$Y]
  newdat<-as.data.frame(cbind(Years,rep(1,length(Years))))
  names(newdat)<-c("Y","FF")
  newdat$Y<-as.factor(newdat$Y)
  newdat$FF<-as.factor(newdat$FF)
  pred<-predict(lm,newdat)
  ind<-rep(NA,length(Year))
  ind[Year%in%Years]<-exp(pred)/(mean(exp(pred)))
  as.data.frame(cbind(Year,ind))

}

#' Combines all iSCAM age composition data across fleets
#'
#' @description iSCAM assessments are often fitted to numerous fleets that have differing
#' age selectivities. iSCAMcomps is a simple way of providing the aggregate catch at age
#' data. It should be noted that this process is important and in a real application would
#' require due diligence (ie peer reviewed data workshop).
#' @param replist S3 class object: the output from a read from an iSCAM data folder
#' @param Year Integer vector: the years of the DLMtool data object ie Data@Year
#' @author T. Carruthers
#' @export iSCAMcomps
iSCAMcomps<-function(replist,Year){

  ny<-length(Year)
  na<-replist$dat$end.age
  CAA<-array(0,c(ny,na))
  compdat<-replist$dat$age.comps
  compN<-replist$dat$age.gears.n

  for(i in 1:length(compdat)){
    comp<-as.data.frame(compdat[[i]])
    cN<-as.numeric(compN[[i]])
    ind<-match(comp$year,Year)
    aind<-match(1:na,names(comp))
    compind<-as.matrix(expand.grid(1:length(ind),aind))
    CAAind<-as.matrix(expand.grid(ind,1:na))
    cNind<-rep(1:length(ind),na)
    CAA[CAAind]<-CAA[CAAind]+ceiling(comp[compind]*cN[cNind])
  }

  CAA
}

LinInt<-function(x){

  nas<-is.na(x)
  ind0<-(1:length(x))[nas]
  ind1<-ind0-1
  ind2<-ind0+1
  x[ind0]<-apply(cbind(x[ind1],x[ind2]),1,mean)
  x
}

#' Reads data from iSCAM file structure into a DLMtool Data object
#'
#' @description A function that uses the file location of a fitted iSCAM
#' model including input files to population the various slots of an
#' data object. iSCAM2OM relies on several functions written by Chris
#' Grandin (DFO PBS).
#' @param iSCAMdir A folder with iSCAM input and output files in it
#' @param Name The name of the operating model
#' @param Source Reference to assessment documentation e.g. a url
#' @param length_timestep How long is a model time step in years
#' (e.g. a quarterly model is 0.25, a monthly model 1/12)
#' @param Author Who did the assessment
#' @author T. Carruthers
#' @importFrom grDevices dev.off gray jpeg png
#' @importFrom coda mcmc
#' @importFrom graphics arrows contour
#' @importFrom stats acf aggregate qnorm window
#' @export iSCAM2Data
iSCAM2Data<-function(iSCAMdir,Name=NULL,Source="No source provided",
                     length_timestep=1,Author="No author provided"){

  message("-- Using function of Chris Grandin (DFO PBS) to extract data from iSCAM file structure --")

  replist<-load.iscam.files(iSCAMdir)

  message("-- End of iSCAM extraction operations --")

  Data <- new("Data")
  if(!is.null(Name)){
    Data@Name<-Name
  }else{
    Data@Name<-replist$path
  }

  catdat<-as.data.frame(replist$dat$catch)
  Data@Cat<-matrix(catdat$value,nrow=1)
  Data@Year<-catdat$year
  ny<-length(Data@Year)
  final3y<-(-2:0)+ny

  inddat<-iSCAMinds(replist$dat$indices,Data@Year,fleeteffect=F) # use linear modelling to get year effect (naive)
  inddat$ind<-LinInt(inddat$ind) # Interpolate NAs
  Data@Ind<-matrix(inddat$ind,nrow=1)
  Data@t<-length(Data@Year)
  Data@AvC<-mean(Data@Cat)
  Data@Dt<-mean(Data@Ind[,final3y])/mean(Data@Ind[,1:3])
  Data@Mort<-replist$mpd$m
  UMSY<-replist$mpd$msy/(replist$mpd$msy+replist$mpd$bmsy)
  FMSY<--log(1-UMSY)
  BMSY<-replist$mpd$bmsy
  MSY<-replist$mpd$msy
  Data@FMSY_M<-FMSY/Data@Mort
  Data@BMSY_B0<-BMSY/replist$mpd$bo
  Data@Cref<-MSY
  Data@Bref<-BMSY
  SSB<-replist$mpd$sbt
  B<-replist$mpd$bt

  SSB0<-replist$mpd$sbo
  depletion<-SSB[length(SSB)-((ny-1):0)]/SSB0
  mult<-mean(inddat$ind/depletion)
  Data@Iref<-Data@BMSY_B0*mult

  Data@vbLinf=replist$dat$linf[1]
  Data@vbK=replist$dat$k[1]
  Data@vbt0=replist$dat$to[1]
  Data@L50= Data@vbLinf*(1-exp(-Data@vbK*(replist$dat$age.at.50.mat-Data@vbt0)))
  A95= -(replist$dat$sd.at.50.mat*log(1/0.95-1)-replist$dat$age.at.50.mat)
  Data@L95=Data@vbLinf*(1-exp(-Data@vbK*(A95-Data@vbt0)))

  Data@MaxAge<-replist$dat$end.age
  FF<-replist$mpd$F
  sel<-FF/apply(FF,1,max)
  selfinal<-apply(sel[final3y,],2,mean)
  AFC<-LinInterp(x=selfinal,y=1:Data@MaxAge,0.05)
  AFS<-LinInterp(x=selfinal,y=1:Data@MaxAge,0.95)
  Data@LFC<-Data@vbLinf*(1-exp(-Data@vbK*(AFC-Data@vbt0)))
  Data@LFS<-Data@vbLinf*(1-exp(-Data@vbK*(AFS-Data@vbt0)))
  Data@CAA<-iSCAMcomps(replist,Data@Year)
  Data@Dep<-mean(depletion[final3y])
  Data@Abun<-mean(B[final3y])
  Data@SpAbun<-mean(SSB[final3y])
  Data@wla=replist$dat$lw.alpha
  Data@wlb=replist$dat$lw.beta
  rec<-replist$mpd$rbar *exp(replist$mpd$delta)*1E6
  SSB<-(replist$mpd$sbt*1000)[1:length(rec)]
  SSBpR<-SSB0/replist$mpd$rbar
  Data@steep<-mean(SRopt(100,SSB,rec,SSBpR,plot=F,type="BH")) # will create a reproducible 1 sample version
  Data@Ref<-Data@Cat[1,ny]
  Data@Ref_type<-"Most recent catches"
  Data@MPrec<-Data@Cat[1,ny]
  Data@MPeff<-1
  Data@LHYear<-ny
  Data@Misc<-list(WARNING="!! this dataset was created automatically using an alpha version of iSCAM2Data and should be treated with caution !!")

  Data

}
tcarruth/MSEtool documentation built on Oct. 19, 2020, 6:09 a.m.