R/model.R

In alko989/s6model: Data-poor stock assessment using a single-sepcies steady-state size-based (s6) model

#' Returns a list with all parameters and functions
#' 
#' Given a Parameters object, it returns as a list all relevant parameters and
#' corresponding functions like pdf and cdf.
#' 
#' 
#' @param p A Parameters object
#' @param FF Numeric. Fishing mortality. This argument is ignored if all optim.* are FALSE 
#' @param calcBRPs Boolean. If true, calculates biological reference points.
#' @param isSurvey boolean, if TRUE a survey selectivity is used for the pdf
#' @param optim.fmsy Logical.
#' @param optim.fmsyr Logical.
#' @param optim.Rrel Logical
#' @return \itemize{
#' \item If optim.* are all FALSE, an invisible list, containing all model parameters, the cdf, the pdf
#' functions. spawning stock biomass (SSB), yield. If calcBRPs is TRUE also Fmsy.
#' \item If optim.fmsy is TRUE, only the yield per Rmax is returned
#' \item If optim.fmsyr is TRUE, only the yield per recruit is returned
#' \item If optim.Rrel is TRUE, only the R/Rmax - 0.5 is returned
#' }
#' @author alko
#' @seealso \code{\link{parameters}}
#' @keywords misc
#' @examples
#' 
#'  getParams()
#' 
#' @export
getParams <- function(p = new("Parameters"),  FF=NULL, calcBRPs=FALSE, isSurvey=FALSE, 
                      optim.fmsy=FALSE, optim.fmsyr = FALSE, optim.Rrel =FALSE) {
  if( ! is(p,"Parameters"))
    stop("Wrong input argument in getParams. Use the Parameters class instead.")  
  Winf <- exp(p@logWinf) * p@scaleWinf
  Fm <- exp(p@logFm) * p@scaleFm
  if(optim.fmsy | optim.fmsyr | optim.Rrel)
    Fm <- FF
  A <- exp(p@logA) * p@scaleA
  n <- exp(p@logn) *p@scalen
  eta_m <- exp(p@logeta_m) * p@scaleeta_m
  eta_S <- exp(p@logeta_S) * p@scaleeta_S
  a <- exp(p@loga) *  p@scalea
  epsilon_a <- exp(p@logepsilon_a) * p@scaleepsilon_a
  epsilon_r <- exp(p@logepsilon_r) * p@scaleepsilon_r
  Wfs <- exp(p@logWfs) * p@scaleWfs
  eta_F <- Wfs/Winf
  if(!(isTRUE(all.equal(p@logeta_F, log(eta_F / p@scaleeta_F))))) warning("Wfs and eta_F do not match! Wfs was used and eta_F was returned correctly")
  p@logeta_F <- log(eta_F / p@scaleeta_F)
  u <-exp(p@logu) * p@scaleu
  M <- p@M
  
  w_r <- w_egg<- 0.001
  Delta <- (log(Winf) - log(w_r)) / (M - 1)
  w <- exp(log(w_r) + (1:M - 1) * Delta)
  
  delta <- diff(w)
  
  psi_F <- (1 + (w / Wfs)^-u )^-1
  
  psi_S <- (1 + (w / (eta_S * Winf))^-u )^-1
  
  psi_m <- (1 + (w / (eta_m * Winf))^-10 )^-1
  m <- a * A * w^(n - 1) +  Fm * psi_F
  m[M] <- Inf
  g <- A*w^n *(1 - (w/Winf)^(1-n) * (epsilon_a + (1 - epsilon_a) * psi_m))
  
  N <- exp(- cumsum((m / g)[-1] * delta)) / g[-1]
  N <- c(1/g[1], N)
  N[M] <- 0
  if(isSurvey)
  {
    fishing <- psi_S * N
  }
  else
  {
    fishing <- psi_F * N
  }
  if( ! (optim.fmsy | optim.fmsyr | optim.Rrel)) {    
    pdfN <-  fishing / sum(fishing * c(delta, 0))
    pdfN.approx <- approxfun(w, pdfN, yleft=0, yright=.Machine$double.xmin)
    cdf <-approxfun(w, cumsum(pdfN * c(delta,0)), yleft=0, yright=1)
  }
  B <- sum((psi_m  * N * w)[-M] * delta)
  Bexpl <- sum((psi_F  * N * w)[-M] * delta)
  Rrel <- 1 - (Winf^(1-n) * w_egg/(epsilon_r * (1 - epsilon_a) * A * B))## * (w_r/w_egg)^(a-1)
  Rp <- epsilon_r * (1 - epsilon_a) * A * Winf ^ (n-1) / w_egg * B
  Y <- Fm * Rrel * sum((psi_F * N * w)[-M] * delta)
  wF <- which.min(abs(Wfs - w)) ## Find the closest weight to Wfs
  
  YR <- Y * 1 / (N[wF] * Rrel * g[wF])# Fm * sum((psi_F * N * w)[-M] * delta)
  
  if(calcBRPs) {
    Fmsy <- optimise(f=getParams, interval=c(0,10), maximum=TRUE, p=p, optim.fmsy=TRUE)$maximum
    FoverFmsy <- Fm/Fmsy
    ##Fmsyr <- optimise(f=getParams, interval=c(0,2), maximum=TRUE, p=p, optim.fmsyr=TRUE)$maximum
    ##FoverFmsyr <- Fm/Fmsyr
    Fcrash <- try(uniroot(f=getParams, interval=c(1e-25,10), p=p, optim.fmsy=TRUE)$root, silent=TRUE)
    Flim <-try(uniroot(f=getParams, interval=c(1e-25,10), p=p,  optim.Rrel = TRUE)$root, silent=TRUE)
  }   
  vb.M <- a * A * Winf^(n-1)*eta_m^(n-1)
  vb.K <- A * Winf^(n-1) / 3
  vb.MK <- vb.M / vb.K
  if(optim.fmsy)
    return(Y)
  if(optim.fmsyr)
    return(YR)
  if(optim.Rrel)
    return(Rrel - 0.5)
  res <- as.list(environment())
  attr(res, "version") <- getVersion()
  return(invisible(res))
}

