Bamsim

BAMSim/01cleanup.bat

@@ -0,0 +1,18 @@
+
+del *.r0*
+del *.p0*
+del *.b0*
+del *.log
+del *.rpt
+del *.obj
+del *.htp
+del *.eva
+del *.bar
+del *.tds
+del *.o
+del tmp_admb
+del variance
+del *.dep
+del *.hes
+del *.tmp
+del *.cov

BAMSim/0funcLenProb.R

@@ -0,0 +1,43 @@
+# =================================================================================
+# Compute matrix of proportion of length at age
+# Matrix has nages rows and nlenbins columns, where columns sum to one
+# KWS and MDD  Jan 2025
+#---------------------------------------------------------------------------------
+
+lenprob.fcn <- function(len.mu, len.cv, len.bins, len.binw) {
+  ## INPUT:
+  # len.mu = mean length at age
+  # len.cv = cv of len at age
+  # len.bins = length bins
+  # len.binw = bin width
+
+  ## OUTPUT:
+  #    lenprob = matrix (nages by nlenbins) of proportion at length by age
+
+  #####################################################################################
+  nages <- length(len.mu)
+  nlenbins <- length(len.bins)
+  len_sd <- array(NA, dim = c(nages))
+  cprob_lenvec <- array(NA, dim = c(nlenbins))
+  lenprob <- array(NA, dim = c(nages, nlenbins))
+  for (i in 1:nages) {
+    len_sd[i] <- len.cv * len.mu[i]
+    zscore_lzero <- (0.0 - len.mu[i]) / len_sd[i]
+    cprob_lzero <- pnorm(zscore_lzero)
+    # first length bin
+    zscore_len <- ((len.bins[1] + 0.5 * len.binw) - len.mu[i]) / len_sd[i]
+    cprob_lenvec[1] <- pnorm(zscore_len)
+    lenprob[i, 1] <- cprob_lenvec[1] - cprob_lzero
+    # most other length bins
+    for (l in 2:(nlenbins - 1)) {
+      zscore_len <- ((len.bins[l] + 0.5 * len.binw) - len.mu[i]) / len_sd[i]
+      cprob_lenvec[l] <- pnorm(zscore_len)
+      lenprob[i, l] <- cprob_lenvec[l] - cprob_lenvec[l - 1]
+    }
+    # last length bin is a plus group
+    zscore_len <- ((lenbins[nlenbins] - 0.5 * len.binw) - len.mu[i]) / len_sd[i]
+    lenprob[i, nlenbins] <- 1.0 - pnorm(zscore_len)
+    lenprob[i, ] <- lenprob[i, ] / (1.0 - cprob_lzero) # renormalize to account for any prob mass below 0
+  }
+  return(lenprob)
+}

BAMSim/0funcMSY.R

