R/predict_mod.R

Defines functions predict_mod

Documented in predict_mod

#' @title Prediction models
#'
#' @description  This function applies Beverton & Holt's yield per recruit model
#'    as well as the Thompson & Bell model. These models predict catch, yield, biomass
#'    and economic values for different
#'    fishing mortality scenarions (in combination with gear changes).
#'
#' @param param a list consisting of following parameters (not all are required):
#' \itemize{
#'   \item \strong{Linf} or \strong{Winf}: infinite length or weight, respectively,
#'      for investigated species in cm [cm],
#'   \item \strong{K}: growth coefficent for investigated species per year [1/year],
#'   \item \strong{t0}: theoretical time zero, at which individuals of this species
#'        hatch,
#'   \item \strong{M}: natural mortality or
#'   \item \strong{Z}: total mortality,
#'   \item \strong{FM}: fishing mortality,
#'   \item \strong{a}: length-weight relationship coefficent (W = a * L^b),
#'   \item \strong{b}: length-weight relationship coefficent (W = a * L^b),
#'   \item \strong{Lr} or \strong{tr}: length or age of recruitment;}
#' additional list objects for the Thompson and Bell model:
#'  \itemize{
#'   \item \strong{midLengths} or \strong{age}: midpoints of the length classes
#'      (length-frequency data) or ages (age composition data),
#'   \item \strong{meanWeight}: vector with mean weight per length group or age class,
#'   \item \strong{meanValue}: vector with mean value per length group or age class,
#' }
#' @param FM_change vector with ascending fishing mortalities (if \code{FM_relative}
#'    is set to TRUE, values can also be relative (only for Thompson and Bell model).
#'    Default are absolute values from 0 to 10). Or
#' @param E_change vector with ascending absolute exploitation rates;
#' @param FM_relative logical; indicating whether \code{FM_change} is relative or
#'    absolute. Default is FALSE (absolute fishing mortalities in \code{FM_change}).
#' @param Lc_change vector with ascending lengths at first capture (Lc), or
#' @param tc_change vector with ascending ages at first capture (tc)
#' @param type indicating which model should be applied: \code{"ypr"} for Beverton
#'    and Holt's yield per recruit model and \code{"ThompBell"} for the Thompson and Bell model
#' @param s_list list with selectivity parameters
#' @param stock_size_1 stock size of smallest size class, if NA values are calculated
#'    relative to a stock size of 1000 individuals
#' @param age_unit in which time unit the data is provided? "month" or "year"
#' @param plus_group if a value is provided, a plus group is created comprising this
#'    size class and all above
#' @param curr.Lc current Lc (length at first capture) if available
#' @param curr.E current exploitation rate if available
#' @param Lmin smallest length group where to start with selection ogive. Not required
#'    for "knife_edge" selection type
#' @param Lincr arbitrary length increment between length groups for estimation of
#'    selection ogive. The smaller the higher the resolution but the slower the model
#'    run. Not required for "knife_edge" selection type
#' @param plot logical; if TRUE results are displayed graphically
#' @param mark logical; if value of choosen points should be displayed in graph (default: TRUE)
#' @param hide.progressbar logical; should progressbar be displayed or hidden? (Default: FALSE)
#'
#' @keywords function prediction ypr
#'
#' @examples
#' #______________________________________
#' # Yiel Per Recruit (YPR) / Beverton and Holt's model
#' #______________________________________
#' # age structured data
#' # Nemipterus marginatus
#' threadfin <- list(Winf = 286, K = 0.37, t0 = -0.2, M = 1.1, tr = 0.4)
#'
#' predict_mod(threadfin, FM_change = seq(0,6,0.1),
#'    tc_change = seq(0.2,1,0.2), type = 'ypr')  #where it is maximal  = MSY
#'
#' # Leiognathus spendens (Pauly, 1980)
#' ponyfish <- list(Winf = 64, K = 1, t0 = -0.2, M = 1.8, tr = 0.2)
#'
#' predict_mod(ponyfish, tc_change = c(0.2,0.3,1.0), type = 'ypr', plot=TRUE)
#'
#' #______________________________________
#' # length structured data
#' # Xiphias gladius (Berkeley and Houde, 1980)
#' swordfish <- list(Linf = 309, K = 0.0949, M = 0.18,
#'                   a = 0.0003, b = 3, Lr = 90)
#'
#' select.list <- list(selecType = 'trawl_ogive', L50 = 120, L75 = 132)
#' #swordfish$midLengths <- seq(60,300,5)
#'
#' output <- predict_mod(param = swordfish, Lc_change = c(100,118,150,180),
#'             s_list = select.list, type = 'ypr', Lmin = 90, Lincr = 8)
#' plot(output)
#'
#' data(hake)
#' hake$Lr <- 35
#' select.list <- list(selecType = 'trawl_ogive', L50 = 50, L75 = 54)
#' output <- predict_mod(param = hake, FM_change = seq(0,3,0.05),
#'                       Lc_change = seq(30,70,1), s_list = select.list,
#'                       type = 'ypr', plot = FALSE, curr.Lc = 50, curr.E = 0.73)
#' plot(output, type = "Isopleth", xaxis1 = "FM", yaxis1 = "Y_R.rel", mark = TRUE)
#'
#' output <- predict_mod(param = hake, E_change = seq(0,1,0.1),
#'                       Lc_change = seq(2,120,2), #s_list = select.list,
#'                       type = 'ypr', plot = FALSE)
#' plot(output, type = "Isopleth", xaxis1 = "E", yaxis1 = "B_R")
#'
#' #______________________________________
#' #      Thompson and Bell model
#' #______________________________________
#' # with age structured data
#' data(shrimps)
#'
#' output <- predict_mod(param = shrimps, FM_change = seq(0.1,20,0.1),
#'      type = "ThompBell", age_unit = "month", plot = TRUE)
#'
#' #______________________________________
#' # with length structured data
#' data(hake)
#' par(mar = c(5, 4, 4, 7))
#' predict_mod(param = hake,FM_change = seq(0.1,3,0.05),
#'      type = 'ThompBell', plot = TRUE)
#'
#' # create list with selectivity information
#' select.list <- list(selecType = 'trawl_ogive', L50 = 50, L75 = 55)
#'
#' output <- predict_mod(param = hake, FM_change = seq(0,2,0.1),
#'      Lc_change = seq(20,70,1),
#'      curr.E = 0.4, curr.Lc = 50,
#'      type = 'ThompBell', s_list = select.list)
#' plot(output, xaxis1 = "FM", yaxis_iso = "Lc", yaxis1 = "B_R", mark = TRUE)
#'
#'
#' @details The Thompson and Bell model incorporates an iteration step simulating the
#'    stock by means
#'    of the \code{\link{stock_sim}} function. In case changes in gear
#'    characteristics -
#'    here measured in terms of Lc or tc, the length or age at first capture,
#'    respectively -
#'    should be explored, a list with selectivity information about the gear has
#'    to be provided and
#'    the prediction models make use of the selectivity \code{\link{select_ogive}}
#'    function.
#'    Sparre and Venema (1998) recommend to treat the last length class always as plus group.
#'    This model is very sensitive to zero observations in the ultimate length
#'    classes. If unrealistic results are returned,
#'    it is recommended to cut length classes with zero observations, group
#'    them in a plus group or to change the interval between
#'    length classes.
#'    Equations which are used in this function assume isometric growth, an
#'    assumption often not met. Further, the assumption that there is no relationship
#'    between the parental stock size and progeny
#'    over a wide range of fishing mortalities or exploitation values,
#'    respectively, is also said to be untrue. By default, the functions assume
#'    knife-edge recruitment and selection of gears (Sparre and Venema, 1998).
#'    If E_change instead of FM_change is used the range is cut at E=0.9, because
#'    higher values of E correspond to unrealistic high values of fishing mortality.
#'    If no selectivity information is given (by use of s_list), knife edge selectivity
#'    with L50 equal to the first argument of Lc_change is assumed.
#'
#'
#' @return A list with the input parameters and dependent on the model type following
#'    list objects:
#' \itemize{
#'   \item \code{type = 'ypr'}
#'   \itemize{
#'      \item \strong{FM}: fishing mortalities,
#'      \item \strong{Lc} or \strong{tc}: lengths or ages at first capture,
#'      \item \strong{list_Lc_runs}: a list with dataframes for each Lc value:
#'      \itemize{
#'        \item \strong{FM_change}: fishing mortalities
#'        \item \strong{E}: expoitation rates
#'        \item \strong{Ty}: mean age in annual yield
#'        \item \strong{LY}: mean length in annual yield
#'        \item \strong{Wy}: mean weight in annual yield
#'        \item \strong{Y_R.rel}: relative yield per recruit (change in catch in
#'            weigth per recruit relative to initial Y/R value)
#'        \item \strong{B_R.rel}: relative biomass per recruit
#'        \item \strong{Y_R}: yield per recurit (catch in weight per recruit)
#'        \item \strong{B_R}: biomass per recruit
#'        \item \strong{B_R.percent}: percentage biomass per recurit in relation to virgin
#'            biomass per recruit
#'      }
#'      \item \strong{df_Es}: a dataframe with references points (columns) for
#'          different Lc values (rows)
#'      \item \strong{df_current}: a dataframe with the exploitation status, yield
#'          and biomass values of current exploitation or selectivity (if E_curr or
#'          Lc_tc_curr provided).
#'   }
#'   \item \code{type = 'ThomBell'}
#'   \itemize{
#'      \item \strong{dt}: delta t,
#'      \item \strong{N}: population number,
#'      \item \strong{dead}: deaths due to natural reasons,
#'      \item \strong{C}: catch,
#'      \item \strong{Y}: yield,
#'      \item \strong{B}: biomass,
#'      \item \strong{V}: value,
#'      \item \strong{totals}: summed up values (total catch, total yield, total value, average biomass),
#'      \item \strong{totC}: total catches for different x factors,
#'      \item \strong{totY}: total yield values for different x factors,
#'      \item \strong{totV}: total values for different x factors,
#'      \item \strong{meanB}: average biomasses for different x factors,
#'      \item \strong{F_change}: fishing mortality changes;
#'   }
#'   \item \code{type = 'ThomBell'} and \code{Lc_change} provided
#'   \itemize{
#'      \item \strong{FM_change}: fishing mortality changes,
#'      \item \strong{Lc_change}: changes in length at first capture,
#'      \item \strong{Lt}: lengths at age,
#'      \item \strong{sel}: probability of capture,
#'      \item \strong{mat_FM_Lc_com.C}: catch matrix for all fishing mortality and Lc/tc combinations,
#'      \item \strong{mat_FM_Lc_com.Y}: yield matrix for all fishing mortality and Lc/tc combinations,
#'      \item \strong{mat_FM_Lc_com.V}: value matrix for all fishing mortality and Lc/tc combinations,
#'      \item \strong{mat_FM_Lc_com.B}: biomass matrix for all fishing mortality and Lc/tc combinations;
#'   }
#' }
#'
#' @importFrom graphics plot
#' @importFrom utils setTxtProgressBar txtProgressBar
#'
#' @references
#' Berkeley, S.A., and Houde, E.D., 1980. Swordfish, Xiphias gladius, dynamics in
#' the Straits of Florida. \emph{ICES C.M.}, 11.
#'
#' Beverton, R.J.H., and Holt, S.J., 1964. Table of yield functions for fishery
#' management. \emph{FAO Fish. Tech. Pap.} 38, 49 p.
#'
#' Beverton, R.J.H., and Holt, S.J., 1966. Manual of methods for fish stock
#' assessment. Pt. 2: Tables of yield functions. \emph{FAO Fisheries Technical Paper},
#' (38)Rev.1:67 p.
#'
#' Boerema, L.K., and J.A. Gulland, 1973. Stock assessment of the Peruvian anchovy
#' (Engraulis ringens) and management of the fishery. \emph{Journal of the Fisheries
#' Board of Canada}, 30(12):2226-2235
#'
#' Garcia, S. and N.P. van Zalinge, 1982. Shrimp fishing in Kuwait: methodology
#' for a joint analysis of the artisanal and industrial fisheries. pp. 119-142 In:
#' Report on the Workshop on assessment of the shrimp stocks of the west coast of
#' the Gulf between Iran and the Arabian Peninsula. Fisheries development in the
#' Gulf. Rome, FAO, FI:DP/RAB/80/015/1, 163 p.
#'
#' Gulland, J.A., 1983. Fish stock assessment: a manual of basic methods.
#' \emph{FAO/Wiley}, New York.
#'
#' Gulland, J.A. and Boerema, L., 1973. Scientific advice on catch levels.
#' \emph{Fish. Bull. (US)} 71:325-335.
#'
#' Jones, R.E. 1957. A much simplified version of the fish yield equation. Doc. No.
#' P. 21. Paper presented at the Lisbon joint meeting of International Commission
#' Northwest Atlantic-Fisheries, International Council for the Exploration of the
#' Sea, and Food and Agriculture Organization of the United Nations. 8 p. [Mimeo].
#'
#' Millar, R.B., and Holst, R., 1997. Estimation of gillnet and hook selectivity using
#' log-linear models. \emph{ICES Journal of Marine Science: Journal du Conseil},
#' 54(3):471-477
#'
#' Pauly, D., 1980. A selection of simple methods for the assessment of tropical fish
#' stocks. \emph{FAO Fisheries Circulars (FAO)}. no. 729.
#'
#' Pauly, D., 1984. Fish population dynamics in tropical waters: a manual for use
#' with programmable calculators. \emph{ICLARM} Stud. Rev. 8, 325 p.
#'
#' Pauly, D. and M. Soriano. 1986. Some practical extensions to Beverton and
#' Holt's relative yield-per-recruit model, p. 491-495. In J.L. Maclean, L.B. Dizon
#' and L.V. Hosillos (eds.) The First Asian Fisheries Forum. Asian Fisheries Society,
#' Manila.
#'
#' Schaefer, M.B., 1954. Some aspects of the dynamics of populations important to the
#' management of the commercial marine fisheries. \emph{Inter-Am. Trop. Tuna Comm.,
#' Bull.} 1(2):27-56.
#'
#' Schaefer, M.B., 1957. A study of the dynamics of the fishery for yellowfin tuna
#' in the eastern tropical Pacific Ocean. \emph{Inter-Am. Trop. Tuna Comm., Bull.}
#' 2:247-268.
#'
#' Sparre, P., and Venema, S.C., 1998. Introduction to tropical fish stock assessment.
#' Part 1. Manual. \emph{FAO Fisheries Technical Paper}, (306.1, Rev. 2). 407 p.
#'
#' @export