##' Simulates catch-at-weight data using the s6model
##'
##' Simulate a data set of individual catch weight or 
##' @param samplesize Integer. Number of simulated individuals
##' @param params Object of class \code{Parameters}. Model parameters
##' @param binsize Numeric. Weight class width. If \code{retDF} is FALSE, \code{binsize} is ignored
##' @param keepZeros Logical. If TRUE the resulting data.frame includes weight classes with zero individuals. Otherwise these weight classes are dropped. Ignored if \code{retDF} is FALSE
##' @param retDF Logical. If TRUE a data.frame is returned with columns Weight (weight classes) and Freq (numbers per weight class)
##' @param ... Extra named arguments passed to \code{getParams}
##' @return Invisible list with
##' \itemize{
##'   \item `sample` containing weights of individuals,
##'   \item `parameters` the parameters as a \code{Parameters} object
##'   \item `Fmsy` the Fmsy of the given parameters
##'   \item `df` A data.frame with columns Weight and Freq with number per weight class. This is returned only if \code{retDF} is TRUE
##' }
##' @author alko
##' @export
simulateData3 <- function(samplesize= 1000, params = parameters(), binsize = 5, keepZeros=TRUE, retDF=TRUE, ...)
{
  applyfun <- if(require(parallel)) mclapply else sapply
  sam <- c()
  with(getParams(params, ...), {
    sam <<- simplify2array(applyfun(runif(samplesize), function(u) uniroot(function(x) {cdf(x) - u}, c(w_r, Winf))$root))
  })
  res <- list(sample = sam, parameters = params, Fmsy = getParams(params,calcBRPs=TRUE)$Fmsy)
  if(retDF) res$df <- sample2df(sam, binsize, keepZeros=keepZeros)
  attr(res, "version") <- getVersion()
  return(invisible(res))  
}