@@ -0,0 +1,163 @@
+
+# =================================================================================
+# Compute Maximum Sustainable Yield (MSY) reference points
+# Output based on Bev-Holt, Ricker, or Null recruitment models
+# KWS March 2006
+# Last update: Feb 2025
+#---------------------------------------------------------------------------------
+
+MSY.func <- function(steep, R0, M, wgt, sp.frac = 0.0, reprod.age, selL, selD, selZ,
+                     SR.switch = 1, SR.bc = 1, maxF = 2.0, step = 0.001, verbose = TRUE) {
+
+  ## INPUT:
+  ##    steep = steepness parameter
+  ##    R0    = virgin recruitment parameter
+  ##    M     = natural mortality, may be constant or vector
+  ##    wgt   = vector of weight at age, supply vector of 1's if values in numbers desired
+  ##    sp.frac = fraction of year when peak spawning occurs, default is 0.0 (i.e., Jan 1)
+  ##    reprod.age = vector of reproductive output at age
+  ##    selL  = selectivity at age for landings
+  ##    selD  = selectivity at age for dead discards; set to vector of 0's if no discards
+  ##    selZ  = selectivity at age for dead fish
+  ##    SR.switch = switch for SR function (1=BH (default), 2=Ricker, 3=Null)
+  ##    OPTIONAL INPUT
+  ##    SR.bc = lognormal bias correction -- e.g., exp(sigma^2/2)
+  ##    maxF  = maximum F examined, default is 2.0
+  ##    step  = accuracy in MSY calculations (default is 0.001)
+  ##
+  ## OUTPUT:
+  ##    MSY   = maximum sustainable yield in weight, units same as wgt input
+  ##    Lmsy_num = landings in numbers at MSY
+  ##    Fmsy  = fishing rate at MSY
+  ##    Dmsy_wgt  = dead discards at MSY, units same as wgt input
+  ##    Dmsy_num  = dead discards in numbers at MSY
+  ##    spr_msy= spawners per recruit at MSY
+  ##    SPRmsy= spawning potential ratio at MSY (spr_msy/spr_virgin)
+  ##    SSBmsy= spawning output at MSY, in units of fecundity input
+  ##    Rmsy  = equilibrium recruitment at MSY
+  ##    Bmsy  = total biomass (male and female) at MSY
+  ##    Emsy  = exploitation rate at MSY (total kills in number / abundance of fish in number)
+  ##    F = vector of F's used to compute MSY
+  ##    Vectors of landings and discards in numbers and weight and vector of SSB as functions of F
+
+  #####################################################################################
+  # Dependencies: 0funcSR.R
+  #####################################################################################
+
+  #  if (! is.numeric(amin)) stop ("Non-numeric value for minimum age!")
+  #  if (! is.numeric(amax)) stop ("Non-numeric value for maximum age!")
+  if (!is.numeric(steep)) stop("Non-numeric value for steepnes!")
+  if (!is.numeric(R0)) stop("Non-numeric value for R0!")
+  if (!is.numeric(M)) stop("Non-numeric value for M!")
+  if (!is.numeric(wgt)) stop("Non-numeric value for weight at age!")
+  if (!is.numeric(selL)) stop("Non-numeric value for L selectivity at age!")
+  if (!is.numeric(selD)) stop("Non-numeric value for D selectivity at age!")
+  if (!is.numeric(selZ)) stop("Non-numeric value for Z selectivity at age!")
+
+  if (verbose) {
+    if (SR.bc == 1) {
+      cat("*** MSY NOTE: Estimates contain no bias correction.\n")
+    } else {
+      cat("*** NOTE: Estimates contain bias correction (mean unbiased).\n")
+    }
+  }
+
+  nages <- length(selL)
+  ## INITIALIZATION
+  if (length(M) > 1) {
+    M_age <- M
+  } # natural mortality at age (may be constant)
+  else {
+    M_age <- rep(M, nages)
+  }
+
+  F <- seq(0.0, maxF, by = step)
+  spr <- rep(0.0, length(F)) # equilibrium spr at F
+  S_eq <- rep(0.0, length(F)) # equilibrium reproductive output (e.g., SSB) at F
+  R_eq <- rep(0.0, length(F)) # equilibrium recruitment at F
+  B_eq <- rep(0.0, length(F)) # equilibrium biomass at F
+  L_num_eq <- rep(0.0, length(F)) # equilibrium landings at F
+  D_num_eq <- rep(0.0, length(F)) # equilibrium dead discards at F
+  L_wgt_eq <- rep(0.0, length(F)) # equilibrium landings at F
+  D_wgt_eq <- rep(0.0, length(F)) # equilibrium dead discards at F
+  E_eq <- rep(0.0, length(F)) # equilibrium exploitation rate at F (landings only)
+
+  L_age <- rep(0.0, nages) # landings at age
+  D_age <- rep(0.0, nages) # dead discards at age
+  F_age <- rep(0.0, nages) # F at age
+  Z_age <- rep(0.0, nages) # Z at age
+
+  ## Compute virgin spr
+  N0_sp <- rep(1.0, times = nages)
+  N0_sp[1] <- N0_sp[1] * exp(-M_age[1] * sp.frac)
+  for (iage in 2:nages) {
+    N0_sp[iage] <- N0_sp[iage - 1] * exp(-1.0 * ((M_age[iage - 1] * (1.0 - sp.frac)) + (M_age[iage] * sp.frac)))
+  }
+  N0_sp[nages] <- N0_sp[nages] / (1. - exp(-1.0 * M_age[nages]))
+
+  spr_F0 <- sum(N0_sp * reprod.age)
+
+  for (i in 1:length(F)) {
+    FL_age <- F[i] * selL
+    FD_age <- F[i] * selD
+    Z_age <- M_age + F[i] * selZ
+
+    N_age <- N_age_sp <- N_age_msy <- rep(1.0, nages) # N at age
+    for (iage in 2:nages) {
+      N_age[iage] <- N_age[iage - 1] * exp(-1.0 * Z_age[iage - 1])
+    }
+    # last age is pooled
+    N_age[nages] <- N_age[nages - 1] * exp(-1. * Z_age[nages - 1]) /
+      (1. - exp(-1.0 * Z_age[nages]))
+    N_age_sp[1:(nages - 1)] <- N_age[1:(nages - 1)] * exp(-sp.frac * Z_age[1:(nages - 1)])
+    N_age_sp[nages] <- N_age_sp[(nages - 1)] * (exp(-(1 - sp.frac) * Z_age[(nages - 1)]) + sp.frac * Z_age[nages]) / (1 - exp(-Z_age[nages]))
+
+    spr[i] <- sum(N_age_sp * reprod.age)
+
+    R_eq[i] <- SR.eq.fcn(
+      SR.switch = SR.switch, BC = SR.bc, h = steep, R0 = R0,
+      Phi.0 = spr_F0, Phi.F = spr[i]
+    )
+    if (R_eq[i] < 0.0000001) R_eq[i] <- 0.0000001
+
+    N_age_msy <- R_eq[i] * N_age
+    N_age_msy_sp <- R_eq[i] * N_age_sp
+
+    for (iage in 1:nages) {
+      L_age[iage] <- N_age_msy[iage] *
+        (FL_age[iage] / Z_age[iage]) * (1. - exp(-1.0 * Z_age[iage]))
+      D_age[iage] <- N_age_msy[iage] *
+        (FD_age[iage] / Z_age[iage]) * (1. - exp(-1.0 * Z_age[iage]))
+    }
+    S_eq[i] <- sum(N_age_msy_sp * reprod.age)
+    B_eq[i] <- sum(N_age_msy * wgt)
+    L_num_eq[i] <- sum(L_age)
+    D_num_eq[i] <- sum(D_age)
+    L_wgt_eq[i] <- sum(L_age * wgt)
+    D_wgt_eq[i] <- sum(D_age * wgt)
+    E_eq[i] <- (sum(L_age) + sum(D_age)) / sum(N_age_msy)
+  } # END F loop
+
+  msy_out <- max(L_wgt_eq)
+  F_msy_out <- F[L_wgt_eq == msy_out]
+  spr_msy_out <- spr[L_wgt_eq == msy_out]
+  SR_msy_out <- spr_msy_out / spr_F0
+  Lnum_msy_out <- L_num_eq[L_wgt_eq == msy_out]
+  Dnum_msy_out <- D_num_eq[L_wgt_eq == msy_out]
+  Dwgt_msy_out <- D_wgt_eq[L_wgt_eq == msy_out]
+  R_msy_out <- R_eq[L_wgt_eq == msy_out]
+  S_msy_out <- S_eq[L_wgt_eq == msy_out]
+  B_msy_out <- B_eq[L_wgt_eq == msy_out]
+  E_msy_out <- E_eq[L_wgt_eq == msy_out]
+
+  if (F_msy_out == maxF) {
+    cat("*** Fmsy reached a bound.\n")
+  }
+
+  return(list(
+    msy = msy_out, Lmsy_num = Lnum_msy_out, Fmsy = F_msy_out, Dmsy_wgt = Dwgt_msy_out, Dmsy_num = Dnum_msy_out,
+    spr_msy = spr_msy_out, SPRmsy = SR_msy_out, SSBmsy = S_msy_out, Rmsy = R_msy_out,
+    Bmsy = B_msy_out, Emsy = E_msy_out,
+    F = F, L_wgt_eq = L_wgt_eq, L_num_eq = L_num_eq, D_wgt_eq = D_wgt_eq, D_num_eq = D_num_eq, SSB_eq = S_eq
+  ))
+}