predict_mod <- function(param, type, FM_change = NA,
                        E_change = NA,
                        FM_relative = FALSE,
                        Lc_change = NULL,
                        tc_change = NULL,
                        s_list = NA,
                        stock_size_1 = NA, age_unit = 'year', curr.E = NA,
                        curr.Lc = NA,
                        plus_group = NA, Lmin = NA, Lincr = NA, plot = FALSE, mark = TRUE,
                        hide.progressbar = FALSE){
  res <- param
  res$FM_relative = FM_relative

  # Beverton and Holt's ypr
  if(type == "ypr"){
    if(FM_relative) stop(noquote("ypr does not work with relative changes in FM, please provide absolute values."))
    M <- res$M
    K <- res$K
    t0 <- ifelse(!is.null(res$t0),res$t0,0)
    a <- res$a  # might be NULL
    b <- res$b  # might be NULL
    Winf <- res$Winf  # might be NULL
    Linf <- res$Linf # might be NULL
    if("Linf" %in% names(res) & "a" %in% names(res) & "b" %in% names(res)){
      Winf <- a * (Linf ^ b)
    }
    #Linf <- ifelse(!is.null(res$Linf),res$Linf, exp(log(Winf/a)/b))
    # REALLY ? maybe without Linf: then message that Winf has to exist
    #if(is.null(Linf) | is.na(Linf)) stop("Either Linf or Winf with a and b has to be provided!")
    #if(is.null(Winf)) Winf <-  a * Linf ^ b  ###exp((log(Linf) - a)/b) # might still be NULL
    # or              Winf <- exp(log(Linf-a)/b)

    if(length(FM_change) == 1 & is.na(FM_change[1]) &
       length(E_change) == 1 & is.na(E_change[1])){
      FM_change <- seq(0,10,0.1)
      print(noquote("No fishing mortality (FM_change) or exploitation rate (E_change) was provided, a default range for fishing mortality of 0 to 10 is used."))
    }

    # transfer E_change into F_change if provided
    # if(length(FM_change) == 1 & is.na(FM_change[1]) & length(E_change) != 1 & !is.na(E_change[1])){
    #   FM_change <- (E_change * M) / (1 - E_change)
    #   FM_change[FM_change == Inf] <- (0.9999 * M) / (1 - 0.9999)
    # }
    if(length(FM_change) == 1 & is.na(FM_change[1]) &
       length(E_change) != 1 & !is.na(E_change[1])){
      E_change <- E_change[E_change <= 0.9]
      FM_change <- (E_change * M) / (1 - E_change)
    }

    # Recruitment  - knife edge
    tr <- res$tr   # might be NULL
    Lr <- res$Lr   # might be NULL
    if(is.null(tr) & is.null(Lr)) stop("Either the age or the length at recruitment (tr or Lr) has to be provided in param!")
    if(!is.null(Linf)){
      if(is.null(tr)) tr <- VBGF(L=Lr,param = list(Linf=Linf,K=K,t0=t0)) # VBGF(L=Lr,Linf=Linf,K=K,t0=t0)
      if(is.null(Lr)) Lr <- VBGF(t=tr,param = list(Linf=Linf,K=K,t0=t0)) # VBGF(t=tr,Linf=Linf,K=K,t0=t0)
    }


    # Selectivity - knife edge or with selctivtiy ogive
    tc <- res$tc   # might be NULL
    Lc <- res$Lc   # might be NULL
    if(is.null(tc) & is.null(Lc)){
      if("L50" %in% s_list) Lc <- s_list$L50
      if("Lc" %in% s_list) Lc <- s_list$Lc
      #if(!("Lc" %in% s_list) & !("L50" %in% s_list))stop("Either the age or the length at first capture (tc or Lc) has to be provided in param! \n Or provide a Lc value in s_list!")
    }
    if(!is.null(Linf)){
      if(is.null(tc) & !is.null(Lc)) tc <- VBGF(L=Lc, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(L=Lc,Linf=Linf,K=K,t0=t0)
      if(is.null(Lc) & !is.null(tc)) Lc <- VBGF(t=tc, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(t=tc,Linf=Linf,K=K,t0=t0)
      if(is.null(tc_change) & !is.null(Lc_change)) tc_change <- VBGF(L=Lc_change, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(L=Lc_change,Linf=Linf,K=K,t0=t0)
      if(is.null(Lc_change) & !is.null(tc_change)) Lc_change <- VBGF(t=tc_change, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(t=tc_change,Linf=Linf,K=K,t0=t0)
    }
    tc <- c(tc,tc_change)
    Lc <- c(Lc,Lc_change)

    if(length(s_list) > 1){
      selecType <- s_list$selecType
    }else{
      selecType <- "knife_edge"
    }


    # HEART
    list_Lc_runs <- vector("list", length(Lc))
    list_Es <- vector("list", length(Lc))

    # show progress bar only if the loop has more than 1 runs
    if (!hide.progressbar) {
     nlk <- length(Lc)
     if(nlk > 1){
       pb <- txtProgressBar(min=1, max=nlk, style=3)
       counter <- 1
     }
    }

    if(is.null(Lc)) Lc_tc <- tc else Lc_tc <- Lc
    for(i in 1:length(Lc_tc)){

      Lci <- Lc[i]
      tci <- tc[i]

      Z <- (M + FM_change)
      E <- FM_change/Z

      # KNIFE EDGE
      if(length(s_list) == 1 ){#| selecType == "knife_edge"){
        input <- list(Linf=Linf,
                      Winf = Winf,
                      K = K,
                      M = M,
                      t0 = t0,
                      tr = tr,
                      tc = tci)
        output <- ypr(input, FM_change)
        B_R <- output$br
        Y_R <- output$yr
        B_R.rel <- output$rbr
        Y_R.rel <- output$ryr
        deri <- output$derivative
      }

      # SELECTION OGIVE
      if(length(s_list) > 1 ){#& selecType != "knife_edge"){
        if("midLengths" %in% names(res)){
          classes <- as.character(res$midLengths)
          # create column without plus group (sign) if present
          classes.num <- do.call(rbind,strsplit(classes, split="\\+"))
          classes.num <- as.numeric(classes.num[,1])
          Lt <- classes.num
          # Lincr <- res$midLengths[2] - res$midLengths[1]
          # Lmin <- res$midLengths[1] - (Lincr/2)
        }
        if(!"midLengths" %in% names(res)){
          if(is.na(Lmin) | is.na(Lincr)) writeLines("No midpoints of length classes are provided. This can be done using Lmin and Lincr or by a \n midLengths element in param. A standard range from 1cm to Linf * 0.98 by 2cm is assumed")
          Lmin <- ifelse(is.na(Lmin), 1, Lmin)
          Lincr <- ifelse(is.na(Lincr), 2, Lincr)
          # mid length vector
          Lt <- seq(Lmin,(Linf*0.98),Lincr)
        } # make Lt fixed except depending on Linf and as detailed as possible, e.g. Lincr = 0.5 (takes a long time)

        # selectivity
        P <- select_ogive(s_list, Lt =  Lt, Lc = Lci)

        input <- list(Linf = ifelse(length(Linf) > 0, Linf, NA),
                      Winf = ifelse(length(Winf) > 0, Winf, NA),
                      K = K,
                      M = M,
                      t0 = t0,
                      tr = tr,
                      tc = tci)
        output <- ypr_sel(input, FM_change, Lt, P)

        # relative yield and biomass per recruit
        B_R.rel <- output$rbr
        Y_R.rel <- output$ryr
        # derivative
        deri <- output$derivative

        #test
        Y_R <- Y_R.rel * Winf * exp(M * (tr - t0))
        B_R <- B_R.rel * Winf * exp(M * (tr - t0))

        # biased because only prints P for largest Lc value
        #if(i == length(Lc)) plot(Lt, P, type = 'l', ylab = 'Prob of capture',
        #                         main = 'Selectivity function')
      }


      # virgin biomass
      if(0 %in% FM_change){
        Bv_R <- B_R[FM_change == 0]
      }else{
        Bv_R <- B_R[FM_change == min(FM_change,na.rm = TRUE)]
        writeLines(paste0("Biomass was not estimated for a fishing mortality (FM) of 0, thus the virgin biomass corresponds to a FM of ",min(FM_change,na.rm = TRUE)))
      }

      #biomass in percetage of virgin biomass
      B_R.percent <- round((B_R / Bv_R ) * 100, digits = 1)

      #mean age in annual yield
      Ty <- (1 / Z) + tci

      #mean length in the annual yield
      S <- exp(-K * (tci - t0))         # the same:    S <- 1 - (Lci/Linf)
      Ly <- Linf * (1 - ((Z*S)/(Z+K)))

      #mean weight in annual yield
      Wy <- (Z) * Winf *
        ((1/Z) - ((3*S)/(Z+K)) +
           ((3*(S^2))/(Z+(2*K))) - ((S^3)/(Z + (3*K))))


      results.PBH <- data.frame(FM = FM_change,
                                E = E)
      if(length(Ty) > 0) results.PBH$Ty <- Ty
      if(length(Ly) > 0) results.PBH$Ly <- Ly
      if(length(Wy) > 0) results.PBH$Wy <- Wy

      results.PBH$Y_R.rel <- Y_R.rel
      results.PBH$B_R.rel <- B_R.rel

      # WHY NECESSARY???
      if(length(Y_R) > 0) results.PBH$Y_R = Y_R
      if(length(B_R) > 0) results.PBH$B_R = B_R
      if(length(B_R.percent) > 0) results.PBH$B_R.percent = B_R.percent

      list_Lc_runs[[i]] <- results.PBH


      # reference points
      Nmsy <- which.max(Y_R.rel)  #  should be the same as which.min(abs(deri)) which is also labelled Nmax
      deri_pot <- deri[1:Nmsy]
      N01 <- which.min(abs(deri_pot - (deri[1] * 0.1)))
      N05 <- which.min(abs(B_R.percent - 50))  #which.min(abs(deri - (deri[1] * 0.5)))

      df_loop_Es <- data.frame(Lc = ifelse(!is.null(Lci),Lci,NA),
                               tc = ifelse(!is.null(tci),tci,NA),
                               F01 = FM_change[N01],
                               Fmsy = FM_change[Nmsy])
      if(length(B_R.percent) > 0) df_loop_Es$F05 <- FM_change[N05]   # WHY NECESSARY????
      # df_loop_Es$Fmax <- FM_change[Nmax]
      df_loop_Es$E01 <- E[N01]
      df_loop_Es$Emsy <- E[Nmsy]
      if(length(B_R.percent) > 0) df_loop_Es$E05 <- E[N05]    # WHY NECESSARY????
      # df_loop_Es$Emax <- E[Nmax]

      list_Es[[i]] <- df_loop_Es

      # update counter and progress bar
      if (!hide.progressbar) {
       if(nlk > 1){
        setTxtProgressBar(pb, counter)
        counter <- counter + 1
       }
      }
    }

    df_Es <- do.call(rbind,list_Es)

    names(list_Lc_runs) <- paste0("Lc_", Lc_tc)   # names(list_tc_runs) <- tc
    ret <- c(res,list(FM_change = FM_change,
                      Lc = Lc,
                      tc = tc,
                      list_Lc_runs = list_Lc_runs,   #   list_tc_runs = list_tc_runs,
                      df_Es = df_Es))   #   df_Es = df_Es,


    if(!is.na(curr.E) & !is.na(curr.Lc)){
      curr.tc <- VBGF(L=curr.Lc, param = list(Linf=Linf,K=K,t0=t0))
      # current exploitation rate
      curr.F = (M * curr.E)/(1-curr.E)  # curr.F <- (M * curr.E)/(1-curr.E)
      tmpList <- list(Linf=Linf,
                      Winf = Winf,
                      K = K,
                      M = M,
                      t0 = t0,
                      tr = tr,
                      tc = curr.tc)
      if(length(s_list) == 1 | selecType == "knife_edge"){
        tmpRES <- ypr(param = tmpList, FM_change = curr.F)
      }
      if(length(s_list) > 1 & selecType != "knife_edge"){
        P <- select_ogive(s_list, Lt =  Lt, Lc = curr.Lc)
        tmpRES <- ypr_sel(param = tmpList, FM_change = curr.F, Lt, P)
        tmpRES$yr <- tmpRES$ryr * Winf * exp(M * (tr - t0))
        tmpRES$br <- tmpRES$rbr * Winf * exp(M * (tr - t0))
      }

      df_currents <- data.frame(curr.Lc = curr.Lc,
                                curr.tc = curr.tc,
                                curr.E = curr.E,
                                curr.F = curr.F,
                                curr.YPR = tmpRES$yr,        #ypr(curr.F, curr.Lc_tc)       #, type = "length"),           # curr.YPR = ypr(curr.F, curr.Lc_tc, type = "age"),
                                curr.YPR.rel = tmpRES$ryr,     #ypr.rel(curr.F, curr.Lc_tc),   #, type = "length"),   # curr.YPR.rel = ypr.rel(curr.F, curr.Lc_tc, type = "age"),
                                curr.BPR = tmpRES$br,         #bpr(curr.F, curr.Lc_tc),           #, type = "length"),           # curr.BPR = bpr(curr.F, curr.Lc_tc, type = "age"),
                                curr.BPR.rel = tmpRES$rbr)     #bpr.rel(curr.F, curr.Lc_tc))   #, type = "length"))   # curr.BPR.rel = bpr.rel(curr.F, curr.Lc_tc, type = "age"))
      ret$currents <- df_currents
    }
  }

  # Thompson and Bell model
  if(type == "ThompBell"){
    meanWeight <- res$meanWeight
    meanValue <- res$meanValue

    #mortalities
    FM <- res$FM
    if(is.null(res$M) & is.null(res$Z)) stop(noquote("Either M or Z (in 'param') has to be provided!"))
    if(!is.null(res$M)){
      nM <- res$M
      Z <- FM + nM
    }else{
      Z <- res$Z
      nM <- Z - FM
    }

    Linf <- res$Linf
    K <- res$K
    t0 <- ifelse("t0" %in% names(res),res$t0,0)

    # Selectivity - knife edge or with selctivtiy ogive
    tc <- res$tc   # might be NULL
    Lc <- res$Lc   # might be NULL
    if(is.null(tc) & is.null(Lc)){
      if("L50" %in% s_list) Lc <- s_list$L50
      if("Lc" %in% s_list) Lc <- s_list$Lc
      #if(!("Lc" %in% s_list) & !("L50" %in% s_list))stop("Either the age or the length at first capture (tc or Lc) has to be provided in param! \n Or provide a Lc value in s_list!")
    }
    if(!is.null(Linf)){
      if(is.null(tc) & !is.null(Lc)) tc <- VBGF(L=Lc, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(L=Lc,Linf=Linf,K=K,t0=t0)
      if(is.null(Lc) & !is.null(tc)) Lc <- VBGF(t=tc, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(t=tc,Linf=Linf,K=K,t0=t0)
      if(is.null(tc_change) & !is.null(Lc_change)) tc_change <- VBGF(L=Lc_change, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(L=Lc_change,Linf=Linf,K=K,t0=t0)
      if(is.null(Lc_change) & !is.null(tc_change)) Lc_change <- VBGF(t=tc_change, param = list(Linf=Linf,K=K,t0=t0)) # VBGF(t=tc_change,Linf=Linf,K=K,t0=t0)
    }
    tc <- c(tc,tc_change)
    Lc <- c(Lc,Lc_change)

    # age based
    if('age' %in% names(res)) classes <- as.character(res$age)
    # length based
    if('midLengths' %in% names(res)) classes <- as.character(res$midLengths)

    # create column without plus group (sign) if present
    classes.num <- do.call(rbind,strsplit(classes, split="\\+"))
    classes.num <- as.numeric(classes.num[,1])

    if(length(FM_change) == 1 & is.na(FM_change[1]) &
       length(E_change) == 1 & is.na(E_change[1])){
      FM_change <- seq(0,10,0.1)
      print(noquote("No fishing mortality (FM_change) or exploitation rate (E_change) \nwas provided, a default range for the absolute \nfishing mortality of 0 to 10 is used."))
    }

    # transfer E_change into F_change if provided
    if(length(FM_change) == 1 & is.na(FM_change[1]) &
       length(E_change) != 1 & !is.na(E_change[1])){
      E_change <- E_change[E_change <= 0.9]
      FM_change <- (E_change * nM) / (1 - E_change)
    }
    if(length(E_change) == 1 & is.na(E_change[1])){
      E_change <- FM_change / (FM_change + nM)
    }

    Lt <- classes.num  # if age in names(res) Lt here is age because classes.num = age


    # Only FM change provided without Lc_tc change
    if((is.null(tc_change) & is.null(Lc_change))){  #  | length(s_list) == 1){

      #if(is.null(res$FM) | length(res$FM) == 1) stop(noquote("Please provide fishing mortality FM (in 'param') as a vector per size class!"))

      if(is.null(res$FM)) stop(noquote("Please provide fishing mortality FM (in 'param')!"))
      if(length(res$FM) == 1){
        if(length(s_list) > 1 | !is.null(Lc[1])){
          print(noquote("Fishing mortality per length class not povided, using selectivity information to derive fishing mortality per length class."))
          if(length(s_list) == 1){
            s_list <- list(selecType = "knife_edge", L50 = Lc[1])
          }
          sel <- select_ogive(s_list, Lt = Lt)
          FM <- res$FM * sel
        }else{
          stop(noquote("Please provide either fishing mortality FM (in 'param') per length class or a Lc value!"))
        }
      }


      #prediction based on f_change
      if(!FM_relative){
        pred_mat <- as.matrix(FM/max(FM, na.rm = TRUE)) %*% FM_change
      }
      if(FM_relative){
        pred_mat <- as.matrix(FM) %*% FM_change
      }


      pred_res_list <- list()
      for(x7 in 1:length(FM_change)){
        param$Z <- pred_mat[,x7] + nM
        param$FM <- pred_mat[,x7]
        resL <- stock_sim(param, age_unit = age_unit,
                          stock_size_1 = stock_size_1, plus_group = plus_group)
        pred_res_list[[x7]] <- resL$totals
      }

      pred_res_df <- do.call(rbind, pred_res_list)
      pred_res_df$FM_change <- FM_change
      pred_res_df$E_change <- E_change

      res2 <- pred_res_df
      res3 <- c(res,res2)

      # reference points
      Bper <- rep(NA,length(pred_res_df$meanB))
      Bper[1] <- 100
      for(ix in 2:length(Bper)){
        Bper[ix] <- pred_res_df$meanB[ix]/pred_res_df$meanB[1] * 100
      }
      N05 <- which.min(abs(Bper - 50))
      Nmsy <- which.max(pred_res_df$totY)

      if(!is.null(Lc[1]) & !is.null(tc[1])){
        df_Es <- data.frame(Lc = Lc,
                            tc = tc,
                            Fmsy = FM_change[Nmsy],
                            F05 = FM_change[N05],
                            Emsy = E_change[Nmsy],
                            E05 = E_change[N05])
      }else{
        df_Es <- data.frame(Fmsy = FM_change[Nmsy],
                            F05 = FM_change[N05],
                            Emsy = E_change[Nmsy],
                            E05 = E_change[N05])
      }


      ret <- c(res3, list(df_Es = df_Es))


      if(!is.na(curr.E)){
        if(!is.na(curr.Lc)){
          curr.tc <- VBGF(L=curr.Lc, param = list(Linf=Linf, K=K, t0=t0))
        }else curr.tc <- NA
        # current exploitation rate
        curr.F = (nM * curr.E)/(1-curr.E)

        if(is.na(curr.Lc)){
          sel <- (FM / max(FM,na.rm=TRUE))
        }else if(!is.na(curr.Lc)){
          s_list <- list(selecType = "knife_edge", L50 = curr.Lc)
          Lt <- res$midLengths
          sel <- select_ogive(s_list, Lt = Lt, Lc = curr.Lc)
        }
        if(length(s_list) != 1){
          Lt <- res$midLengths
          sel <- select_ogive(s_list, Lt = Lt)
        }

        mati <- sel * curr.F
        param.loop <- res
        param.loop$FM <- mati
        param.loop$Z <- mati + nM
        res2 <- stock_sim(param = param.loop, age_unit = age_unit,
                          stock_size_1 = stock_size_1, plus_group=plus_group)
        mati2 <- res2$totals

        df_currents <- data.frame(curr.Lc = curr.Lc,
                                  curr.tc = curr.tc,
                                  curr.E = curr.E,
                                  curr.F = curr.F,
                                  curr.C = mati2$totC,
                                  curr.Y = mati2$totY,
                                  curr.V = mati2$totV,
                                  curr.B = mati2$meanB)
        ret$currents <- df_currents
      }
    }

    # FM and Lc_tc change provided
    if(!is.null(tc_change) | !is.null(Lc_change)){
      # instead of s_list the outcome of one of the other select functions?

      if(length(s_list) == 1){
        s_list <- list(selecType = "knife_edge", L50 = Lc[1])
      }
      sel <- select_ogive(s_list, Lt = Lt) #classes.num

      sel.list <- list()
      for(x19 in 1:length(Lc)){
        sel.list[[x19]] <- select_ogive(s_list, Lt = Lt, Lc = Lc[x19]) #classes.num
      }
      Lc_mat <- do.call(cbind,sel.list)
      colnames(Lc_mat) <- Lc

      Lc_mat_FM <- Lc_mat   #max(FM, na.rm=TRUE)  # with one it should correspond to actual fishing mortality not to change in mortality (x factor)

      #list with FM_Lc_matrices per FM_change
      FM_Lc_com_mat.list <- list()
      if(!FM_relative){
        for(x20 in 1:length(colnames(Lc_mat_FM))){
          FM_Lc_com_mat.list[[x20]] <- as.matrix(Lc_mat_FM[,x20]) %*% FM_change
          colnames(FM_Lc_com_mat.list[[x20]]) <- FM_change
        }
      }
      if(FM_relative){
        for(x20 in 1:length(colnames(Lc_mat_FM))){
          FM_Lc_com_mat.list[[x20]] <- as.matrix(Lc_mat_FM[,x20] * FM) %*% FM_change
          colnames(FM_Lc_com_mat.list[[x20]]) <- FM_change
        }
      }

      param.loop <- res

      pred.FM_Lc_com_res_loopC_list <- vector("list",length(FM_Lc_com_mat.list))
      pred.FM_Lc_com_res_loopY_list <- vector("list",length(FM_Lc_com_mat.list))
      pred.FM_Lc_com_res_loopB_list <- vector("list",length(FM_Lc_com_mat.list))
      pred.FM_Lc_com_res_loopV_list <- vector("list",length(FM_Lc_com_mat.list))

      if (!hide.progressbar) {
        nlk <- prod(length(FM_Lc_com_mat.list),dim(FM_Lc_com_mat.list[[1]])[2])
        pb <- txtProgressBar(min=1, max=nlk, style=3)
        counter <- 1
      }

      for(x21 in 1:length(FM_Lc_com_mat.list)){  #loop for length of list == Lc changes
        mati <- FM_Lc_com_mat.list[[x21]]

        pred.FM_Lc_com_res_loop1_list <- list()
        for(x22 in 1:dim(mati)[2]){

          param.loop$FM <- mati[,x22]
          param.loop$Z <- mati[,x22] + nM
          res2 <- stock_sim(param = param.loop, age_unit = age_unit,
                            stock_size_1 = stock_size_1, plus_group=plus_group)
          pred.FM_Lc_com_res_loop1_list[[x22]] <- res2$totals

          # update counter and progress bar
          if (!hide.progressbar) {
          setTxtProgressBar(pb, counter)
          counter <- counter + 1
          }
        }
        prev_mat <- do.call(rbind, pred.FM_Lc_com_res_loop1_list)
        prev_matC <- prev_mat[,'totC']
        prev_matY <- prev_mat[,'totY']
        prev_matB <- prev_mat[,'meanB']
        prev_matV <- prev_mat[,'totV']

        pred.FM_Lc_com_res_loopC_list[[x21]] <- prev_matC
        pred.FM_Lc_com_res_loopY_list[[x21]] <- prev_matY
        pred.FM_Lc_com_res_loopB_list[[x21]] <- prev_matB
        pred.FM_Lc_com_res_loopV_list[[x21]] <- prev_matV
      }

      #for catch
      mat_FM_Lc_com.C <- do.call(rbind, pred.FM_Lc_com_res_loopC_list)
      rownames(mat_FM_Lc_com.C) <- Lc
      colnames(mat_FM_Lc_com.C) <- FM_change

      #for yield
      mat_FM_Lc_com.Y <- do.call(rbind, pred.FM_Lc_com_res_loopY_list)
      rownames(mat_FM_Lc_com.Y) <- Lc
      colnames(mat_FM_Lc_com.Y) <- FM_change

      #for biomass
      mat_FM_Lc_com.B <- do.call(rbind, pred.FM_Lc_com_res_loopB_list)
      rownames(mat_FM_Lc_com.B) <- Lc
      colnames(mat_FM_Lc_com.B) <- FM_change

      #for value
      mat_FM_Lc_com.V <- do.call(rbind, pred.FM_Lc_com_res_loopV_list)
      rownames(mat_FM_Lc_com.V) <- Lc
      colnames(mat_FM_Lc_com.V) <- FM_change

      # transvers matrices for plotting (the opposite arrangement from book)
      mat_FM_Lc_com.C <- t(mat_FM_Lc_com.C)
      mat_FM_Lc_com.Y <- t(mat_FM_Lc_com.Y)
      mat_FM_Lc_com.B <- t(mat_FM_Lc_com.B)
      mat_FM_Lc_com.V <- t(mat_FM_Lc_com.V)

      # reference points
      mat_FM_Lc_com.Bper <- matrix(NA,ncol=dim(mat_FM_Lc_com.B)[2],
                                           nrow=dim(mat_FM_Lc_com.B)[1])
      mat_FM_Lc_com.Bper[1,] <- 100
      for(ix in 2:dim(mat_FM_Lc_com.B)[1]){
        mat_FM_Lc_com.Bper[ix,] <- mat_FM_Lc_com.B[ix,]/mat_FM_Lc_com.B[1,] *100
      }
      N05 <- apply(mat_FM_Lc_com.Bper, MARGIN = 2,
                   FUN = function(x) which.min(abs(x - 50)))

      Nmsy <- apply(mat_FM_Lc_com.Y, MARGIN = 2, FUN = which.max)

      if((!is.null(Lc[1]) & !is.null(tc[1])) | (!is.na(Lc[1]) & !is.na(tc[1])) ){
        df_Es <- data.frame(Lc = Lc,
                            tc = tc,
                            Fmsy = FM_change[Nmsy],
                            F05 = FM_change[N05],
                            Emsy = E_change[Nmsy],
                            E05 = E_change[N05])
      }else{
        df_Es <- data.frame(Fmsy = FM_change[Nmsy],
                            F05 = FM_change[N05],
                            Emsy = E_change[Nmsy],
                            E05 = E_change[N05])
      }


      ret <- c(res,
               list(FM_change = FM_change,
                    # FM_relative = FM_relative,
                    E_change = E_change,
                    Lc_change = Lc_change,
                    tc_change = tc_change,
                    Lt = Lt,
                    sel = sel,
                    mat_FM_Lc_com.C = mat_FM_Lc_com.C,
                    mat_FM_Lc_com.Y = mat_FM_Lc_com.Y,
                    mat_FM_Lc_com.V = mat_FM_Lc_com.V,
                    mat_FM_Lc_com.B = mat_FM_Lc_com.B,
                    df_Es = df_Es))


      if(!is.na(curr.E)){
        if(!is.na(curr.Lc)){
          curr.tc <- VBGF(L=curr.Lc, param = list(Linf=Linf, K=K, t0=t0))
        }else curr.tc <- NA

        # current exploitation rate
        curr.F = (nM * curr.E)/(1-curr.E)

        if(is.na(curr.Lc)){
          sel <- FM / max(FM, na.rm = TRUE)
        }else if(!is.na(curr.Lc) | length(s_list) == 1){
          s_list <- list(selecType = "knife_edge", L50 = curr.Lc)
          sel <- select_ogive(s_list, Lt = Lt, Lc = curr.Lc)
        }else if(!is.na(curr.Lc) | length(s_list) != 1){
          sel <- select_ogive(s_list, Lt = Lt, Lc = curr.Lc)
        }
        mati <- sel * curr.F
        param.loop <- res
        param.loop$FM <- mati
        param.loop$Z <- mati + nM
        res2 <- stock_sim(param.loop, age_unit,
                          stock_size_1, plus_group=plus_group)
        mati2 <- res2$totals

        df_currents <- data.frame(curr.Lc = curr.Lc,
                                  curr.tc = curr.tc,
                                  curr.E = curr.E,
                                  curr.F = curr.F,
                                  curr.C = mati2$totC,
                                  curr.Y = mati2$totY,
                                  curr.V = mati2$totV,
                                  curr.B = mati2$meanB)
        ret$currents <- df_currents
      }
    }
  }

  # return results and plot
  class(ret) <- "predict_mod"
  if(plot) plot(ret, mark = mark)
  return(ret)
}

Try the TropFishR package in your browser

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

TropFishR documentation built on June 29, 2018, 5 p.m.