##' Convert a vector sample to data.frame with counts per weight class
##'
##' Takes a vector containing individual catch weights and returns a data.frame with numbers per weight class
##' @param sam Numeric vector. Individual weights
##' @param binsize Numeric. Weight class width
##' @param keepZeros Logical. If TRUE the resulting data.frame includes weight classes with zero individuals. Otherwise these weight classes are dropped.
##' @return A data.frame with columns Weight and Freq with number per weight class.
##' @author alko
##' @export
sample2df <- function(sam, binsize, keepZeros=TRUE) {
  df <- as.data.frame(table(cut(sam, seq(0,max(sam, na.rm = TRUE) + binsize, binsize),
                                labels=seq(binsize/2, max(sam, na.rm = TRUE) + binsize/2, binsize) )), stringsAsFactors=FALSE)
  names(df) <- c("Weight","Freq")
  if (! keepZeros) df <- df[df$Freq > 0, ]
  df$Weight <- as.numeric(df$Weight)
  attr(df, "binsize") <- binsize
  df
}

rparam <- function(value, range.sd, range.cv, lb= -Inf, ub = Inf, unif=FALSE)
{
  if(unif)
    return(runif(1,lb, ub))
  if(range.sd == 0) {
    if(range.cv == 0) {
      return(value)
    } else {
      res <- rnorm(1,value, value*range.cv)
      while( (res < lb) || (res > ub) ) {
        res <- rnorm(1,value, value*range.cv)
      }
      return(res)
    }
  } else {
    res <- exp(rnorm(1, log(value), range.sd))
    while((res < lb) || (res > ub) ) {
      res <- rlnorm(1, log(value), range.sd)
    }
    return(res)    
  }
}

##' Random model parameters
##'
##' Random model parameters with constraints in Fmsy and relative recruitment (Rrel = R/Rmax)
##' @param parameter.names Character vector. Names of the parameters
##' @param parameter.value Numeric vector. Parameter mean values
##' @param parameter.sd Numeric vector. Standard deviation of the parameters
##' @param parameter.cv Numeric vector. Coefficient of variation of log transformed parameters
##' @param parameter.lbound Numeric vector. Lower bound of the distributions
##' @param parameter.ubound Numeric vector. Upper bound of the distributions
##' @param parameter.unif Logical vector. Use uniformly distributed parameters
##' @param Rrel.gt Numeric. Relative recruitment constraint. It allows parameters that lead
##' to Rrel at least equal to Rrel.gt
##' @param Fmsy.gt Numeric. Fmsy constraint. It allows parameters that lead
##' to Fmsy at least equal to Fmsy.gt
##' @return \code{Parameters}
##' @author alko
##' @export
getRandomParameters <-
  function(parameter.names=c("A", "n" ,"eta_m","eta_F", "a" ,"Fm","Winf","epsilon_a", "epsilon_r"),
           parameter.value =c(4.5,0.75 , 0.25  ,  0.05 , 0.35,0.25,  1e4 ,    0.8    ,     0.1    ),
           parameter.sd    =c( 0 ,  0  ,   0   ,  0.5  ,  0  , 0  ,   0  ,     0     ,      0     ),
           parameter.cv    =c(0.5,  0  ,  0.3  ,   0   , 0.5 , 0  ,   0  ,    0.1    ,     0.5    ),
           parameter.lbound=c( 0 ,  0  ,  0.01 ,   0   ,0.01 , 0  ,   0  ,     0     ,      0     ),
           parameter.ubound=c(Inf,  0  ,  Inf  ,  Inf  , Inf , 0  ,   0  ,    0.99   ,      Inf   ),
           parameter.unif  =c(F,F,F,F,F,F,F,F,F), Rrel.gt=-Inf, Fmsy.gt=0) {
    while(TRUE) {
      par.vals <- sapply(seq(along.with=parameter.value),
                         function(x) rparam(parameter.value[x], 
                                            parameter.sd[x], 
                                            parameter.cv[x],
                                            parameter.lbound[x],
                                            parameter.ubound[x],
                                            parameter.unif[x]))
      res <- parameters(parameter.names, par.vals, FALSE)
      if( getParams(res)$Rrel >=  Rrel.gt & getParams(res, calcBRPs=T)$Fmsy >= Fmsy.gt) return(res)
    }
  }