BAMSim/0funcMisc.R

@@ -0,0 +1,67 @@
+# =================================================================================
+# Miscellaneous functions.
+# Currently includes functions for geometric mean, power function, baranov catch eqn, and Lorenzen M
+# KWS Feb 2025
+#---------------------------------------------------------------------------------
+
+
+#####################################################################################
+geomean <- function(x) {
+  ## compute geometric mean of a time series x
+  nx <- length(x)
+  y <- prod(x)^(1 / nx)
+  return(y)
+}
+
+#####################################################################################
+pow <- function(x, a, b) {
+  ## INPUT PARAMETERS:
+  ##    x = vector or scalar of independent variable
+  ##    a = coefficient
+  ##    b = exponent
+  ## OUTPUT PARAMETERS:
+  ##    y = values of fcn at x
+
+  y <- a * x^b
+  return(y)
+}
+
+#####################################################################################
+baranov <- function(F.age, Z.age, N.age, wgt.age = NULL) {
+  ## INPUT PARAMETERS:
+  # F.age = vector of fishing mortality rate at age
+  # Z.age = vector of total mortality rate at age including dead discards
+  # N.age = vector of abundance at age
+  # wgt.age = optional vector of weight at age
+  ## OUTPUT PARAMETERS:
+  # A list with total landings in numbers and optionally in weight
+  # This will be "removals" with the addition of dead discards - MDD
+
+  L.age.n <- L.age.w <- rep(NA, times = length(N.age))
+
+  L.age.n <- F.age * N.age * (1.0 - exp(-Z.age)) / Z.age
+  L.n <- sum(L.age.n)
+
+  if (!is.null(wgt.age)) {
+    L.age.w <- wgt.age * L.age.n
+    L.w <- sum(L.age.w)
+  } else {
+    L.w <- NA
+  }
+
+  return(list(L.N = L.n, L.W = L.w))
+} # end fcn baranov
+
+#####################################################################################
+M.lorenzen <- function(length, a = 0, b = -1) {
+  ## FCN based on Lorenzen et al. 2022. Fisheries Research 106327
+  ## INPUT PARAMETERS:
+  # length = vector of length at age in cm
+  # a = intercept parameter
+  # b = slope parameter
+  ## OUTPUT
+  # M at age, which likely needs to be scaled for stock specific values
+  Mage <- exp(a + b * log(length))
+  return(Mage)
+}
+#####################################################################################