##' Random parameters with fixed Winf
##'
##' Shorthand function for a specific asymptotic weight
##' @param winf Numeric. Asymptotic weight
##' @return \code{Parameters}
##' @author alko
##' @rdname getRandomParameters
##' @export
getRandomParameters.fixedWinf <- function(winf, Rrel.gt=-Inf, Fmsy.gt=0) {
  parameter.names <- c("A", "n" ,"eta_m","eta_F", "a" ,"Fm","Winf","epsilon_a", "epsilon_r")
  parameter.value <- c(4.5,0.75 , 0.25  ,  0.05 , 0.35,0.25,  winf ,    0.8    ,     0.1    )
  getRandomParameters(parameter.names, parameter.value, Rrel.gt=Rrel.gt, Fmsy.gt=Fmsy.gt)
}

# ##' Multicore lapply function with progress bar
# ##'
# ##' Wrapper around parallel::mclapply function with text progress bar
# ##' @param X  a vector (atomic or list) or an expressions vector.  Other
# ##' objects (including classed objects) will be coerced by `as.list`.
# ##' @param FUN the function to be applied to (`mclapply`) each element of `X` or (`mcmapply`) in parallel to `...`.
# ##' @param ... Optional arguments to FUN
# ##' @param progressbar Logical. If TRUE a text progress bar is shown.
# ##' @return Result from mclapply
# ##' @note If mclapply is not available, lapply is used instead.
# ##' @author alko
# ##' @export
# tmclapply <- function(X, FUN, ..., progressbar=TRUE){
#   aplfun <- if(require(parallel)) mclapply else lapply
#   start <- Sys.time()
#   if(progressbar)
#     pb <- txtProgressBar(min = 0, max = 100, style=3)
#   results <- local({
#     f <- fifo(tempfile(), open="w+b", blocking=TRUE)
#     if (inherits(parallel:::mcfork(), "masterProcess")) {
#       progress <- 0.0
#       while(progress < 1 && !isIncomplete(f)) {
#         msg <- readBin(f, "double")
#         progress <- progress + as.numeric(msg)
#         if(progressbar)
#           setTxtProgressBar(pb, progress * 100)
#         tt <- (1 - progress)*(difftime(Sys.time(), start, units="mins"))/ progress
#         cat(" ETC:", as.integer(tt), "min(s) and", round((tt - as.integer(tt)) * 60, 0) ,"secs")
#         if( ! progressbar) cat("\r")
#       } 
#       parallel:::mcexit()
#     }
#     res <- aplfun(X, function(x) {
#       rr <- FUN(x)
#       writeBin(1/length(X), f)
#       rr
#     })
#     close(f)
#     if(progressbar) {
#       setTxtProgressBar(pb,100)
#       close(pb)
#     }
#     res
#   })
#   cat(difftime(Sys.time(), start, units="mins"), "mins\n")
#   results
# } 

##' Calculates the Fmsy reference point
##'
##' @param params An object of class \code{Parameters} 
##' @return numeric F that leads to MSY
##' @author alko
##' @export
calcFmsy <- function(params=NULL) {
  if(!require(TMB)) stop("TMB package not installed.")
  if(is.null(params)) return (NULL)
  if(is(params, "Parameters")) {
    params <- as.list(params)
  }
  if( ! is(params,"list"))
    stop("params is of class ", class(params))
  def <- list(n=0.75, epsilon_a=0.8, epsilon_r=0.1, A=4.47, eta_m=0.25, a=0.27, M=1000, u = 10)
  def <- replace(def, names(params), unlist(params))
  obj <- MakeADFun(def, list(logF = log(0.2)), DLL="calcFmsy")
  obj$env$tracemgc <- FALSE
  obj$env$silent <- TRUE; newtonOption(obj=obj,trace=0); config(trace.optimize = 0,DLL="calcFmsy")
  opt <- try(do.call("optim", obj))
  res <- try(sdreport(obj)$val)
  if(is(res, "try-error"))
    return(NULL)
  names(res) <- NULL
  res
}

##' Change the bin size of a weight frequency data.frame
##' 
##' 
##' @param df data.frame with colums 'Weight' and 'Freq'
##' @param binsize numeric, the bin size in grams
##' @param keepZeros logical, if TRUE the bins with zero observation are not removed
##' @export
##' 
##' @examples 
##' ## Simulate a data set with bin size equal to 100 gr
##' dat <- simulateData3(binsize=100)
##' 
##' ## Change the bin size to 200 gr
##' dat <- changeBinsize(dat$df, binsize = 200)
##' @return data.frame with weight frequencies binned with bin size \code{binsize}
changeBinsize <- function(df, binsize = 10, keepZeros = TRUE) {
  sample2df(rep(df$Weight, df$Freq), binsize, keepZeros=keepZeros)
}

##' @export
changeBinsize2 <- function(df, binsize = 10, keepZeros = TRUE, weight.col = "Weight", freq.col = "Freq", verbose = options()$verbose) {
  if(is(df, "list")) { 
    return(lapply (df, changeBinsize2, binsize = binsize, keepZeros = keepZeros, 
                   weight.col = weight.col, freq.col = freq.col, verbose = verbose)) 
  }
  if(dim(df)[1] == 0 ) return( structure(data.frame(Weight = numeric(0), Freq = integer(0)), binsize = binsize))
  cuts <- seq(0, max(df[weight.col], na.rm = TRUE) + binsize, binsize)
  labs <- head(cuts + binsize/2, -1)
  res <- data.frame(Weight = cut(df[[weight.col]], breaks = cuts, labels = labs), Freq = df[[freq.col]])
  res <- rbind(res, data.frame(Weight = labs, Freq = 0))
  res <- aggregate(Freq ~ Weight, data = res, FUN = sum )
  res$Weight <- as.numeric(as.character(res$Weight))
  res$Freq <- as.integer(round(as.numeric(as.character(res$Freq)), digits = 0))
  rl <- rle(res$Freq)
  res <- head(res, nrow(res) - if(tail(rl$values, 1) == 0) tail(rl$lengths, 1) else 0)
  if (! keepZeros) {
    res <- res[df$Freq > 0, ]
  }
  res <- structure(res, binsize = binsize)
  if(verbose) cat("Done")
  res
}

#' @export
changeBinsize3 <- function(df, binsize = 10, keepZeros = TRUE, weight.col = "Weight", freq.col = "Freq", verbose = options()$verbose) {
  if(is(df, "list")) { 
    return(lapply (df, changeBinsize2, binsize = binsize, keepZeros = keepZeros, 
                   weight.col = weight.col, freq.col = freq.col, verbose = verbose)) 
  }
  if(dim(df)[1] == 0 ) return( structure(data.frame(Weight = numeric(0), Freq = integer(0)), binsize = binsize))
  res <- df %>% 
    group_by_(WeightClass = ~cut(weight.col, breaks = seq(0, max(weight.col) + binsize, by = binsize) )) %>% 
    summarise_(Freq = ~sum(freq.col)) %>% 
    mutate_(Weight = ~WeightClass %>% 
             sapply( . %>% as.character %>% getmid)) %>% 
    select(Weight, Freq, WeightClass) %>% 
    `attr<-`("binsize", binsize)
  if (! keepZeros) {
    res <- res[df$Freq > 0, ]
  }
  if(verbose) cat("Done")
  res
}


##' @export
##' @title Find newest file in folder matching pattern
##' @description Given a folder and a pattern, rerurns the file that was modified latest.
##'
##' @param path character, relative or absolute path
##' @param pattern the pattern to search for, used by \code{link{dir}}
##'
##' @seealso \code{\link{dir}}
findLatest <- function(path = ".", pattern = "") {
  files <- dir(path, pattern, ignore.case = TRUE, full.names = TRUE)
  if(length(files) == 0) return(NULL)
  if(length(files) == 1) return(files)
  df <- lapply(files, function(f){
    data.frame(fn = f, mtime = file.info(f)$mtime, stringsAsFactors = FALSE)
  })
  df <- do.call(rbind.data.frame, df)
  latest <- df[order(df$mtime, decreasing = TRUE), ][1,1]  
  latest
}