BAMSim/0funcSPR.R

@@ -0,0 +1,79 @@
+# =================================================================================
+# Compute Spawning Potential Ratio (SPR) reference points
+# Computes values for F30, F35, F40, F45, F50, and a user-supplied value
+# Computes corresponding values of spawning output if provided equilibrium recruitment value
+# KWS April 2024
+# Last update: Feb 2025
+#---------------------------------------------------------------------------------
+
+SPR.func <- function(nages, max.F = 1.0, R.eq = 0.0, SPR.input = 0.4, sp.frac = 0.0,
+                     reprod.age, M.age, selex.age) {
+
+  ## INPUT:
+  # nages is number of ages, starting with age 1
+  # max.F is the maximum F value over which to compute SPR, default is 1.0
+  # R.eq allows computation of biomass benchmarks, if specified
+  # SPR.input is a user-supplied value if desired, as a decimal value, default is 0.4
+  # sp.frac is the fraction of year when peak spawning occurs, default is 0.0 (i.e., Jan 1)
+  # reprod.age is the reproductive contribution of each age (e.g., female sex ratio X female maturity X female weight)
+  # M.age is age dependent natural mortality
+  # selex.age is selectivity at age representing total mortality from fishing (including landings and discards)
+  #
+  ## OUTPUT:
+  # F.pr = vector of F's used to compute SPR
+  # SPR = vector of SPR as a function of F
+  # s.per.rec = spawners per recruit as a function of F
+  # SSB.eq = reproductive output as a function of F, in units of reprod.age. If R.eq is not specified, this defaults to zero
+  # FX = F providing SPR=X
+  # SSB.FX = SSB corresponding FX. Applies only if R.eq is specified, else defaults to zero.
+
+  F.pr <- SPR <- sp.per.rec <- seq(0, max.F, by = 0.001) # F values over which to compute SPR
+  SSB.eq <- seq(0, length.out = length(F.pr)) # equilibrium SSB associated with each F, conditional on R.eq
+  F.age <- Z.age <- rep(0, nages) # F and Z at age
+  N.pr <- N.pr.sp <- rep(1, nages) # N at age at start of year and at peak spawn time
+
+  Fchoice <- SSB.F30 <- SSB.F35 <- SSB.F40 <- SSB.F45 <- SSB.F50 <- SSB.Fchoice <- NA
+
+  for (i in 1:length(F.pr)) {
+    F.age <- selex.age * F.pr[i]
+    Z.age <- M.age + F.age
+
+    for (iage in 2:nages) {
+      N.pr[iage] <- N.pr[iage - 1] * exp(-Z.age[iage - 1])
+    }
+    N.pr[nages] <- N.pr[nages] / (1 - exp(-Z.age[nages])) # plus group
+    N.pr.sp[1:(nages - 1)] <- N.pr[1:(nages - 1)] * exp(-sp.frac * Z.age[1:(nages - 1)])
+    N.pr.sp[nages] <- N.pr.sp[(nages - 1)] * (exp(-(1 - sp.frac) * Z.age[(nages - 1)]) + sp.frac * Z.age[nages]) / (1 - exp(-Z.age[nages]))
+    sp.per.rec[i] <- sum(N.pr.sp * reprod.age)
+
+    SSB.eq[i] <- sp.per.rec[i] * R.eq
+
+    SPR[i] <- sp.per.rec[i] / sp.per.rec[1]
+  } # end i (F) loop
+
+  if (SPR.input != 0.0) {
+    Fchoice <- F.pr[which.min(abs(SPR - SPR.input))]
+  }
+  F30 <- F.pr[which.min(abs(SPR - 0.3))]
+  F35 <- F.pr[which.min(abs(SPR - 0.35))]
+  F40 <- F.pr[which.min(abs(SPR - 0.40))]
+  F45 <- F.pr[which.min(abs(SPR - 0.45))]
+  F50 <- F.pr[which.min(abs(SPR - 0.50))]
+
+  SSB.F30 <- SSB.F35 <- SSB.F40 <- SSB.F45 <- SSB.F50 <- SSB.Fchoice <- 0.0
+  if (R.eq != 0.0) {
+    SSB.F30 <- SSB.eq[F.pr == F30]
+    SSB.F35 <- SSB.eq[F.pr == F35]
+    SSB.F40 <- SSB.eq[F.pr == F40]
+    SSB.F45 <- SSB.eq[F.pr == F45]
+    SSB.F50 <- SSB.eq[F.pr == F50]
+    if (SPR.input != 0.0) {
+      SSB.Fchoice <- SSB.eq[F.pr == Fchoice]
+    }
+  }
+  return(list(
+    F.pr = F.pr, SPR = SPR, s.per.rec = sp.per.rec, SSB.eq = SSB.eq,
+    F30 = F30, F35 = F35, F40 = F40, F45 = F45, F50 = F50, Fproxy = Fchoice,
+    SSB.F30 = SSB.F30, SSB.F35 = SSB.F35, SSB.F40 = SSB.F40, SSB.F45 = SSB.F45, SSB.F50 = SSB.F50, SSB.Fproxy = SSB.Fchoice
+  ))
+} # end SPR.func