R/VRAP2_functions.R
In MichaelFolkes/VRAP2: A Fork of the NMFS's VRAP Package for RER and populaton viability computations
#' @title Title
#'
#' @param input
#'
#' @return
#' @keywords internal
#'
#' @examples
AEQcalc <- function(input){
# ***** SUB AEQcalc *****
# Compute the AEQs for each age group
# TAB: what are AEQs? Adult Equivalents!
# Sub AEQcalc()
# Dim TmpA As Double 'TAB: added line
# Dim TmpS As Double 'TAB: added line
# Dim Age As Integer 'TAB: added line
#
# TmpA <- 0
# TmpS <- 0
# For Age <- MaxAge% To MinAge% Step -1
# AEQ(Age) <- MatRate(Age) + TmpS * (1 - MatRate(Age)) * TmpA
# TmpA <- AEQ(Age)
# TmpS <- 1 - NatMort(Age)
# Next Age
#
# End Sub
MinAge <- input$MinAge
MaxAge <- input$MaxAge
MatRate <- input$MatRate
NatMort <- input$NatMort
res <- rep(NA, 6)
AEQ <- rep(0,MaxAge);
TmpA <- 0
TmpS <- 0
for(Age in MaxAge:MinAge){
AEQ[Age] <- MatRate[Age] + TmpS * (1 - MatRate[Age]) * TmpA
# res <- rbind(res, round(c(Age, MatRate[Age], TmpS,TmpA, NatMort[Age],AEQ[Age] ),3))
TmpA <- AEQ[Age]
TmpS <- 1 - NatMort[Age]
}
# colnames(res) <- c("Age", "MatRate", "TmpS", "TmpA", "NatMort","AEQ")
# print(res[-1,])
#return(round(AEQ,3) )
return(AEQ)
}#END AEQcalc
#' @title Title
#'
#' @param input
#'
#' @return
#' @keywords internal
#'
#' @examples
AEQcalc_modified <- function(input){
MinAge <- input$MinAge
MaxAge <- input$MaxAge
MatRate <- input$MatRate
NatMort <- input$NatMort
AEQ <- rep(0,MaxAge);
ad <- rep(NA,MaxAge)
sur <- rep(NA,MaxAge)
TmpA <- 0
TmpS <- 0
res <- data.frame(TmpA=numeric(MaxAge), TmpS=numeric(MaxAge), AEQ=numeric(MaxAge), MatRate=numeric(MaxAge))
for(Age in MaxAge:MinAge){
AEQ[Age] <- MatRate[Age] + TmpS * (1 - MatRate[Age]) * TmpA
res$TmpA[Age] <- TmpA
res$TmpS[Age] <- TmpS
res$AEQ[Age] <- AEQ[Age]
res$MatRate[Age] <- MatRate[Age]
TmpA <- AEQ[Age]
TmpS <- 1 - NatMort[Age]
}
return(res)
}#END AEQcalc_modified
#' @title Title
#'
#' @param p
#' @param lastx
#' @param x
#'
#' @return
#' @keywords internal
#'
#' @examples
Autocorrel <- function(p, lastx, x){
# ***** FUNCTION Autocorrel *****
# Compute autocorrelated variable, p <- autocorrelation, lastx <- last value of variable x
# NJS: created 7/9/02
# Function Autocorrel(p As Double, lastx As Double, x As Double) As Double
#
# Autocorrel <- p * lastx + (1 - p) * x
#
# End Function
#
val <- p * lastx + (1 - p) * x;
return(val)
}#END Autocorrel
#' @title Title
#'
#' @param Buffer
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
BufferInit <- function(Buffer, inputs){
# ***** BufferInit *****
# ADJUST TARGET RATE FOR BUFFER EXPLOITATION
#These are the internal function variables
# Defines MxR, SRb, QetR
# Changes DR, an original input
rtn.list <- list()
if(inputs$StepFunc=="POP"){
PBff <- Buffer #population capacity buffer; SRb x this
EBff <- 1 #exploitation buffer; er x this
}
if(inputs$StepFunc=="ER"){
PBff <- 1 #population capacity buffer; SRb x this
EBff <- Buffer #exploitation buffer; er x this
}
BufSRb <- PBff * inputs$BSRb #adjust capacity upward
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
MxR <- BufSRb * inputs$AveEnv #never used except printed out; Max Recruits adj by average environment
QetR <- inputs$BSRa * inputs$AveEnv * inputs$DL2
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
MxR <- BufSRb * inputs$AveEnv #never used except printed out;
QetR <- inputs$AveEnv / ((1 / BufSRb) + (1 / inputs$BSRa) * (1 / inputs$DL2))
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
MxR <- inputs$BSRa * inputs$AveEnv * BufSRb / 2.71828 #never used except printed out;
QetR <- inputs$AveEnv * inputs$BSRa * inputs$DL2 * exp(-inputs$DL2 / BufSRb)
}
#cat(paste("Max Recruits", MxR))
rtn.list$DR <- inputs$DR
#DR is one of the original inputs; as is DL2; QetR is local variable
if(inputs$DR == 1){
if(QetR == 0){
rtn.list$DR <- 0
}else{
rtn.list$DR <- inputs$DL2 / QetR #TAB: don't like this, changing DR which is one of the original inputs
}
}
#NumBreakPoints is an original input
#This is just the ER to use for a particular sim (Buffer); only changed if StepFunc="ER"
rtn.list$BufTargetU <- c()
for(Break in 1:(inputs$NumBreakPoints + 1)){
if(inputs$NumBreakPoints > 1){ #we have 3 or more Buffer targets
if(Break == 1){
rtn.list$BufTargetU[1] <- min(inputs$TargetU[2] * EBff, inputs$TargetU[1]) #use the smaller of these 2
}else{
rtn.list$BufTargetU[Break] <- inputs$TargetU[Break] * EBff
}
}else{ # we have 1 or 2 Buffer targets
rtn.list$BufTargetU[Break] <- inputs$TargetU[Break] * EBff
}
}
rtn.list$BufSRb <- BufSRb #adjusted if StepFunc <- Pop
return(rtn.list)
}#END BufferInit
#' @title Title
#'
#' @param TempCohort
#' @param Cohort
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
CompAgeCohort <- function(TempCohort, Cohort, inputs){
# ***** CompAgeCohort *****
# Sub CompAgeCohort(TempCohort() As Double)
# Dim Age As Integer 'counter for ages
#
# For Age <- MaxAge% To MinAge% + 1 Step -1
# Cohort(Age) <- TempCohort(Age - 1)
# Next Age
#
# recruitment I think
# Cohort(MinAge%) <- Cohort(MinAge% - 1)
#
# End Sub
#TempCohort is Cohort - mortality and fish that return to spawn
#Cohort age t in next year is cohort age t-1 this year
# for(Age in inputs$MaxAge:(inputs$MinAge + 1)){
# Cohort[Age] <- TempCohort[Age - 1]
# }
Cohort[inputs$MaxAge:(inputs$MinAge + 1)] <- TempCohort[inputs$MaxAge:(inputs$MinAge + 1)-1]
#This 'ages' recruitment in Cohort[1]
# EEH: this only works if MinAge=2, because Cohort[1] was set in CompRecruits (not Cohort[minage-1])
# EEH: why isn't TempCohort[inputs$MinAge - 1] <- Cohort[inputs$MinAge - 1]?
Cohort[inputs$MinAge] <- Cohort[inputs$MinAge - 1]
return(Cohort)
}#END CompAgeCohort
#' @title Title
#'
#' @param Alpha
#' @param Beta
#'
#' @return
#' @keywords internal
#'
#' @examples
CompBetaVariate <- function(Alpha, Beta){
# ***** CompBetaVariate *****
# This subroutine generates a beta r.v.
# It does so by using the gamma distribution
# Mean <- alpha/(alpha+beta)
# Variance <- (a*b)/((a+b)^2*(a+b+1))
# TAB: not sure this is correct way of getting beta distr
# since Beta(a,b) <- (Gamma(a)*Gamma(b))/(Gamma(a+b))
# ***********************************************
# Function CompBetaVariate(Alpha As Double, Beta As Double)
# Dim G1 As Double
# Dim G2 As Double
#
# If CONSTRUN <- True Then
# CompBetaVariate <- Alpha / (Alpha + Beta) 'expected value
# Else
# If Alpha <= 1 Then 'TAB: changed ALP to Alpha in this line
# G1 <- Gam1(Alpha, 1)
# Else
# G1 <- Gam2(Alpha, 1)
# End If
#
# If Beta <= 1 Then 'TAB: changed BET to Beta in this line
# G2 <- Gam1(Beta, 1)
# Else
# G2 <- Gam2(Beta, 1)
# End If
#
# CompBetaVariate <- G1 / (G1 + G2)
# End If
# End Function
#EEH for testing purposes. Replace with rbeta(alpha,beta) later
#local variables
# Dim G1 As Double
# Dim G2 As Double
# if(get("CONSTRUN", pkg_globals)){ #no random variables
# val <- Alpha / (Alpha + Beta) # expected value
# }else{
G1 <- rgamma(1, Alpha, scale=1)
G2 <- rgamma(1, Beta, scale=1)
val <- G1 / (G1 + G2)
# }
return(val)
}#END CompBetaVariate
#' @title Title
#'
#' @param Regime
#' @param Year
#' @param inputs
#' @param BufTargetU
#' @param Cohort
#' @param AEQ
#' @param YearStats
#'
#' @return
#' @keywords internal
#'
#' @examples
CompEscpmnt <- function(Regime, Year, inputs, BufTargetU, Cohort, AEQ, YearStats){
#If you have breakpoints, then Regime refers to which of these to use
# INITIALIZE PARAMETERS, VARIABLES, AND ARRAYS
Converge <- "No"
MaxAge=inputs$MaxAge
# SET EXPLOITATION RATE TARGET AFTER ADJUSTMENT FOR MANAGEMENT ERROR.
# SET INITIAL FISHERY SCALE To 1.0 UNLESS TARGET U=0
ErrorScale <- 1
if(inputs$MgmtError == "YES") {
#ErrorScale <- rgamma(1, inputs$GammaMgmtA, scale=inputs$GammaMgmtB)
#keep sampling until realized ER is <1
#this is a break from the original method of converting all er>=1 to .99
#as R lacks a do loop, this does the same.
#(need to sample error before evaluating it, and 'while' can't do that)
#idea taken from:
#https://stackoverflow.com/questions/4357827/do-while-loop-in-r
repeat{
ErrorScale <- eval(parse(text=inputs$MgmtError.FUN))
if(BufTargetU[Regime] * ErrorScale< 1 ) {break} }
}#END if(inputs$MgmtError == "YES") {
if(BufTargetU[Regime] > 0){
BufTargetUError <- BufTargetU[Regime] * ErrorScale
if(BufTargetUError >= 1){
#cat("Warning - Target HR Reduced to 0.99\n")
BufTargetUError <- 0.99 #not allowed to take all fish?
}
FishScale <- 1 #there is fishing
}else{
TempCohort <- Cohort*(1-inputs$MatRate)
Escpmnt <- Cohort * inputs$MatRate
Escpmnt[Escpmnt < 1] <- 0
AEQMort <- rep(0, inputs$MaxAge)
TotAEQMort <- 0
TotEscpmnt <- sum(Escpmnt)
Converge <- "Yes" #don't need to go through the while loop
}#END if(BufTargetU[Regime] > 0){
n <- 0
while(Converge == "No"){
n <- n + 1
# COMPUTE PRETERMINAL MORTALITY AND UPDATE COHORT
tmp=FishScale * inputs$PTU
tmp[tmp>1] <- 1 #can't take more fish than are there
PTMort <- tmp * Cohort
TempCohort <- Cohort - PTMort
# for(Age in inputs$MinAge:inputs$MaxAge){
# if(FishScale * inputs$PTU[Age] < 1){ #FishScale is 1 or 0
# PTMort[Age] <- FishScale * Cohort[Age] * inputs$PTU[Age] #take a fraction of fish
# }else{
# PTMort[Age] <- Cohort[Age] #take all fish at that age
# }
# TempCohort[Age] <- Cohort[Age] - PTMort[Age] #pre-term mort
# } #end for loop age
# COMPUTE MATURE RUN AND UPDATE COHORT
MatRun <- TempCohort * inputs$MatRate
TempCohort <- TempCohort - MatRun
# for(Age in inputs$MinAge:inputs$MaxAge){
# MatRun[Age] <- TempCohort[Age] * inputs$MatRate[Age]
# TempCohort[Age] <- TempCohort[Age] - MatRun[Age]
# } #end for loop Age
# COMPUTE MATURE MORTALITY AND ESCAPEMENT
# MatU, MatMort, Escpmnt, MatRun are all 1:MaxAge vectors
tmp <- FishScale * inputs$MatU
tmp[tmp>1] <- 1 #can't take more fish than are there
MatMort <- tmp * MatRun
Escpmnt <- MatRun - MatMort #spawners
Escpmnt[Escpmnt < 1] <- 0
# for(Age in inputs$MinAge:inputs$MaxAge){
# if(FishScale * inputs$MatU[Age] < 1){
# MatMort[Age] <- FishScale * MatRun[Age] * inputs$MatU[Age] #these are the fish that return and are taken
# }else{ #all mature fish taken
# MatMort[Age] <- MatRun[Age]
# }
# Escpmnt[Age] <- MatRun[Age] - MatMort[Age] #spawners
# if(Escpmnt[Age] < 1) Escpmnt[Age] <- 0
# }# end Age for loop
# COMPUTE AEQ TOTAL MORTALITY EXPLOITATION RATE
AEQMort <- AEQ * PTMort + MatMort
TotAEQMort <- sum(AEQMort)
TotEscpmnt <- sum(Escpmnt)
# for(Age in inputs$MinAge:inputs$MaxAge){
# AEQMort[Age] <- (AEQ[Age] * PTMort[Age]) + MatMort[Age]
# TotAEQMort <- TotAEQMort + AEQMort[Age]
# TotEscpmnt <- TotEscpmnt + Escpmnt[Age]
# #!!! MinAge fish (age 2) not counted
# if(Age > inputs$MinAge) TotAdultEscpmnt <- TotAdultEscpmnt + Escpmnt[Age]
# } #end for loop
# COMPARE ACTUAL AND TARGET; SCALE EXPLOITATION RATES IF NECESSARY
# if any of these met, exit
if(n > 100 | BufTargetU[Regime] <= 0 | (TotAEQMort + TotEscpmnt) < 1){
Converge <- "Yes"
}else{ #none are met, so determine if converged
ActualU <- TotAEQMort / (TotAEQMort + TotEscpmnt)
PercentError <- abs((ActualU - BufTargetUError) / BufTargetUError)
if(PercentError > inputs$ConvergeCrit){
FishScale <- FishScale * BufTargetUError / ActualU
}else{
Converge <- "Yes"
}
}
} #end of while
#not needed inside while loop
#TotAdultEscpmnt <- sum(Escpmnt[(inputs$MinAge+1):inputs$MaxAge])
YearStats$AEQMort[Year,] <- AEQMort
YearStats$Escpmnt[Year,] <- Escpmnt
YearStats$TotAdultEscpmnt[Year] <- sum(Escpmnt[(inputs$MinAge+1):inputs$MaxAge])
YearStats$TotAEQMort[Year] <- TotAEQMort
YearStats$TotEscpmnt[Year] <- TotEscpmnt
YearStats$TempCohort <- TempCohort
YearStats$tmpdat[Year] <- ErrorScale
#mf addition of other output to YearStats:
YearStats$TargetUError[Year] <- BufTargetU[Regime] * ErrorScale
YearStats$Cohort[Year,] <- Cohort
if(exists("FishScale")){
YearStats$ptu.scaled[Year,] <- FishScale * inputs$PTU
YearStats$MatU.scaled[Year,] <- FishScale * inputs$MatU
} else{
YearStats$ptu.scaled[Year,] <- NA
YearStats$MatU.scaled[Year,] <- NA
}
return(YearStats)
}#END CompEscpmnt
#' @title Title
#'
#' @param inputs
#' @param CohortBeforeNatMort
#'
#' @return
#' @keywords internal
#'
#' @examples
CompNatMort <- function(inputs, CohortBeforeNatMort){
# ***** CompNatMort *****
# Let the number of fish in each age class decrease according to
# the natural mortality in that age class.
# Sub CompNatMort()
# Dim Age% 'counter for ages
#
# For Age% <- MinAge% - 1 To MaxAge%
# Cohort(Age%) <- Cohort(Age%) * (1 - NatMort(Age%))
# Next Age%
#
# End Sub
#updates global array Cohort
# CohortAfterNatMort <- CohortBeforeNatMort
# for(Age in (inputs$MinAge - 1):inputs$MaxAge){
# CohortAfterNatMort[Age] <- CohortBeforeNatMort[Age] * (1 - inputs$NatMort[Age])
# }
CohortAfterNatMort <- CohortBeforeNatMort * (1 - inputs$NatMort)
return(CohortAfterNatMort)
}#END CompNatMort
#' @title Title
#'
#' @param YearStats
#' @param Year
#' @param inputs
#' @param repvars
#' @param staticvars
#' @param BufSRb
#'
#' @return
#' @keywords internal
#'
#' @examples
CompRecruits <- function(YearStats, Year, inputs, repvars, staticvars, BufSRb){
#if StepFunc="Pop" then SRb (capacity) is being changed.
if(inputs$StepFunc=="POP"){
SRb=BufSRb
}else{
SRb=inputs$BSRb
}
#cummulative stats
LastRanFlow <- repvars$LastRanFlow
LastRanError <- repvars$LastRanError
LastRanMarine <- repvars$LastRanMarine
#Horrid; renaming TotEscpmnt to Escpmnt
if(inputs$EscChoice == "YES"){
Escpmnt <- YearStats$TotAdultEscpmnt[Year]
}else{ Escpmnt <- YearStats$TotEscpmnt[Year] }
if(Escpmnt < 1){ Escpmnt <- 0 }
# COMPUTE AEQ RECRUITMENT
if(toupper(inputs$SRType) == "HOC2"){
if(inputs$BSRa * Escpmnt > SRb){
AEQRecruits <- SRb
}else{
AEQRecruits <- inputs$BSRa * Escpmnt
}
}
if(toupper(inputs$SRType) == "RIC2"){
AEQRecruits <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb)
}
if(toupper(inputs$SRType) == "BEV2"){
AEQRecruits <- 0
if(Escpmnt > 0) AEQRecruits <- 1/((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
if(any(inputs$SRType==c("HOC2","RIC2","BEV2")) & inputs$SurvScale == "YES"){
#StockScalar <- GammaSample(inputs$SRErrorA, inputs$SRErrorB)
#StockScalar <- rgamma(1, inputs$SRErrorA, scale=inputs$SRErrorB)
StockScalar <- eval(parse(text=inputs$SRError.FUN))
AEQRecruits <- StockScalar * AEQRecruits
}
if(any(inputs$SRType==c("HOC3", "RIC3", "BEV3"))){
if(inputs$TrndCycF == "Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
#RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- eval(parse(text=inputs$FlowError.FUN))
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
} #Autoc
if(inputs$TrndCycF == "Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
#EEH: In old code this was GammaFlowA which meant the global was overwritten
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF == "Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0){ RanFlow <- 0 }
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC3")
if(inputs$BSRa * Escpmnt > SRb){ r <- SRb * FWS
}else{ r <- inputs$BSRa * Escpmnt * FWS }
if(inputs$SRType == "RIC3"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * FWS
}
if(inputs$SRType == "BEV3"){
r <- 0;
if(Escpmnt > 0) r <- FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * repvars$LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
}#end 3 parameter cases
if(any(inputs$SRType==c("HOC4", "RIC4", "BEV4"))){
#Set Marine Surv Params
if(inputs$TrndCycM == "Autoc"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
#RanMarine <- GammaSample(inputs$GammaMarA, inputs$GammaMarB)
RanMarine <- rgamma(1, inputs$GammaMarA, scale=inputs$GammaMarB)
RanMarine <- Autocorrel(inputs$TCM1, LastRanMarine, RanMarine)
}
}
if(inputs$TrndCycM == "Trend"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
Mean <- Trend(inputs$TCM1, Year, inputs$MarAve, inputs$TCM2)
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- Mean
}else{
GammaMarTrA <- Mean * Mean / Var
GammaMarTrB <- Var / Mean
#RanMarine <- GammaSample(GammaMarTrA, GammaMarTrB)
RanMarine <- rgamma(1, GammaMarTrA, scale=GammaMarTrB)
}
}
}
if(inputs$TrndCycM == "Cycle"){
Mean <- Cycle(inputs$TCM1, inputs$TCM2, inputs$TCM3, Year)
# should use TCF4 as scalar in place of MarAve
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- inputs$MarAve
}else{
GammaMarCyA <- Mean * Mean / Var
GammaMarCyB <- Var / Mean
#RanMarine <- GammaSample(GammaMarCyA, GammaMarCyB)
RanMarine <- rgamma(1, GammaMarCyA, scale=GammaMarCyB)
}
}
if(RanMarine < 0) RanMarine <- 0
LastRanMarine <- RanMarine
MS <- RanMarine^inputs$BSRc
if(MS < 0) MS <- 0
#Set the flow parameters
if(inputs$TrndCycF =="Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
}
if(inputs$TrndCycF =="Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF =="Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- inputs$FlowAve
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0) RanFlow <- 0
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC4"){
if(inputs$BSRa * Escpmnt > SRb){
r=SRb * MS * FWS
}else{ r <- inputs$BSRa * Escpmnt * MS * FWS }
}
if(inputs$SRType == "RIC4"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * MS * FWS
}
if(inputs$SRType == "BEV4"){
r <- 0;
if(Escpmnt > 0) r <- MS * FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
} #end param 4 cases
if(inputs$depen == "YES" & Escpmnt < inputs$DL1 + 1){
if(Escpmnt < inputs$DL2 + 1){
fac <- Escpmnt / inputs$DL2 * inputs$DR
}else{
fac <- ((1 - (inputs$DL1 - Escpmnt) / (inputs$DL1 - inputs$DL2)) * (1 - inputs$DR)) + inputs$DR }
if(fac < 0) fac <- 0
AEQRecruits <- fac * AEQRecruits
}
# COMPUTE RECRUITS that are age 1;
# RecruitsAtAge1 is from Recruits() and is "total fraction of age 1 ind that eventually return"
# AEQRecruits/RecruitsAtAge1 <- Age 1 or Cohort[1]
CohortAge1 <- AEQRecruits / staticvars$RecruitsAtAge1
# ADD MARINE SURVIVAL IF STOCK RECRUIT FUNCTION IS SPAWNER TO SMOLT
if(inputs$MarSurv == "YES"){
BetaVariate <- CompBetaVariate(inputs$BetaMarA, inputs$BetaMarB)
CohortAge1 <- CohortAge1 * BetaVariate
}
repvars$Cohort[1] <- CohortAge1
repvars$LastRanFlow <- LastRanFlow
repvars$LastRanError <- LastRanError
repvars$LastRanMarine <- LastRanMarine
#mf addition:
repvars$AEQRecruits <- AEQRecruits
return(repvars)
}#END CompRecruits
#' @title Title
#'
#' @param YearStats
#' @param Year
#' @param inputs
#' @param repvars
#' @param staticvars
#' @param BufSRb
#'
#' @return
#' @keywords internal
#'
#' @examples
CompRecruits2 <- function(YearStats, Year, inputs, repvars, staticvars, BufSRb){
# SRb=inputs$BSRb
#cummulative stats
LastRanFlow <- repvars$LastRanFlow
LastRanError <- repvars$LastRanError
LastRanMarine <- repvars$LastRanMarine
#Horrid; renaming TotEscpmnt to Escpmnt
if(inputs$EscChoice == "YES"){
Escpmnt.sum <- YearStats$TotAdultEscpmnt[Year]
}else{ Escpmnt.sum <- YearStats$TotEscpmnt[Year] }
if(Escpmnt.sum < 1){ Escpmnt.sum <- 0 }
#COMPUTE AEQ RECRUITMENT
AEQRecruits <- eval(parse(text=inputs$SR.FUN))
if(inputs$SRType == "HOC2"){
if(inputs$BSRa * Escpmnt > SRb){
AEQRecruits <- SRb
}else{
AEQRecruits <- inputs$BSRa * Escpmnt
}
}
if(inputs$SRType == "RIC2"){
AEQRecruits <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb)
}
if(inputs$SRType == "BEV2"){
AEQRecruits <- 0
if(Escpmnt > 0) AEQRecruits <- 1/((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
if(any(inputs$SRType==c("HOC2","RIC2","BEV2")) & inputs$SurvScale == "YES"){
#StockScalar <- GammaSample(inputs$SRErrorA, inputs$SRErrorB)
#StockScalar <- rgamma(1, inputs$SRErrorA, scale=inputs$SRErrorB)
StockScalar <- eval(parse(text = inputs$SRError.FUN))
AEQRecruits <- StockScalar * AEQRecruits
}
if(any(inputs$SRType==c("HOC3", "RIC3", "BEV3"))){
if(inputs$TrndCycF == "Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
#RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- eval(parse(text = inputs$FWError.FUN))
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
} #Autoc
if(inputs$TrndCycF == "Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
#EEH: In old code this was GammaFlowA which meant the global was overwritten
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF == "Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0){ RanFlow <- 0 }
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC3")
if(inputs$BSRa * Escpmnt > SRb){ r <- SRb * FWS
}else{ r <- inputs$BSRa * Escpmnt * FWS }
if(inputs$SRType == "RIC3"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * FWS
}
if(inputs$SRType == "BEV3"){
r <- 0;
if(Escpmnt > 0) r <- FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * repvars$LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
}#end 3 parameter cases
if(any(inputs$SRType==c("HOC4", "RIC4", "BEV4"))){
#Set Marine Surv Params
if(inputs$TrndCycM == "Autoc"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
#RanMarine <- GammaSample(inputs$GammaMarA, inputs$GammaMarB)
RanMarine <- eval(parse(text = inputs$MarError.FUN))
RanMarine <- Autocorrel(inputs$TCM1, LastRanMarine, RanMarine)
}
}
if(inputs$TrndCycM == "Trend"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
Mean <- Trend(inputs$TCM1, Year, inputs$MarAve, inputs$TCM2)
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- Mean
}else{
GammaMarTrA <- Mean * Mean / Var
GammaMarTrB <- Var / Mean
#RanMarine <- GammaSample(GammaMarTrA, GammaMarTrB)
RanMarine <- rgamma(1, GammaMarTrA, scale=GammaMarTrB)
}
}
}
if(inputs$TrndCycM == "Cycle"){
Mean <- Cycle(inputs$TCM1, inputs$TCM2, inputs$TCM3, Year)
# should use TCF4 as scalar in place of MarAve
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- inputs$MarAve
}else{
GammaMarCyA <- Mean * Mean / Var
GammaMarCyB <- Var / Mean
#RanMarine <- GammaSample(GammaMarCyA, GammaMarCyB)
RanMarine <- rgamma(1, GammaMarCyA, scale=GammaMarCyB)
}
}
if(RanMarine < 0) RanMarine <- 0
LastRanMarine <- RanMarine
MS <- RanMarine^inputs$BSRc
if(MS < 0) MS <- 0
#Set the flow parameters
if(inputs$TrndCycF =="Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
}
if(inputs$TrndCycF =="Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF =="Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- inputs$FlowAve
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0) RanFlow <- 0
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC4"){
if(inputs$BSRa * Escpmnt > SRb){
r=SRb * MS * FWS
}else{ r <- inputs$BSRa * Escpmnt * MS * FWS }
}
if(inputs$SRType == "RIC4"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * MS * FWS
}
if(inputs$SRType == "BEV4"){
r <- 0;
if(Escpmnt > 0) r <- MS * FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
} #end param 4 cases
if(inputs$depen == "YES" & Escpmnt < inputs$DL1 + 1){
if(Escpmnt < inputs$DL2 + 1){
fac <- Escpmnt / inputs$DL2 * inputs$DR
}else{
fac <- ((1 - (inputs$DL1 - Escpmnt) / (inputs$DL1 - inputs$DL2)) * (1 - inputs$DR)) + inputs$DR }
if(fac < 0) fac <- 0
AEQRecruits <- fac * AEQRecruits
}
# COMPUTE RECRUITS that are age 1;
# RecruitsAtAge1 is from Recruits() and is "total fraction of age 1 ind that eventually return"
# AEQRecruits/RecruitsAtAge1 <- Age 1 or Cohort[1]
CohortAge1 <- AEQRecruits / staticvars$RecruitsAtAge1
# ADD MARINE SURVIVAL IF STOCK RECRUIT FUNCTION IS SPAWNER TO SMOLT
if(inputs$MarSurv == "YES"){
BetaVariate <- CompBetaVariate(inputs$BetaMarA, inputs$BetaMarB)
CohortAge1 <- CohortAge1 * BetaVariate
}
repvars$Cohort[1] <- CohortAge1
repvars$LastRanFlow <- LastRanFlow
repvars$LastRanError <- LastRanError
repvars$LastRanMarine <- LastRanMarine
return(repvars)
}#END CompRecruits2
#' @title Title
#'
#' @param BufNum
#' @param inputs
#' @param BufSRb
#' @param YearStats
#' @param SummaryStats
#'
#' @return
#' @keywords internal
#'
#' @examples
CompStatsEEH <- function(BufNum, inputs, BufSRb, YearStats, SummaryStats){
#BufNum is BufMax, an integer
NYears <- inputs$NYears
NRuns <- inputs$NRuns
SRErrorA <- inputs$SRErrorA
SRErrorA <- inputs$SRErrorB
SRa <- inputs$BSRa
if(inputs$StepFunc=="POP"){
SRb <- BufSRb
}else{
SRB <- inputs$BSRb
}
SRc <- inputs$BSRc
SRd <- inputs$BSRd
SurvScale <- inputs$SurvScale
#In VB code Escapmnt (over ages) is redefined to be TotEscpmnt (sum over ages)
#kept terminology closer to what is being donw
TotEscpmnt <- YearStats$TotEscpmnt #a matrix of Escpmnt for each year
# COMPUTE STATISTICS
# Sum up over all reps and divide by reps
# COMPUTE AVERAGE MORTALITY OVER REPETITIONS
#divide by NRuns because we want the average over all runs
#Note in VB code, AEQMort is initially by age and then in SaveYearData, it is reassigned to be TotAEQMort.
#Here I keep AEQMort and TotAEQMort separate
SummaryStats$AvgAEQMort[BufNum] <- SummaryStats$AvgAEQMort[BufNum]+mean(YearStats$TotAEQMort)/inputs$NRuns
# COMPUTE AVERAGE CALENDAR YEAR HARVEST RATE
#CalendarHR <- YearStats$CalendarHR
#divide by NRuns because we want the average over all runs
SummaryStats$AvgCaHR[BufNum] <- SummaryStats$AvgCaHR[BufNum] + mean(YearStats$CalendarHR)/inputs$NRuns
# COMPUTE PROPORTION OF TIMES ESCAPEMENT WAS LESS THAN LOWER LEVEL
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$AvgECrit[BufNum] <- SummaryStats$AvgECrit[BufNum] + sum(TotEscpmnt < inputs$ECrit)/(inputs$NYears*inputs$NRuns)
# FIND MINIMUM, MAXIMUM, AND AVERAGE ESCAPEMENT
SummaryStats$MinEscpmnt[BufNum,] <- pmin(SummaryStats$MinEscpmnt[BufNum,],TotEscpmnt,na.rm=TRUE)
SummaryStats$MaxEscpmnt[BufNum,] <- pmax(SummaryStats$MaxEscpmnt[BufNum,],TotEscpmnt,na.rm=TRUE)
#divide these by NRuns because we want the average over all runs
SummaryStats$AvgEscpmnt[BufNum, ] <- SummaryStats$AvgEscpmnt[BufNum, ] + TotEscpmnt/inputs$NRuns
# COMPUTE PROPORTION OF TIMES ENDING ESCAPEMENT WENT BELOW QET (DL2) (depensation)
#uses last EndAv years
assessmentYears <- (inputs$NYears - (inputs$EndAv - 1)):inputs$NYears
# divide by NRuns at end to get frac of reps below
if(any(TotEscpmnt[assessmentYears] < (inputs$DL2 + 1))){
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$PropExt[BufNum] <- SummaryStats$PropExt[BufNum] + 1/inputs$NRuns
}
# COMPUTE PROPORTION OF TIMES ENDING ESCAPEMENT EXCEEDED UPPER LEVEL
EscpmntPositive <- TotEscpmnt
EscpmntPositive[EscpmntPositive<1] <- 1
#assessmentYears defined above; uses geometric mean over last EndAv years
GeomeanEscpmnt <- exp(mean(log(EscpmntPositive)[assessmentYears]))
if(GeomeanEscpmnt > inputs$ERecovery){
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$PropRec[BufNum] <- SummaryStats$PropRec[BufNum] + 1/inputs$NRuns
}
# COMPUTE BROOD YEAR HARVEST RATES
BYrAEQMort <- rep(0, inputs$NYears - (inputs$MaxAge-inputs$MinAge)) # ReDim BYrAEQMort(-1 To NYears% - 5)
BYrEscpmnt <- rep(0, inputs$NYears - (inputs$MaxAge-inputs$MinAge)) #ReDim BYrEscpmnt(-1 To NYears% - 5)
# AGE 2
Age <- 2
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,2]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,2]
# AGE 3
Age<- 3
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,3]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,3]
# AGE 4
Age <- 4
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,4]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,4]
# AGE 5
Age <- 5
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,5]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,5]
TempBYrHR <- BYrAEQMort / (BYrAEQMort + BYrEscpmnt)
TempBYrHR[BYrAEQMort < 1E-16] <- 0
#divide these by NRuns because we want the average over all runs
SummaryStats$BufAvgBYrHR[BufNum] <- SummaryStats$BufAvgBYrHR[BufNum] + mean(TempBYrHR)/inputs$NRuns
SummaryStats$AvgBYrHR[BufNum, ] <- SummaryStats$AvgBYrHR[BufNum, ] + TempBYrHR/inputs$NRuns
#don't divide these by NRuns because Max/Min
SummaryStats$MaxBYrHR[BufNum, ] <- pmax(SummaryStats$MaxBYrHR[BufNum, ], TempBYrHR, na.rm=TRUE)
SummaryStats$MinBYrHR[BufNum, ] <- pmax(SummaryStats$MinBYrHR[BufNum, ], TempBYrHR, na.rm=TRUE)
SummaryStats$AveRanMarine[BufNum] <- SummaryStats$AveRanMarine[BufNum] + mean(YearStats$RanMarine)/inputs$NRuns
SummaryStats$AveRanFlow[BufNum] <- SummaryStats$AveRanFlow[BufNum] + mean(YearStats$RanFlow)/inputs$NRuns
return(SummaryStats)
}#END CompStatsEEH
#' @title Title
#'
#' @param TotAdultEscpmnt
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
CompStockStatus <- function(TotAdultEscpmnt, inputs){
# ***** CompStockStatus *****
# Sub CompStockStatus(TotAdultEscpmnt As Double, Status%)
# Dim Break%
# Dim EscpmntCheck$
#
# Break% <- 1
# EscpmntCheck$ <- "True"
# While Break% <= NumBreakPoints% And EscpmntCheck$ <- "True"
# If TotAdultEscpmnt > EscpmntBreakPoint(Break%) Then
# EscpmntCheck$ <- "True"
# Break% <- Break% + 1
# Else
# EscpmntCheck$ <- "False"
# End If
# Wend
#
# Status% <- Break%
#
# End Sub
Break <- 1
EscpmntCheck <- "True"
while(Break <= inputs$NumBreakPoints & EscpmntCheck == "True"){
if(TotAdultEscpmnt > inputs$EscpmntBreakPoint[Break]){
EscpmntCheck <- "True"
Break <- Break + 1
}else{
EscpmntCheck <- "False"
}
}
Status <- Break
return(Status)
}#END CompStockStatus
#' @title Title
#'
#' @param a
#' @param p
#' @param s
#' @param y
#'
#' @return
#' @keywords internal
#'
#' @examples
Cycle <- function(a, p, s, y){
# ***** FUNCTION Cycle *****
# Compute cyclic variable, a <- amplitude, p <- period, s <- start, y <- year, x <- mean value of variable
# NJS: created 7/9/02 corrected 9/16/03
#
# Function Cycle(a As Double, p As Double, s As Double, y, x As Double) As Double
# a is amplitude, p is period, s is starting point, y time period
# what is x doing here? It is average value and is not needed here.
# Dim cy As Double
# cy <- Sin(2# * 3.141592654 * (y + s - 1) / p)
# 'in good survival, cycle ranges from 1 to a (amplitude)
# in bad survival, cycle ranges from 1/a to 1 (this might be lower than expected)
# If cy >= 0 Then
# cy <- (cy * (a - 1)) + 1
# Else
# cy <- (cy * (1 - (1 / a))) + 1
# End If
# Cycle <- cy + x ' use if x is changed to scalar
# Cycle <- cy
#
# End Function
# a is amplitude, p is period, s is starting point, y time period
# EEH: changed to allow vectors of y
cy <- sin(2 * pi * (y + s - 1) / p)
# in good survival, cycle ranges from 1 to a (amplitude)
# in bad survival, cycle ranges from 1/a to 1 (this might be lower than expected)
newcy <- cy
newcy[cy >= 0]<- (cy[cy >= 0] * (a - 1)) + 1
newcy[cy < 0] <- (cy[cy < 0] * (1 - (1 / a))) + 1
return(newcy)
}#END Cycle
#' @title Title
#'
#' @param Alpha
#' @param Beta
#'
#' @return
#' @keywords internal
#'
#' @examples
GammaSample <- function(Alpha, Beta){
# ***** GammmaSample *****
# Function generates a random gamma deviate with shape
# paramater alpha and scale parameter beta
# ******************************************************************
# if(get("CONSTRUN", pkg_globals)){ #no random variables
# val <- Alpha * Beta # expected value
# }else{
val <- rgamma(1, Alpha, scale=Beta)
# }
return(val)
}#END GammaSample
#' @title Title
#'
#' @param InFile
#'
#' @return
#' @keywords internal
#'
#' @examples
GetInput <- function(InFile){
# ***** GetInput *****
# Read in a .rav file and assign all the variables
# Returns the list of all inputs
#
#The rav file has , as end of input/separator
inputs <- list()
inputs$InFile <- InFile
readit <- function(skip, n){ read.table(InFile,nrows=1,sep=",",stringsAsFactors=FALSE,skip=skip)[1,n] }
# GET TITLE FOR RUN
inputs$Title <- readit(0,1) #line 1
# --------------- RUN PARAMETERS SECTION
# INPUT RANDOM NUMBER SEED, NUMBER OF CYCLES AND REPETITIONS,
# MINIMUM AND MAXIMUM AGE
inputs$RanSeed <- readit(1,1) #line 2, RanSeed
inputs$NRuns <- readit(2,1) #line 2, NRuns
inputs$NYears <- readit(3,1) #line 3, NYears
inputs$MinAge <- readit(4,1) #line 4, NYears
inputs$MaxAge <- readit(4,2) #line 4, NYears
inputs$ConvergeCrit <- readit(5,1) #line 5, ConvergeCrit
inputs$Debugg <- readit(6,1) #line 6, Debugg
inputs$Debugg <- toupper(inputs$Debugg)
if(!(inputs$Debugg %in% c("YES", "NO"))) stop("Unknown debug selection (yes/no only)")
# ----- END OF RUN PARAMETERS SECTION -
# -- STOCK-RECRUIT SECTION ---
# 'GET FORM OF SPAWNER RECRUIT FUNCTION AND PARAMETERS
# njs updated these explanations
# BSRa <- productivity parameter
# BSRb <- density dependent prarameter
# BSRc <- marine survival paramater - M^c
# BSRd <- freshwater survival parameter - exp(dF)(d should be entered as negative)
# HOC2 - Hockey stick R=Min(aS,b) a <- producitvity b=MaxRecruits
# HOC3 - R <- Min(aS,) exp(dF) (freshwater index - may be used to predict smolts)
# HOC4 - R<- Min(aS,) M^c exp(dF)
# RIC2 - Ricker R=aS exp(-bS) a <- productvity b=1/capacity
# RIC3 - R<- aS exp(-bS+dF) (freshwater index - may be used to predict smolts)
# RIC4 - Ricker with marine survival and freshwater survival
# R=aS M^c exp(-bS+dF)
# BEV2 - Beverton-Holt R=1/(b+a/S) a <- 1/productivity b=1/MaxRecruits
# BEV3 - R=1/(b+a/S) exp(dF)(freshwater index - may be used to predict smolts)
# BEV4 - BH with marine survival and freshwater survival
# R=1/(a+b/S)M^c exp(dF)
inputs$SRType <- readit(7,1) #line 8, SRType
inputs$SRType <- toupper(inputs$SRType)
SRType=inputs$SRType
allowedSRType <- c("HOC2", "HOC3", "HOC4", "BEV2", "BEV3", "BEV4", "RIC2", "RIC3", "RIC4")
if(!(SRType %in% allowedSRType)) stop("Unknown stock recruit type")
if(SRType %in% c("HOC2", "RIC2", "BEV2")){
inputs$BSRa <- readit(8,1) #line 4, NYears
inputs$BSRb <- readit(8,2) #line 4, NYears
inputs$BSRc <- 0;
inputs$BSRd <- 0;
inputs$AveEnv <- 1
}
#Skip the next 6 lines and jump to line 16
if(SRType %in% c("HOC3", "RIC3", "BEV3")){
inputs$BSRa <- readit(8,1)
inputs$BSRb <- readit(8,2)
inputs$BSRc <- 0;
inputs$BSRd <- readit(8,3)
#skip next 3 lines
inputs$FlowAve <- readit(12,1)
inputs$FlowCV <- readit(12,2)
inputs$TrndCycF <- readit(13,1)
if(!(inputs$TrndCycF %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for Flow")
inputs$TCF1 <- readit(14,1)
inputs$TCF2 <- readit(14,2)
inputs$TCF3 <- readit(14,3)
inputs$BSRc <- 0
inputs$FlowSD <- inputs$FlowCV * inputs$FlowAve
if(inputs$FlowSD > 0){
inputs$GammaFlowA <- inputs$FlowAve * inputs$FlowAve / (inputs$FlowSD * inputs$FlowSD)
inputs$GammaFlowB <- (inputs$FlowSD * inputs$FlowSD) / inputs$FlowAve
}else{
inputs$GammaFlowA <- 0
inputs$GammaFlowB <- 0
}
inputs$AveEnv <- exp(inputs$BSRd * inputs$FlowAve)
}
if(SRType %in% c("HOC4", "RIC4", "BEV4")){
inputs$BSRa <- readit(8,1)
inputs$BSRb <- readit(8,2)
inputs$BSRc <- readit(8,3);
inputs$BSRd <- readit(8,4)
inputs$MarAve <- readit(9,1)
inputs$MarCV <- readit(9,2)
inputs$TrndCycM <- readit(10,1)
if(!(inputs$TrndCycM %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for marine survival")
inputs$TCM1 <- readit(11,1)
inputs$TCM2 <- readit(11,2)
inputs$TCM3 <- readit(11,3)
inputs$FlowAve <- readit(12,1)
inputs$FlowCV <- readit(12,2)
inputs$TrndCycF <- readit(13,1)
if(!(inputs$TrndCycF %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for Flow")
inputs$TCF1 <- readit(14,1)
inputs$TCF2 <- readit(14,2)
inputs$TCF3 <- readit(14,3)
inputs$MarSD <- inputs$MarCV * inputs$MarAve
if(inputs$MarSD > 0){
inputs$GammaMarA <- inputs$MarAve * inputs$MarAve / (inputs$MarSD * inputs$MarSD)
inputs$GammaMarB <- (inputs$MarSD * inputs$MarSD) / inputs$MarAve
}else{
inputs$GammaMarA <- 0
inputs$GammaMarB <- 0
}
inputs$FlowSD <- inputs$FlowCV * inputs$FlowAve
if(inputs$FlowSD > 0){
inputs$GammaFlowA <- inputs$FlowAve * inputs$FlowAve / (inputs$FlowSD * inputs$FlowSD)
inputs$GammaFlowB <- (inputs$FlowSD * inputs$FlowSD) / inputs$FlowAve
}else{
inputs$GammaFlowA <- 0
inputs$GammaFlowB <- 0
}
inputs$AveEnv <- inputs$MarAve^inputs$BSRc * exp(inputs$BSRd * inputs$FlowAve)
}
# depensation with DL1=depensation esc; DL2 <- QET; DR % of predicted R realized at QET
inputs$depen <- readit(15,1)
inputs$depen <- toupper(inputs$depen)
if(!(inputs$depen %in% c("YES", "NO"))) stop("Unknown depensation selection (yes/no only)")
inputs$DL1 <- readit(16,1)
inputs$DL2 <- readit(16,2)
inputs$DR <- readit(16,3)
inputs$EscChoice <- readit(17,1)
inputs$EscChoice <- toupper(inputs$EscChoice)
if(!(inputs$EscChoice %in% c("YES", "NO"))) stop("Unknown escapement choice selection (yes/no only)")
# GET PARAMETERS TO ADD VARIABILITY TO STOCK RECRUIT FUNCTION
# SRErrorA and B are used for Hoc2, Ric2 and Bev2, SRErrorA, SRErrorB and ResCorPar (optional)
# are used for Hoc3, Ric3, Bev3, Hoc4, Ric4, and Bev4
inputs$SurvScale <- readit(18,1)
inputs$SurvScale <- toupper(inputs$SurvScale)
if(!(inputs$SurvScale %in% c("YES", "NO"))) stop("Unknown SR variability selection (yes/no only)")
if(inputs$SurvScale == "YES"){
inputs$SRErrorA <- readit(19,1)
inputs$SRErrorB <- readit(19,2)
inputs$ResCorParameter <- readit(19,3)
}
# GET PARAMETERS FOR SMOLT TO ADULT (MARINE) SURVIVAL IF STOCK RECRUIT FUNCTION
# IS FOR SMOLTS FROM SPAWNERS (FRESHWATER PRODUCTION)
inputs$MarSurv <- readit(20,1)
inputs$MarSurv <- toupper(inputs$MarSurv)
if(!(inputs$MarSurv %in% c("YES", "NO"))) stop("Unknown marine survival selection (yes/no only)")
if(inputs$MarSurv == "YES"){
inputs$BetaMarA <- readit(21,1)
inputs$BetaMarB <- readit(21,2)
inputs$CorrMS <- readit(21,3)
}
# - END OF STOCK RECRUIT SECTION ---
# FISHERY MANAGEMENT PARAMETERS --
# INPUT THE NUMBER OF BREAKPOINTS AND DIMENSION ARRAYS
inputs$NumBreakPoints <- readit(22,1)
# IDENTIFY WHICH LEVEL IS THE BASE REGIME
inputs$BaseRegime <- readit(23,1)
inputs$EscpmntBreakPoint <- rep(NA,inputs$NumBreakPoints+1)
inputs$TargetU <- rep(NA,inputs$NumBreakPoints+1)
#begin dynamically keeping track of line number
cline=24
# INPUT BREAKPOINTS AND TARGET EXPLOITATION RATES
for(BreakPoint in 1:(inputs$NumBreakPoints + 1)){
if(BreakPoint <= inputs$NumBreakPoints){
inputs$EscpmntBreakPoint[BreakPoint] <- readit(cline,1)
inputs$TargetU[BreakPoint] <- readit(cline,2)
cline=cline+1
}else{
inputs$TargetU[BreakPoint] <- readit(cline,1)
cline=cline+1
}
}
# INPUT PARAMETERS FOR MANAGEMENT ERROR
inputs$MgmtError <- readit(cline,1); cline=cline+1
inputs$MgmtError <- toupper(inputs$MgmtError)
if(!(inputs$MgmtError %in% c("YES", "NO"))) stop("Unknown manag error selection (yes/no only)")
# Read in dummy values if management error is No
if(inputs$MgmtError == "YES"){
inputs$GammaMgmtA <- readit(cline,1);
inputs$GammaMgmtB <- readit(cline,2); cline=cline+1
}else{
inputs$GammaMgmtA <- 0
inputs$GammaMgmtB <- 0
#still need to iterate by a step so the import isn't out of sync.
cline=cline+1
}
# ------- END OF FISHERY MANAGMENT INPUTS
inputs$ECrit <- readit(cline,1); cline=cline+1
inputs$ERecovery <- readit(cline,1);
inputs$EndAv <- readit(cline,2); cline=cline+1
# INPUT STEP SIZE AND RANGE FOR TARGET EXPLOITATION RATES OR STARTING ESCAPEMENT
# This program outputs information for a range of either exploitation rates or
# starting escapement levels.
# The BUFFER levels input below represent percentages of the base target
# exploitation rate or the starting escapement level.
# (e.g., a Buffer of .05 means that the ER will be 5% of the base target rate.)
# You may enter a number greater than 1.0, which means that the ER will be
# greater than the base target rate. Under normal runs with 0 breakpoints,
# the base target rate is not used other than to determine start, end, and step
# levels.
inputs$StepFunc <- readit(cline,1); cline=cline+1
inputs$StepFunc <- toupper(inputs$StepFunc)
if(!(inputs$StepFunc %in% c("POP", "ER"))) stop("Unknown step selection")
inputs$BufferStep <- readit(cline,1); cline=cline+1
inputs$BufferStart <- readit(cline,1);
inputs$BufferEnd <- readit(cline,2); cline=cline+1
#integer of the number of steps of ER or Pop Capacity to do # is 1:BufMax
inputs$BufMax <- round((inputs$BufferEnd - inputs$BufferStart) / inputs$BufferStep + 1)
# INPUT LOWER AND UPPER ESCAPEMENT TEST LEVELS
# The LOWER ESCAPEMENT LEVEL (ECrit) is the escapement level used by the
# program to test how often the observed escapements fall below this level.
# It may represent a "critical" level below which the spawner-recruit function
# destabilizes, and the stock increases risk of extinction, or it could be any
# level one just wanted to monitor frequency of achieving less than or equal
# to that level.
# The UPPER ECSAPEMENT LEVEL (ERecovery) is the escapement level used to compare
# against the geometric mean of the last n (EndAv) years of the run. It may be a
# management escapement goal, an interim recovery level, or some other target level.
# ------ BASE STOCK DATA SECTION
#
# INPUT INITIAL POPULATION SIZE BY AGE
# The initial population size by age, for all ages, is used to seed the model.
# In year 1, the management actions are applied to this population, and a portion of
# each size class escapes. By running the program over a range of starting population
# sizes, one can determine what minimum population size (or escapement) is needed
# to guarentee a probability of population viability.
# NJS Although the input allows for different MaxAge% from 5, the internal workings of the
# are fixed to 5 in many places.
inputs$CohortStart <- rep(0,inputs$MaxAge)
for(Age in (inputs$MinAge - 1):inputs$MaxAge){
inputs$CohortStart[Age] <- readit(cline,1); cline=cline+1
}
# INPUT NATURAL MORTALITY RATES
inputs$NatMort <- rep(0,inputs$MaxAge)
for(Age in (inputs$MinAge - 1):inputs$MaxAge){
inputs$NatMort[Age] <- readit(cline,1); cline=cline+1
}
# INPUT MATURATION RATES BY AGE (AEQ will be calculated)
inputs$MatRate <- rep(0,inputs$MaxAge)
for(Age in inputs$MinAge:inputs$MaxAge){
inputs$MatRate[Age] <- readit(cline,1); cline=cline+1
}
# INPUT PRETERMINAL AND MATURE EXPLOITATION RATES BY AGE
# These will be used to proportion target ER by age and fishery
# by def, PTU and MatU are 0 before MinAge
inputs$PTU <- rep(0,inputs$MaxAge)
inputs$MatU <- rep(0,inputs$MaxAge)
for(Age in inputs$MinAge:inputs$MaxAge){
inputs$PTU[Age] <- readit(cline,1);
inputs$MatU[Age] <- readit(cline,2); cline=cline+1
}
return(inputs)
}#END GetInput
#' @title Title
#'
#' @param InFile
#' @param OutFileBase
#' @param NRuns
#' @param silent
#' @param ...
#'
#' @return
#' @keywords internal
#'
#' @examples
Main <- function(InFile=NULL, OutFileBase=NULL, NRuns=-1, silent=FALSE,...){
#require(stringr)
#if not called with input file, then user is prompted to input one
if(is.null(InFile)) InFile <- file.choose()
if(!file.exists(InFile)) stop("Specified input file does not exist.")
if(is.null(OutFileBase)){
OutFileBase <- tools::file_path_sans_ext(basename(InFile))
#excluding this should allow the exclustion of stringr package:
# tmp <- strsplit(InFile,"\\\\")[[1]]
# InFileBase=tmp[length(tmp)]
# tmp <- strsplit(InFileBase,"/")[[1]]
# InFileBase <- tmp[length(tmp)]
# if(str_detect(InFileBase,"[.]")){
# OutFileBase <- strsplit(InFileBase,"[.]")[[1]][1];
# }else{
# OutFileBase <- InFileBase;
# }
}
#Two lists will be passed in and out of functions
# inputs <- list() #is everything from the .rav file
# staticvars <- list() #is anything computed from that; static
# READ INPUT DATA AND CALCULATE AEQs
# direct from .rav file or simple calculation from rav file inputs
#inputs <- GetInput(InFile)
time.start <- Sys.time()
cat("\n", InFile, time.start, "\n")
inputs.list <- read_rav2(InFile)
#input.vars has just the variable values and none of the metadata
inputs <- inputs.list$inputs.vars
if(NRuns>0) inputs$NRuns <- NRuns;
#add the output file names to the inputs
inputs <- SetOutFileNames(OutFileBase, inputs)
out <- RunSims(inputs, silent,...)
# SAVE SUMMARY RESULTS .sum
if(!silent) cat("\nSaving summary...\n")
data.summary <- SaveSummary(out$inputs, out$SummaryStats, out$staticvars)
# SAVE ESCAPEMENT DATA .esc
if(!silent) cat("Saving escapement data...\n")
data.escape.avg <- SaveEscpmntData(out$inputs, out$SummaryStats)
# SAVE BROOD YEAR EXPLOITATION RATE DATA .byr
if(!silent) cat("Saving BYr year data...\n")
data.byer.avg <- SaveBYrData(out$inputs, out$SummaryStats)
vrap.output <- list(list(inputs=inputs, output=out, data.summary=data.summary, data.escape.avg=data.escape.avg, data.byer.avg=data.byer.avg ))
names(vrap.output) <- tools::file_path_sans_ext(InFile)
saveRDS(vrap.output,file = paste0(tools::file_path_sans_ext(InFile), ".rds"))
# return(list(inputs=inputs, output=out, data.summary=data.summary, data.escape.avg=data.escape.avg, data.byer.avg=data.byer.avg ))
cat(Sys.time()-time.start,"\n")
}#END Main
#' @title Title
#'
#' @param prop
#' @param prev
#'
#' @return
#' @keywords internal
#'
#' @examples
progressBar <- function(prop=0, prev=0){
if (prev < 50) {
if (prop > 1) {
prop <- prop/100
}
bars <- round(prop * 50, digits <- 0) - prev
prev <- bars + prev
if (prop == 0) {
cat(" |2% |20% |40% |60% |80% |100%\n",
"Progress: ", sep <- "")
}
else {
while (bars > 0) {
cat("|")
bars <- bars - 1
}
}
if (prev == 50) {
cat("\n")
}
return(prev)
}
prev
}#END progressBar
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
Recruits <- function(inputs){
# ***** Recruits *****
# Function Recruits() As Double
# Compute factor to convert calculated spawner equivalent
# production to age cohort (source is PSC Chinook Model).
# #
# Tmp <- 0
# X9 <- 1 - NatMort(1)
# For Age% <- MinAge% To MaxAge%
# X9 <- X9 * (1 - NatMort(Age%))
# Tmp <- Tmp + X9 * MatRate(Age%)
# X9 <- X9 * (1 - MatRate(Age%))
# Next Age%
# Recruits <- Tmp
#
# End Function
# The SR function gets us the AEQRecruits from spawners in year t.
# We needs to then translate that to age 1 indiv in pop (Cohort[1])
# We know AEQRecruit. How many Age 1 individuals does that translate to? Age1 * (1- total fraction lost) <- AEQRecruits
# So Age1 <- AEQRecruits/(1-total fraction lost)
# Tmp here is total fraction of age 1 ind that eventually return
# AEQRecruits/Tmp <- Age 1 or Cohort[1]
Tmp <- 0
X9 <- 1 - inputs$NatMort[1] #survival ag 1
for(Age in inputs$MinAge:inputs$MaxAge){
X9 <- X9 * (1 - inputs$NatMort[Age])
Tmp <- Tmp + X9 * inputs$MatRate[Age]
X9 <- X9 * (1 - inputs$MatRate[Age])
}
#AEQRecruits / compvars$RecruitsAtAge1
return(Tmp) #Recruits at age 1
}#END Recruits
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
RepInit <- function(inputs){
rtn.list <- list()
#Set default values used when RanFlow and RanMarine are used (2 and 3 param SR models)
rtn.list$LastRanFlow <- 0
rtn.list$LastRanMarine <- 0
rtn.list$Cohort <- inputs$CohortStart
# SET INITIAL SEED FOR AUTOCORRELATED RESIDUALS
if(inputs$SurvScale == "YES"){
rtn.list$LastRanError <- rnorm(1, 0, sd=sqrt(inputs$SRErrorB))
}
# SET INITIAL VALUES FOR LAST RANMARINE AND RANFLOW
# EEH: VB code was misisng HOC3 and BEV3 from this
if(any(inputs$SRType==c("HOC3", "BEV3", "RIC3"))){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
rtn.list$LastRanFlow <- inputs$FlowAve
}else{
if(inputs$TCF2 > 0){ rtn.list$LastRanFlow <- inputs$TCF2;
}else{ rtn.list$LastRanFlow <- rgamma(1,inputs$GammaFlowA, scale=inputs$GammaFlowB); }
}
}
if(any(inputs$SRType==c("HOC4", "BEV4", "RIC4"))){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
rtn.list$LastRanMarine <- inputs$MarAve
}else{
if(inputs$TCM2 > 0){ rtn.list$LastRanMarine <- inputs$TCM2
}else{ rtn.list$LastRanMarine <- rgamma(1,inputs$GammaMarA, scale=inputs$GammaMarB) }
}
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
rtn.list$LastRanFlow <- inputs$FlowAve
}else{
if(inputs$TCF2 > 0){ rtn.list$LastRanFlow <- inputs$TCF2;
}else{ rtn.list$LastRanFlow <- rgamma(1,inputs$GammaFlowA, scale=inputs$GammaFlowB); }
}
}
return(rtn.list)
}#END RepInit
#' @title Title
#'
#' @param inputs
#' @param silent
#'
#' @return
#' @keywords internal
#'
#' @examples
RunSims <- function(inputs, silent){
# RunSims takes the input list and runs the VRAP functions
#Set up list that will hold static computed variables. These don't change with each rep or buffer
staticvars <- list()
# CALCULATE AEQs
# staticvars are computed variables; they are static
staticvars$AEQ <- AEQcalc(inputs)
# COMPUTE FACTOR TO TRANSLATE AEQ RECRUITMENT TO AGE 1
staticvars$RecruitsAtAge1 <- Recruits(inputs)
ptm <- proc.time()
# PROGRAM EXECUTION SECTION
Buffer <- inputs$BufferStart
#Set up the SummaryStats
SummaryStats <- SetupSummaryStats(inputs)
#Set up the YearStats
YearStats <- list()
YearStats$AEQMort <- matrix(0, inputs$NYears, inputs$MaxAge)
YearStats$Escpmnt <- matrix(0, inputs$NYears, inputs$MaxAge)
YearStats$TotAEQMort <- matrix(NA, inputs$NYears)
YearStats$TotAdultEscpmnt <- matrix(NA, inputs$NYears)
YearStats$TotEscpmnt <- matrix(NA, inputs$NYears)
YearStats$RanFlow <- matrix(NA, inputs$NYears)
YearStats$RanMarine <- matrix(NA, inputs$NYears)
#MF addition:
#inputs$output.raw is a string of the variable names to be saved for output
#if some variables are named in output.raw then build rawoutput to store them:
if(exists("output.raw", where=inputs)){
if(!is.na(inputs$output.raw)){ rawoutput <- vector(mode = "list", length = inputs$BufMax*inputs$NRuns)}
}
YearStats$Recruits <- rep(NA, inputs$NYears)
YearStats$AEQRecruits <- rep(NA, inputs$NYears)
YearStats$TargetUError <- rep(NA, inputs$NYears)
YearStats$ptu.scaled <- matrix(NA, inputs$NYears, 5)
YearStats$MatU.scaled <- matrix(NA, inputs$NYears, 5)
YearStats$Cohort <- matrix(NA, inputs$NYears, 5)
YearStats$tmpdat <- rep(NA, inputs$NYears)
if(!silent) prev <- progressBar()
# For each ER or Pop Cap level, go loop through NRuns, and for each NRun, loop through Year
for(BufNum in 1:inputs$BufMax){
# INITIALIZE BUFFER SPECIFIC PARAMETERS AND ARRAYS
out <- BufferInit(Buffer, inputs)
inputs$DR <- out$DR
#This is the ER or Pop Cap to use for this loop
#used in CompEscpmnt
BufTargetU <- out$BufTargetU
#This is the SRb to use when StepFunc <- "POP"; not used when StepFunc <- "ER"
#used in CompRecruits and CompStats
BufSRb <- out$BufSRb
# REPETITION LOOP
for(Rep in 1:inputs$NRuns){
if(!silent) prev <- progressBar((inputs$NRuns*(BufNum-1)+Rep)/(inputs$BufMax*inputs$NRuns),prev)
# INITIALIZE REPETITION SPECIFIC PARAMETERS AND ARRAYS
# repvar is a list of variables that change each year: Cohort, LastRanError, LastRanFlow, LastRanMarine
# These are updated with each year
# starts Cohort at the init # at each age in the input file
# Cohort seems to mean population #s at each age
# YearStats is a list of variables that are saved for each year
# SummaryStats passed in to fix first year Ran to be same across reps
repvars <- RepInit(inputs)
if(Rep == 1){
SummaryStats$FirstRanMarine <- repvars$LastRanMarine
SummaryStats$FirstRanFlow <- repvars$LastRanFlow
}
# BEGIN Year LOOP
for(Year in 1:inputs$NYears){
#Save this for Summary at end
YearStats$RanMarine[Year] <- repvars$LastRanMarine
YearStats$RanFlow[Year] <- repvars$LastRanFlow
# NATURAL MORTALITY PROCESS
# Update Cohort (pop size) at Age with Natural Mortality
repvars$Cohort <- CompNatMort(inputs, repvars$Cohort)
# COMPUTE ESCAPEMENT USING BASE LEVEL EXPLOITATION RATE
#Escapment is spawners
# cat("\n",Rep, class(staticvars$AEQ), Year, BufTargetU)
YearStats <- CompEscpmnt(inputs$BaseRegime, Year, inputs, BufTargetU, repvars$Cohort, staticvars$AEQ, YearStats)
YearStats$TempCohort <- round(YearStats$TempCohort,0)
repvars$TempCohort <- YearStats$TempCohort
#commented out since I want to test without the regime feature first
# # CHECK STOCK STATUS
# # Going to be a number between 1 and NumBreaks
# Status <- CompStockStatus(YearStats$TotAdultEscpmnt, inputs)
#
# # ADJUST REGIME IF WARRANTED
# if(Status !<- inputs$BaseRegime){
# YearStats=CompEscpmnt(Status, Year, inputs, BufTargetU, repvars$Cohort, staticvars$AEQ, YearStats)
# }
# AGE COHORT
# now update cohort (non-spawner pop size at age) for next year using the updated "cohort" from CompEscpmnt
repvars$Cohort <- CompAgeCohort(repvars$TempCohort, repvars$Cohort, inputs)
# COMPUTE RECRUITS FROM ESCAPEMENT
# Uses SR (Escpmnt to AEQRecruits) function with error added to that to get AEQRecruits from this years spawners (Escpmnt)
# Needs to then translate that to age 1 indiv in pop (Cohort[1])
# We know AEQRecruit. How many Age 1 individuals does that translate to? Age1 * (1- total fraction lost) <- AEQRecruits
# So Age1 <- AEQRecruits/(1-total fraction lost)
# Set Cohort[1], LastRanError, LastRanFlow, LastRanMarine
repvars <- CompRecruits(YearStats, Year, inputs, repvars, staticvars, BufSRb)
#repvars <- CompRecruits2(YearStats, Year, inputs, repvars, staticvars, BufSRb)
# SAVE YEAR DATA
# stores the year data that will be needed for the statistics at the end
YearStats <- SaveYearData(Year, YearStats)
#mf: I don't like this use of the term recruits. This is summing, the returns by age to get total returns for a return year:
YearStats$Recruits[Year] <- sum(repvars$Cohort, na.rm = TRUE)
YearStats$TargetU <- BufTargetU
#mf addition:
YearStats$AEQRecruits[Year] <- repvars$AEQRecruits
} #for loop for Year
#this cumulates some stats over each rep
SummaryStats<- CompStatsEEH(BufNum, inputs, BufSRb, YearStats, SummaryStats)
#mf:here is the output of all data:
if(exists("rawoutput")){
output.varnames <- unlist(strsplit(inputs$output.raw, split = " +"))
#convert appropriate data types to integer:
integer.data <- c("AEQMort", "Escpmnt", "TotAEQMort", "TotAdultEscpmnt", "TotEscpmnt", "Recruits", "AEQRecruits", "TempCohort")
for(i in names(YearStats)){
if(i %in% integer.data) storage.mode(YearStats[[i]]) <- "integer"}
YearStats$rep <- Rep
rawoutput[[(BufNum-1) *inputs$NRuns + Rep]] <- YearStats[output.varnames]
}#end if(exists(rawoutput)){
} #for loop for Rep
# next ER level or Pop Cap (SRb) level
Buffer <- Buffer + inputs$BufferStep
} #for loop for BufNum
#if not saving any detailed output
if(!exists("rawoutput"))rawoutput <- NULL
comp.time <- proc.time() - ptm
return(list(inputs=inputs, SummaryStats=SummaryStats, staticvars=staticvars, time=comp.time, rawoutput=rawoutput))
}#END RunSims
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveBYrData <- function(inputs, SummaryStats, writefile=FALSE){
file <- inputs$OutFileByr
output <- ""
BufNum <- 0
for(Buffer in seq(inputs$BufferStart,inputs$BufferEnd,inputs$BufferStep)){
BufNum <- BufNum + 1
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
for(Byr in -1:(inputs$NYears - inputs$MaxAge)){
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff,digits=2), nsmall=2,width=6));
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(Byr,nsmall=0,width=7));
#plus 2 because R indexing starts at 1
output <- c(output, format(round(SummaryStats$MinBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, format(round(SummaryStats$AvgBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, format(round(SummaryStats$MaxBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, "\r\n");
} #for loop for Byr
} #for loop for Buffer
if(writefile) cat(file=file, output)
invisible(output)
}#END SaveBYrData
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveEscpmntData <- function(inputs, SummaryStats, writefile=FALSE){
file <- inputs$OutFileEsc;
output <- ""
BufNum <- 0
for(Buffer in seq(inputs$BufferStart,inputs$BufferEnd,inputs$BufferStep)){
BufNum <- BufNum + 1
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
for(Year in 1:inputs$NYears){
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff,digits=2), nsmall=2,width=6));
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(Year,nsmall=0,width=7));
#plus 2 because R indexing starts at 1
output <- c(output, format(round(SummaryStats$MinEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, format(round(SummaryStats$MaxEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, "\r\n");
} #for loop for year
} #for loop for buffer
if(writefile) cat(file=file, output)
invisible(output)
}#END SaveEscpmntData
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param staticvars
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveSummary <- function(inputs, SummaryStats, staticvars, writefile=FALSE){
file <- inputs$OutFileSum;
output=""
# PRINT HEADER INFORMATION
output <- c(output, format(paste("RAPVIABILITY (R) Version "),width=54), "Date:", format(Sys.Date()),"\r\n");
output <- c(output, format("Title:",width=14), inputs$Title, "\r\n");
output <- c(output, format("Input File:",width=19), inputs$InFile, "\r\n");
output <- c(output, "Copy of Input File:", inputs$InFile, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Basic Simulation Input Parameters:","\r\n");
output <- c(output, " ","# of Years=", inputs$NYears, " # of Reps=", inputs$NRuns, " HR Conv.Crit=", inputs$ConvergeCrit, " Seed", inputs$RanSeed, "\r\n");
output <- c(output, " Range start ", inputs$BufferStart, " end ", inputs$BufferEnd, " by ", inputs$BufferStep, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Stock Recruit Function Input Parameters:", "\r\n");
output <- c(output, " ", "Function Type: ", inputs$SRType, "\r\n");
if(inputs$SRType == "HOC2"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=maxrecruits=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "HOC3"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)* exp(d*FWindex)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, " b=maxrecruits=", inputs$BSRb, sep=""),width=34), " d=FW parameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "HOC4"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)* exp(d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits=", inputs$BSRb, " c=MSparameter=", inputs$BSRc, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "RIC2"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=capacity=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "RIC3"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b + d*FWindex)", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=capacity=", inputs$BSRb, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "RIC4"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b + d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=capacity=", inputs$BSRb, " c=MSparameter=", inputs$BSRc, " d=FW parameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "BEV2"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] ", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=maxrecruits=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "BEV3"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] * exp(d*FWindex)", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits", inputs$BSRb, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "BEV3"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] * exp(d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits", inputs$BSRb, " c=MSparameter", inputs$BSRc, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
output <- c(output, "\r\n");
if(inputs$SurvScale == "YES"){
if(inputs$SRType %in% c("HOC2", "RIC2", "BEV2")){
output <- c(output, "Stock-Recruit Error Parameters (gamma distr.) [R=f(s)*e]:", "\r\n");
output <- c(output, format(paste(" ", "A=", inputs$SRErrorA, sep=""), width=29), " B=", inputs$SRErrorB, "\r\n");
output <- c(output, "\r\n");
}
if(inputs$SRType %in% c("HOC3", "RIC3", "BEV3", "HOC4", "RIC4", "BEV4")){
output <- c(output, "Stock-Recruit Error Parameters [R=f(S)*exp(e)]:", "\r\n");
output <- c(output, format(paste(" ", "MSE=", inputs$SRErrorB, sep=""), width=29), " ResCor=", inputs$ResCorParameter, "\r\n");
output <- c(output, "\r\n");
}
}
if(inputs$depen == "YES"){
output <- c(output, "Depensation at escap:", inputs$DL1, " QET:", inputs$DL2, "fraction of depensation at QET", inputs$DR, "\r\n");
output <- c(output, "\r\n");
}
if(inputs$MarSurv == "YES"){
if( inputs$SRType %in% c("HOC2", "RIC2", "BEV2", "HOC3", "RIC3", "BEV3")){
output <- c(output, "Smolt to Adult Survival Rate Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Beta_A=", inputs$BetaMarA, sep=""), width=34), " Beta_B=", inputs$BetaMarB, "\r\n");
output <- c(output, "\r\n");
}else{
output <- c(output, "Should not be using this for ", inputs$SRType, "\r\n");
output <- c(output, "\r\n");
}
}
BufMax <- (inputs$BufferEnd - inputs$BufferStart) / inputs$BufferStep + 1
if( inputs$SRType %in% c("HOC3", "RIC3", "BEV3")){
output <- c(output, "Freshwater Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaFlowA, sep=""), width=34), " Gamma B=", inputs$GammaFlowB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$GammaFlowA * inputs$GammaFlowB, sep=""), width=34), " var=", inputs$GammaFlowB * inputs$GammaFlowB * inputs$GammaFlowA, "\r\n");
output <- c(output, " ", inputs$TrndCycM, inputs$TCF1, inputs$TCF2, inputs$TCF3, "\r\n");
output <- c(output, " ", "First Flow index generated <- ", SummaryStats$FirstRanFlow, "\r\n");
output <- c(output, " ", "Average Flow index generated <- ", SummaryStats$AveRanFlow / inputs$NYears / inputs$NRuns / inputs$BufMax, "\r\n");
output <- c(output, "\r\n");
}
if( inputs$SRType %in% c("HOC4", "RIC4", "BEV4")){
output <- c(output, "Marine Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaMarA, sep=""), width=34), " Gamma B=", inputs$GammaMarB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$MarAve, sep=""), width=34), " st dev=", inputs$MarSD, "\r\n");
output <- c(output, " ", inputs$TrndCycM, inputs$TCM1, inputs$TCM2, inputs$TCM3, "\r\n");
output <- c(output, " ", "First Marine index generated <- ", SummaryStats$FirstRanMarine, "\r\n");
output <- c(output, " ", "Average Marine index generated <- ", mean(SummaryStats$AveRanMarine), "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Freshwater Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaFlowA, sep=""), width=34), " Gamma B=", inputs$GammaFlowB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$FlowAve, sep=""), width=34), " st dev=", inputs$FlowSD, "\r\n");
output <- c(output, " ", inputs$TrndCycF, inputs$TCF1, inputs$TCF2, inputs$TCF3, "\r\n");
output <- c(output, " ", "First Flow index generated <- ", SummaryStats$FirstRanFlow, "\r\n");
output <- c(output, " ", "Average Flow index generated <- ", mean(SummaryStats$AveRanFlow), "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, "Fishery Regime Parameters:", "\r\n");
if(inputs$NumBreakPoints >= 1){
for(Break in 1:inputs$NumBreakPoints){
output <- c(output, " ", "Breakpoint", format(Break, digits=2), "\r\n");
output <- c(output, " ", "HR Below Breakpoint=", inputs$TargetU[Break], "\r\n");
output <- c(output, " ", "Escapement=", inputs$EscpmntBreakPoint[Break], "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, " ", "HR Above BreakPoint=", inputs$TargetU[inputs$NumBreakPoints + 1], "\r\n");
output <- c(output, "\r\n");
}else{
output <- c(output, " ", "Base ER <- ", inputs$TargetU[inputs$BaseRegime], "\r\n");
}
if(inputs$MgmtError == "YES"){
output <- c(output, "Management Variability Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaMgmtA, sep=""), width=34), " Gamma B=", inputs$GammaMgmtB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$GammaMgmtA * inputs$GammaMgmtB, sep=""), width=34), " var=", inputs$GammaMgmtB * inputs$GammaMgmtB * inputs$GammaMgmtA, "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, "AEQ for age class", "\r\n");
for(Age in inputs$MinAge:inputs$MaxAge){
output <- c(output, " ", "Age", Age, " AEQ =", staticvars$AEQ[Age], "\r\n");
}
output <- c(output, "Recruits At Age 1<- ", staticvars$RecruitsAtAge1, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Regime Evaluation Parameters: QET <- ", inputs$DL2, "\r\n");
output <- c(output, " ", "Lower Escapement Level (LEL)=", inputs$ECrit, "\r\n");
output <- c(output, " ", "Upper Escapement Level (UEL)=", inputs$ERecovery, "\r\n");
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output <- c(output, " ", "Max Return (under average variability) =", inputs$BSRa * inputs$AveEnv * inputs$BSRb / 2.71828, "\r\n");
}else{
output <- c(output, " ", "Max Return (under average variability) =", inputs$BSRb * inputs$AveEnv, "\r\n");
}
output <- c(output, "\r\n");
output <- c(output, "SUMMARY STATISTICS", "\r\n");
output <- c(output, "All statistics are averaged over repetitions", "\r\n");
output <- c(output, " . . . . __________Escapement___________", "\r\n");
output <- c(output, ". b Total-Exploit.-Rate . . #fish %runs %yrs %runs 1st LastYrs pop_size.", "\r\n");
output <- c(output, ". param. TgtER CYrER BYrER Mort. extnct <LEL end>UEL Year Ave. at_equil. ", "\r\n");
Buffer <- inputs$BufferStart
for(BufNum in 1:inputs$BufMax){
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff, digits=2),nsmall=2,width=7));
output <- c(output, format(round(SummaryStats$AvgCaHR[BufNum], digits=3),nsmall=3,width=7));
output <- c(output, format(round(SummaryStats$BufAvgBYrHR[BufNum], digits=3),nsmall=3,width=7));
output <- c(output, format(round(SummaryStats$AvgAEQMort[BufNum]),nsmall=0,width=6));
output <- c(output, format(round(100 * SummaryStats$PropExt[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(100 * SummaryStats$AvgECrit[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(100 * SummaryStats$PropRec[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, 1]),nsmall=0,width=7));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, inputs$NYears]),nsmall=0,width=6));
if(inputs$StepFunc == "POP"){
# njs MxR adjusted for average environmental conditions
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
output <- c(output, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime])),nsmall=0,width=8));
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output <- c(output, format(round(PBff * inputs$BSRb * log(inputs$AveEnv * inputs$BSRa * (1 - inputs$TargetU[inputs$BaseRegime]))),nsmall=0,width=8));
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
output <- c(output, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime] - ((1 / inputs$BSRa) * inputs$AveEnv))),nsmall=0,width=8));
}
}
output <- c(output, "\r\n");
Buffer <- Buffer + inputs$BufferStep
} #for(BufNum in 1:inputs$BufMax)
### this creates a table of just the summary columns:
output.tabulated <- data.frame(beta=numeric(), TgtER=numeric(), CYrER=numeric(), BYrER=numeric(), fishMort=numeric(), percent.runs.extinct=numeric(), percent.years.belowLEL=numeric(), percent.years.aboveUEL=numeric(), first.year=numeric(), lastyears.avg=numeric())
Buffer <- inputs$BufferStart
for(BufNum in 1:inputs$BufMax){
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
output.row <- cbind(format(round(inputs$BSRb * PBff),nsmall=0,width=7),
format(round(inputs$TargetU[inputs$BaseRegime] * EBff, digits=2),nsmall=2,width=7),
format(round(SummaryStats$AvgCaHR[BufNum], digits=3),nsmall=3,width=7),
format(round(SummaryStats$BufAvgBYrHR[BufNum], digits=3),nsmall=3,width=7),
format(round(SummaryStats$AvgAEQMort[BufNum]),nsmall=0,width=6),
format(round(100 * SummaryStats$PropExt[BufNum],digits=1),nsmall=1,width=6),
format(round(100 * SummaryStats$AvgECrit[BufNum],digits=1),nsmall=1,width=6),
format(round(100 * SummaryStats$PropRec[BufNum],digits=1),nsmall=1,width=6),
format(round(SummaryStats$AvgEscpmnt[BufNum, 1]),nsmall=0,width=7),
format(round(SummaryStats$AvgEscpmnt[BufNum, inputs$NYears]),nsmall=0,width=6))
if(inputs$StepFunc == "POP"){
# njs MxR adjusted for average environmental conditions
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime])),nsmall=0,width=8))
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * log(inputs$AveEnv * inputs$BSRa * (1 - inputs$TargetU[inputs$BaseRegime]))),nsmall=0,width=8))
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime] - ((1 / inputs$BSRa) * inputs$AveEnv))),nsmall=0,width=8))
}
}#if(inputs$StepFunc == "POP"){
output.row <- as.numeric(output.row)
output.tabulated <- rbind(output.tabulated, output.row)
Buffer <- Buffer + inputs$BufferStep
} #for(BufNum in 1:inputs$BufMax)
colnames(output.tabulated) <- c("beta", "TgtER", "CYrER", "BYrER", "fishMort", "percent.runs.extinct", "percent.years.belowLEL", "percent.years.aboveUEL", "first.year", "lastyears.avg")
### end table creation
if(writefile) cat(file=file, output)
invisible(list(output=output,output.tabulated=output.tabulated ))
}#END SaveSummary
#' @title Title
#'
#' @param Year
#' @param YearStats
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveYearData <- function(Year, YearStats){
if(YearStats$TotAEQMort[Year] + YearStats$TotEscpmnt[Year] > 0.000000001){
HR <- YearStats$TotAEQMort[Year] / (YearStats$TotAEQMort[Year] + YearStats$TotEscpmnt[Year])
}else{
HR <- 0
}
#Update AEQMort to hold total
#in VB code Escpmnt (a vector for each age) is redefined to be TotEscpmnt (sum over ages)
#unclear why
# YearStats$AEQMort[Year] <- YearStats$TotAEQMort
# YearStats$Escpmnt[Year] <- YearStats$TotEscpmnt
YearStats$CalendarHR[Year] <- HR
# YearStats$Age2AEQMort[Year] <- YearStats$AEQMort[Year,2]
# YearStats$Age3AEQMort[Year] <- YearStats$AEQMort[Year,3]
# YearStats$Age4AEQMort[Year] <- YearStats$AEQMort[Year,4]
# YearStats$Age5AEQMort[Year] <- YearStats$AEQMort[Year,5]
# YearStats$Age2Escpmnt[Year] <- YearStats$Escpmnt[Year,2]
# YearStats$Age3Escpmnt[Year] <- YearStats$Escpmnt[Year,3]
# YearStats$Age4Escpmnt[Year] <- YearStats$Escpmnt[Year,4]
# YearStats$Age5Escpmnt[Year] <- YearStats$Escpmnt[Year,5]
return(YearStats)
}#END SaveYearData
#' @title Title
#'
#' @param BaseName
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
SetOutFileNames <- function(BaseName, inputs){
# ***** SetOutFileNames ******
# Used by GetOutFiles and also by GetCommandLine
# ***********************************************************************
# NoExtensionName 'stores path and name, but not extension
# Position 'position of final "\" in pathname
# PathName 'string containing the path (including last "\")
# InFileNameOnly 'name (excl path) of input file
# SET OUTPUT FILES
inputs$OutFileSum <- paste(BaseName,".sum",sep="")
inputs$OutFileEsc <- paste(BaseName,".esc",sep="")
inputs$OutFileByr <- paste(BaseName,".byr",sep="")
inputs$OutFileLog <- paste(BaseName,".log",sep="")
inputs$InFileCopy <- paste(BaseName,".rav",sep="")
return(inputs)
}#END SetOutFileNames
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
SetupSummaryStats <- function(inputs){
SummaryStats <- list()
#most get set to 0 because I add info/NRuns each rep
#mins and maxs get NA because I use min(a,b,na.rm=TRUE) to ignore initial setting
SummaryStats[["AvgEscpmnt"]] <- matrix(0, inputs$BufMax, inputs$NYears)
for(el in c("MinEscpmnt","MaxEscpmnt"))
SummaryStats[[el]] <- matrix(NA, inputs$BufMax, inputs$NYears)
SummaryStats[["AvgBYrHR"]] <- matrix(0, inputs$BufMax, inputs$NYears - (inputs$MaxAge - inputs$MinAge))
for(el in c("MinBYrHR","MaxBYrHR"))
SummaryStats[[el]] <- matrix(NA, inputs$BufMax, inputs$NYears - (inputs$MaxAge - inputs$MinAge))
for(el in c("AvgAEQMort","AvgECrit","AvgCaHR", "BufAvgBYrHR", "PropExt", "PropRec", "BufAvgBYrHR", "AveRanFlow", "AveRanMarine"))
SummaryStats[[el]] <- rep(0, inputs$BufMax)
return(SummaryStats)
}#END SetupSummaryStats
#' @title Title
#'
#' @param t
#' @param y
#' @param x
#' @param z
#'
#' @return
#' @keywords internal
#'
#' @examples
Trend <- function(t, y, x, z){
# ***** FUNCTION Trend *****
# Compute variable x with trend t for year y
# NJS: created 7/9/02
# Function Trend(t As Double, y, x As Double, z) As Double
# t is trend rate, y is time increment, x is first value, z is type of trend
#
# If z <- 0 Then
# Trend <- x * (1 + t) ^ y
# ElseIf z <- 1 Then
# Trend <- x + (y * t)
# Else
# If Trend < 0 Then Trend <- 0
#
# Print "Unknown trend/cycle function"
# Stop
# End If
#
# End Function
# t is trend rate, y is time increment, x is first value, z is type of trend 0 or 1
# EEH: changed to work with vectors of y
if(!(z == 0 | z==1)) stop("Unknown trend/cycle function")
if(z == 0){
val <- x * (1 + t) ^ y
}else{
val <- x + (y * t)
}
val[val < 0] <- 0
return(val)
}#END Trend
MichaelFolkes/VRAP2 documentation built on June 9, 2019, 6:27 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' @title Title
#'
#' @param input
#'
#' @return
#' @keywords internal
#'
#' @examples
AEQcalc <- function(input){
# ***** SUB AEQcalc *****
# Compute the AEQs for each age group
# TAB: what are AEQs? Adult Equivalents!
# Sub AEQcalc()
# Dim TmpA As Double 'TAB: added line
# Dim TmpS As Double 'TAB: added line
# Dim Age As Integer 'TAB: added line
#
# TmpA <- 0
# TmpS <- 0
# For Age <- MaxAge% To MinAge% Step -1
# AEQ(Age) <- MatRate(Age) + TmpS * (1 - MatRate(Age)) * TmpA
# TmpA <- AEQ(Age)
# TmpS <- 1 - NatMort(Age)
# Next Age
#
# End Sub
MinAge <- input$MinAge
MaxAge <- input$MaxAge
MatRate <- input$MatRate
NatMort <- input$NatMort
res <- rep(NA, 6)
AEQ <- rep(0,MaxAge);
TmpA <- 0
TmpS <- 0
for(Age in MaxAge:MinAge){
AEQ[Age] <- MatRate[Age] + TmpS * (1 - MatRate[Age]) * TmpA
# res <- rbind(res, round(c(Age, MatRate[Age], TmpS,TmpA, NatMort[Age],AEQ[Age] ),3))
TmpA <- AEQ[Age]
TmpS <- 1 - NatMort[Age]
}
# colnames(res) <- c("Age", "MatRate", "TmpS", "TmpA", "NatMort","AEQ")
# print(res[-1,])
#return(round(AEQ,3) )
return(AEQ)
}#END AEQcalc
#' @title Title
#'
#' @param input
#'
#' @return
#' @keywords internal
#'
#' @examples
AEQcalc_modified <- function(input){
MinAge <- input$MinAge
MaxAge <- input$MaxAge
MatRate <- input$MatRate
NatMort <- input$NatMort
AEQ <- rep(0,MaxAge);
ad <- rep(NA,MaxAge)
sur <- rep(NA,MaxAge)
TmpA <- 0
TmpS <- 0
res <- data.frame(TmpA=numeric(MaxAge), TmpS=numeric(MaxAge), AEQ=numeric(MaxAge), MatRate=numeric(MaxAge))
for(Age in MaxAge:MinAge){
AEQ[Age] <- MatRate[Age] + TmpS * (1 - MatRate[Age]) * TmpA
res$TmpA[Age] <- TmpA
res$TmpS[Age] <- TmpS
res$AEQ[Age] <- AEQ[Age]
res$MatRate[Age] <- MatRate[Age]
TmpA <- AEQ[Age]
TmpS <- 1 - NatMort[Age]
}
return(res)
}#END AEQcalc_modified
#' @title Title
#'
#' @param p
#' @param lastx
#' @param x
#'
#' @return
#' @keywords internal
#'
#' @examples
Autocorrel <- function(p, lastx, x){
# ***** FUNCTION Autocorrel *****
# Compute autocorrelated variable, p <- autocorrelation, lastx <- last value of variable x
# NJS: created 7/9/02
# Function Autocorrel(p As Double, lastx As Double, x As Double) As Double
#
# Autocorrel <- p * lastx + (1 - p) * x
#
# End Function
#
val <- p * lastx + (1 - p) * x;
return(val)
}#END Autocorrel
#' @title Title
#'
#' @param Buffer
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
BufferInit <- function(Buffer, inputs){
# ***** BufferInit *****
# ADJUST TARGET RATE FOR BUFFER EXPLOITATION
#These are the internal function variables
# Defines MxR, SRb, QetR
# Changes DR, an original input
rtn.list <- list()
if(inputs$StepFunc=="POP"){
PBff <- Buffer #population capacity buffer; SRb x this
EBff <- 1 #exploitation buffer; er x this
}
if(inputs$StepFunc=="ER"){
PBff <- 1 #population capacity buffer; SRb x this
EBff <- Buffer #exploitation buffer; er x this
}
BufSRb <- PBff * inputs$BSRb #adjust capacity upward
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
MxR <- BufSRb * inputs$AveEnv #never used except printed out; Max Recruits adj by average environment
QetR <- inputs$BSRa * inputs$AveEnv * inputs$DL2
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
MxR <- BufSRb * inputs$AveEnv #never used except printed out;
QetR <- inputs$AveEnv / ((1 / BufSRb) + (1 / inputs$BSRa) * (1 / inputs$DL2))
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
MxR <- inputs$BSRa * inputs$AveEnv * BufSRb / 2.71828 #never used except printed out;
QetR <- inputs$AveEnv * inputs$BSRa * inputs$DL2 * exp(-inputs$DL2 / BufSRb)
}
#cat(paste("Max Recruits", MxR))
rtn.list$DR <- inputs$DR
#DR is one of the original inputs; as is DL2; QetR is local variable
if(inputs$DR == 1){
if(QetR == 0){
rtn.list$DR <- 0
}else{
rtn.list$DR <- inputs$DL2 / QetR #TAB: don't like this, changing DR which is one of the original inputs
}
}
#NumBreakPoints is an original input
#This is just the ER to use for a particular sim (Buffer); only changed if StepFunc="ER"
rtn.list$BufTargetU <- c()
for(Break in 1:(inputs$NumBreakPoints + 1)){
if(inputs$NumBreakPoints > 1){ #we have 3 or more Buffer targets
if(Break == 1){
rtn.list$BufTargetU[1] <- min(inputs$TargetU[2] * EBff, inputs$TargetU[1]) #use the smaller of these 2
}else{
rtn.list$BufTargetU[Break] <- inputs$TargetU[Break] * EBff
}
}else{ # we have 1 or 2 Buffer targets
rtn.list$BufTargetU[Break] <- inputs$TargetU[Break] * EBff
}
}
rtn.list$BufSRb <- BufSRb #adjusted if StepFunc <- Pop
return(rtn.list)
}#END BufferInit
#' @title Title
#'
#' @param TempCohort
#' @param Cohort
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
CompAgeCohort <- function(TempCohort, Cohort, inputs){
# ***** CompAgeCohort *****
# Sub CompAgeCohort(TempCohort() As Double)
# Dim Age As Integer 'counter for ages
#
# For Age <- MaxAge% To MinAge% + 1 Step -1
# Cohort(Age) <- TempCohort(Age - 1)
# Next Age
#
# recruitment I think
# Cohort(MinAge%) <- Cohort(MinAge% - 1)
#
# End Sub
#TempCohort is Cohort - mortality and fish that return to spawn
#Cohort age t in next year is cohort age t-1 this year
# for(Age in inputs$MaxAge:(inputs$MinAge + 1)){
# Cohort[Age] <- TempCohort[Age - 1]
# }
Cohort[inputs$MaxAge:(inputs$MinAge + 1)] <- TempCohort[inputs$MaxAge:(inputs$MinAge + 1)-1]
#This 'ages' recruitment in Cohort[1]
# EEH: this only works if MinAge=2, because Cohort[1] was set in CompRecruits (not Cohort[minage-1])
# EEH: why isn't TempCohort[inputs$MinAge - 1] <- Cohort[inputs$MinAge - 1]?
Cohort[inputs$MinAge] <- Cohort[inputs$MinAge - 1]
return(Cohort)
}#END CompAgeCohort
#' @title Title
#'
#' @param Alpha
#' @param Beta
#'
#' @return
#' @keywords internal
#'
#' @examples
CompBetaVariate <- function(Alpha, Beta){
# ***** CompBetaVariate *****
# This subroutine generates a beta r.v.
# It does so by using the gamma distribution
# Mean <- alpha/(alpha+beta)
# Variance <- (a*b)/((a+b)^2*(a+b+1))
# TAB: not sure this is correct way of getting beta distr
# since Beta(a,b) <- (Gamma(a)*Gamma(b))/(Gamma(a+b))
# ***********************************************
# Function CompBetaVariate(Alpha As Double, Beta As Double)
# Dim G1 As Double
# Dim G2 As Double
#
# If CONSTRUN <- True Then
# CompBetaVariate <- Alpha / (Alpha + Beta) 'expected value
# Else
# If Alpha <= 1 Then 'TAB: changed ALP to Alpha in this line
# G1 <- Gam1(Alpha, 1)
# Else
# G1 <- Gam2(Alpha, 1)
# End If
#
# If Beta <= 1 Then 'TAB: changed BET to Beta in this line
# G2 <- Gam1(Beta, 1)
# Else
# G2 <- Gam2(Beta, 1)
# End If
#
# CompBetaVariate <- G1 / (G1 + G2)
# End If
# End Function
#EEH for testing purposes. Replace with rbeta(alpha,beta) later
#local variables
# Dim G1 As Double
# Dim G2 As Double
# if(get("CONSTRUN", pkg_globals)){ #no random variables
# val <- Alpha / (Alpha + Beta) # expected value
# }else{
G1 <- rgamma(1, Alpha, scale=1)
G2 <- rgamma(1, Beta, scale=1)
val <- G1 / (G1 + G2)
# }
return(val)
}#END CompBetaVariate
#' @title Title
#'
#' @param Regime
#' @param Year
#' @param inputs
#' @param BufTargetU
#' @param Cohort
#' @param AEQ
#' @param YearStats
#'
#' @return
#' @keywords internal
#'
#' @examples
CompEscpmnt <- function(Regime, Year, inputs, BufTargetU, Cohort, AEQ, YearStats){
#If you have breakpoints, then Regime refers to which of these to use
# INITIALIZE PARAMETERS, VARIABLES, AND ARRAYS
Converge <- "No"
MaxAge=inputs$MaxAge
# SET EXPLOITATION RATE TARGET AFTER ADJUSTMENT FOR MANAGEMENT ERROR.
# SET INITIAL FISHERY SCALE To 1.0 UNLESS TARGET U=0
ErrorScale <- 1
if(inputs$MgmtError == "YES") {
#ErrorScale <- rgamma(1, inputs$GammaMgmtA, scale=inputs$GammaMgmtB)
#keep sampling until realized ER is <1
#this is a break from the original method of converting all er>=1 to .99
#as R lacks a do loop, this does the same.
#(need to sample error before evaluating it, and 'while' can't do that)
#idea taken from:
#https://stackoverflow.com/questions/4357827/do-while-loop-in-r
repeat{
ErrorScale <- eval(parse(text=inputs$MgmtError.FUN))
if(BufTargetU[Regime] * ErrorScale< 1 ) {break} }
}#END if(inputs$MgmtError == "YES") {
if(BufTargetU[Regime] > 0){
BufTargetUError <- BufTargetU[Regime] * ErrorScale
if(BufTargetUError >= 1){
#cat("Warning - Target HR Reduced to 0.99\n")
BufTargetUError <- 0.99 #not allowed to take all fish?
}
FishScale <- 1 #there is fishing
}else{
TempCohort <- Cohort*(1-inputs$MatRate)
Escpmnt <- Cohort * inputs$MatRate
Escpmnt[Escpmnt < 1] <- 0
AEQMort <- rep(0, inputs$MaxAge)
TotAEQMort <- 0
TotEscpmnt <- sum(Escpmnt)
Converge <- "Yes" #don't need to go through the while loop
}#END if(BufTargetU[Regime] > 0){
n <- 0
while(Converge == "No"){
n <- n + 1
# COMPUTE PRETERMINAL MORTALITY AND UPDATE COHORT
tmp=FishScale * inputs$PTU
tmp[tmp>1] <- 1 #can't take more fish than are there
PTMort <- tmp * Cohort
TempCohort <- Cohort - PTMort
# for(Age in inputs$MinAge:inputs$MaxAge){
# if(FishScale * inputs$PTU[Age] < 1){ #FishScale is 1 or 0
# PTMort[Age] <- FishScale * Cohort[Age] * inputs$PTU[Age] #take a fraction of fish
# }else{
# PTMort[Age] <- Cohort[Age] #take all fish at that age
# }
# TempCohort[Age] <- Cohort[Age] - PTMort[Age] #pre-term mort
# } #end for loop age
# COMPUTE MATURE RUN AND UPDATE COHORT
MatRun <- TempCohort * inputs$MatRate
TempCohort <- TempCohort - MatRun
# for(Age in inputs$MinAge:inputs$MaxAge){
# MatRun[Age] <- TempCohort[Age] * inputs$MatRate[Age]
# TempCohort[Age] <- TempCohort[Age] - MatRun[Age]
# } #end for loop Age
# COMPUTE MATURE MORTALITY AND ESCAPEMENT
# MatU, MatMort, Escpmnt, MatRun are all 1:MaxAge vectors
tmp <- FishScale * inputs$MatU
tmp[tmp>1] <- 1 #can't take more fish than are there
MatMort <- tmp * MatRun
Escpmnt <- MatRun - MatMort #spawners
Escpmnt[Escpmnt < 1] <- 0
# for(Age in inputs$MinAge:inputs$MaxAge){
# if(FishScale * inputs$MatU[Age] < 1){
# MatMort[Age] <- FishScale * MatRun[Age] * inputs$MatU[Age] #these are the fish that return and are taken
# }else{ #all mature fish taken
# MatMort[Age] <- MatRun[Age]
# }
# Escpmnt[Age] <- MatRun[Age] - MatMort[Age] #spawners
# if(Escpmnt[Age] < 1) Escpmnt[Age] <- 0
# }# end Age for loop
# COMPUTE AEQ TOTAL MORTALITY EXPLOITATION RATE
AEQMort <- AEQ * PTMort + MatMort
TotAEQMort <- sum(AEQMort)
TotEscpmnt <- sum(Escpmnt)
# for(Age in inputs$MinAge:inputs$MaxAge){
# AEQMort[Age] <- (AEQ[Age] * PTMort[Age]) + MatMort[Age]
# TotAEQMort <- TotAEQMort + AEQMort[Age]
# TotEscpmnt <- TotEscpmnt + Escpmnt[Age]
# #!!! MinAge fish (age 2) not counted
# if(Age > inputs$MinAge) TotAdultEscpmnt <- TotAdultEscpmnt + Escpmnt[Age]
# } #end for loop
# COMPARE ACTUAL AND TARGET; SCALE EXPLOITATION RATES IF NECESSARY
# if any of these met, exit
if(n > 100 | BufTargetU[Regime] <= 0 | (TotAEQMort + TotEscpmnt) < 1){
Converge <- "Yes"
}else{ #none are met, so determine if converged
ActualU <- TotAEQMort / (TotAEQMort + TotEscpmnt)
PercentError <- abs((ActualU - BufTargetUError) / BufTargetUError)
if(PercentError > inputs$ConvergeCrit){
FishScale <- FishScale * BufTargetUError / ActualU
}else{
Converge <- "Yes"
}
}
} #end of while
#not needed inside while loop
#TotAdultEscpmnt <- sum(Escpmnt[(inputs$MinAge+1):inputs$MaxAge])
YearStats$AEQMort[Year,] <- AEQMort
YearStats$Escpmnt[Year,] <- Escpmnt
YearStats$TotAdultEscpmnt[Year] <- sum(Escpmnt[(inputs$MinAge+1):inputs$MaxAge])
YearStats$TotAEQMort[Year] <- TotAEQMort
YearStats$TotEscpmnt[Year] <- TotEscpmnt
YearStats$TempCohort <- TempCohort
YearStats$tmpdat[Year] <- ErrorScale
#mf addition of other output to YearStats:
YearStats$TargetUError[Year] <- BufTargetU[Regime] * ErrorScale
YearStats$Cohort[Year,] <- Cohort
if(exists("FishScale")){
YearStats$ptu.scaled[Year,] <- FishScale * inputs$PTU
YearStats$MatU.scaled[Year,] <- FishScale * inputs$MatU
} else{
YearStats$ptu.scaled[Year,] <- NA
YearStats$MatU.scaled[Year,] <- NA
}
return(YearStats)
}#END CompEscpmnt
#' @title Title
#'
#' @param inputs
#' @param CohortBeforeNatMort
#'
#' @return
#' @keywords internal
#'
#' @examples
CompNatMort <- function(inputs, CohortBeforeNatMort){
# ***** CompNatMort *****
# Let the number of fish in each age class decrease according to
# the natural mortality in that age class.
# Sub CompNatMort()
# Dim Age% 'counter for ages
#
# For Age% <- MinAge% - 1 To MaxAge%
# Cohort(Age%) <- Cohort(Age%) * (1 - NatMort(Age%))
# Next Age%
#
# End Sub
#updates global array Cohort
# CohortAfterNatMort <- CohortBeforeNatMort
# for(Age in (inputs$MinAge - 1):inputs$MaxAge){
# CohortAfterNatMort[Age] <- CohortBeforeNatMort[Age] * (1 - inputs$NatMort[Age])
# }
CohortAfterNatMort <- CohortBeforeNatMort * (1 - inputs$NatMort)
return(CohortAfterNatMort)
}#END CompNatMort
#' @title Title
#'
#' @param YearStats
#' @param Year
#' @param inputs
#' @param repvars
#' @param staticvars
#' @param BufSRb
#'
#' @return
#' @keywords internal
#'
#' @examples
CompRecruits <- function(YearStats, Year, inputs, repvars, staticvars, BufSRb){
#if StepFunc="Pop" then SRb (capacity) is being changed.
if(inputs$StepFunc=="POP"){
SRb=BufSRb
}else{
SRb=inputs$BSRb
}
#cummulative stats
LastRanFlow <- repvars$LastRanFlow
LastRanError <- repvars$LastRanError
LastRanMarine <- repvars$LastRanMarine
#Horrid; renaming TotEscpmnt to Escpmnt
if(inputs$EscChoice == "YES"){
Escpmnt <- YearStats$TotAdultEscpmnt[Year]
}else{ Escpmnt <- YearStats$TotEscpmnt[Year] }
if(Escpmnt < 1){ Escpmnt <- 0 }
# COMPUTE AEQ RECRUITMENT
if(toupper(inputs$SRType) == "HOC2"){
if(inputs$BSRa * Escpmnt > SRb){
AEQRecruits <- SRb
}else{
AEQRecruits <- inputs$BSRa * Escpmnt
}
}
if(toupper(inputs$SRType) == "RIC2"){
AEQRecruits <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb)
}
if(toupper(inputs$SRType) == "BEV2"){
AEQRecruits <- 0
if(Escpmnt > 0) AEQRecruits <- 1/((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
if(any(inputs$SRType==c("HOC2","RIC2","BEV2")) & inputs$SurvScale == "YES"){
#StockScalar <- GammaSample(inputs$SRErrorA, inputs$SRErrorB)
#StockScalar <- rgamma(1, inputs$SRErrorA, scale=inputs$SRErrorB)
StockScalar <- eval(parse(text=inputs$SRError.FUN))
AEQRecruits <- StockScalar * AEQRecruits
}
if(any(inputs$SRType==c("HOC3", "RIC3", "BEV3"))){
if(inputs$TrndCycF == "Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
#RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- eval(parse(text=inputs$FlowError.FUN))
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
} #Autoc
if(inputs$TrndCycF == "Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
#EEH: In old code this was GammaFlowA which meant the global was overwritten
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF == "Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0){ RanFlow <- 0 }
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC3")
if(inputs$BSRa * Escpmnt > SRb){ r <- SRb * FWS
}else{ r <- inputs$BSRa * Escpmnt * FWS }
if(inputs$SRType == "RIC3"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * FWS
}
if(inputs$SRType == "BEV3"){
r <- 0;
if(Escpmnt > 0) r <- FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * repvars$LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
}#end 3 parameter cases
if(any(inputs$SRType==c("HOC4", "RIC4", "BEV4"))){
#Set Marine Surv Params
if(inputs$TrndCycM == "Autoc"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
#RanMarine <- GammaSample(inputs$GammaMarA, inputs$GammaMarB)
RanMarine <- rgamma(1, inputs$GammaMarA, scale=inputs$GammaMarB)
RanMarine <- Autocorrel(inputs$TCM1, LastRanMarine, RanMarine)
}
}
if(inputs$TrndCycM == "Trend"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
Mean <- Trend(inputs$TCM1, Year, inputs$MarAve, inputs$TCM2)
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- Mean
}else{
GammaMarTrA <- Mean * Mean / Var
GammaMarTrB <- Var / Mean
#RanMarine <- GammaSample(GammaMarTrA, GammaMarTrB)
RanMarine <- rgamma(1, GammaMarTrA, scale=GammaMarTrB)
}
}
}
if(inputs$TrndCycM == "Cycle"){
Mean <- Cycle(inputs$TCM1, inputs$TCM2, inputs$TCM3, Year)
# should use TCF4 as scalar in place of MarAve
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- inputs$MarAve
}else{
GammaMarCyA <- Mean * Mean / Var
GammaMarCyB <- Var / Mean
#RanMarine <- GammaSample(GammaMarCyA, GammaMarCyB)
RanMarine <- rgamma(1, GammaMarCyA, scale=GammaMarCyB)
}
}
if(RanMarine < 0) RanMarine <- 0
LastRanMarine <- RanMarine
MS <- RanMarine^inputs$BSRc
if(MS < 0) MS <- 0
#Set the flow parameters
if(inputs$TrndCycF =="Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
}
if(inputs$TrndCycF =="Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF =="Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- inputs$FlowAve
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0) RanFlow <- 0
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC4"){
if(inputs$BSRa * Escpmnt > SRb){
r=SRb * MS * FWS
}else{ r <- inputs$BSRa * Escpmnt * MS * FWS }
}
if(inputs$SRType == "RIC4"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * MS * FWS
}
if(inputs$SRType == "BEV4"){
r <- 0;
if(Escpmnt > 0) r <- MS * FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
} #end param 4 cases
if(inputs$depen == "YES" & Escpmnt < inputs$DL1 + 1){
if(Escpmnt < inputs$DL2 + 1){
fac <- Escpmnt / inputs$DL2 * inputs$DR
}else{
fac <- ((1 - (inputs$DL1 - Escpmnt) / (inputs$DL1 - inputs$DL2)) * (1 - inputs$DR)) + inputs$DR }
if(fac < 0) fac <- 0
AEQRecruits <- fac * AEQRecruits
}
# COMPUTE RECRUITS that are age 1;
# RecruitsAtAge1 is from Recruits() and is "total fraction of age 1 ind that eventually return"
# AEQRecruits/RecruitsAtAge1 <- Age 1 or Cohort[1]
CohortAge1 <- AEQRecruits / staticvars$RecruitsAtAge1
# ADD MARINE SURVIVAL IF STOCK RECRUIT FUNCTION IS SPAWNER TO SMOLT
if(inputs$MarSurv == "YES"){
BetaVariate <- CompBetaVariate(inputs$BetaMarA, inputs$BetaMarB)
CohortAge1 <- CohortAge1 * BetaVariate
}
repvars$Cohort[1] <- CohortAge1
repvars$LastRanFlow <- LastRanFlow
repvars$LastRanError <- LastRanError
repvars$LastRanMarine <- LastRanMarine
#mf addition:
repvars$AEQRecruits <- AEQRecruits
return(repvars)
}#END CompRecruits
#' @title Title
#'
#' @param YearStats
#' @param Year
#' @param inputs
#' @param repvars
#' @param staticvars
#' @param BufSRb
#'
#' @return
#' @keywords internal
#'
#' @examples
CompRecruits2 <- function(YearStats, Year, inputs, repvars, staticvars, BufSRb){
# SRb=inputs$BSRb
#cummulative stats
LastRanFlow <- repvars$LastRanFlow
LastRanError <- repvars$LastRanError
LastRanMarine <- repvars$LastRanMarine
#Horrid; renaming TotEscpmnt to Escpmnt
if(inputs$EscChoice == "YES"){
Escpmnt.sum <- YearStats$TotAdultEscpmnt[Year]
}else{ Escpmnt.sum <- YearStats$TotEscpmnt[Year] }
if(Escpmnt.sum < 1){ Escpmnt.sum <- 0 }
#COMPUTE AEQ RECRUITMENT
AEQRecruits <- eval(parse(text=inputs$SR.FUN))
if(inputs$SRType == "HOC2"){
if(inputs$BSRa * Escpmnt > SRb){
AEQRecruits <- SRb
}else{
AEQRecruits <- inputs$BSRa * Escpmnt
}
}
if(inputs$SRType == "RIC2"){
AEQRecruits <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb)
}
if(inputs$SRType == "BEV2"){
AEQRecruits <- 0
if(Escpmnt > 0) AEQRecruits <- 1/((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
if(any(inputs$SRType==c("HOC2","RIC2","BEV2")) & inputs$SurvScale == "YES"){
#StockScalar <- GammaSample(inputs$SRErrorA, inputs$SRErrorB)
#StockScalar <- rgamma(1, inputs$SRErrorA, scale=inputs$SRErrorB)
StockScalar <- eval(parse(text = inputs$SRError.FUN))
AEQRecruits <- StockScalar * AEQRecruits
}
if(any(inputs$SRType==c("HOC3", "RIC3", "BEV3"))){
if(inputs$TrndCycF == "Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
#RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- eval(parse(text = inputs$FWError.FUN))
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
} #Autoc
if(inputs$TrndCycF == "Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
#EEH: In old code this was GammaFlowA which meant the global was overwritten
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF == "Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0){ RanFlow <- 0 }
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC3")
if(inputs$BSRa * Escpmnt > SRb){ r <- SRb * FWS
}else{ r <- inputs$BSRa * Escpmnt * FWS }
if(inputs$SRType == "RIC3"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * FWS
}
if(inputs$SRType == "BEV3"){
r <- 0;
if(Escpmnt > 0) r <- FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * repvars$LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
}#end 3 parameter cases
if(any(inputs$SRType==c("HOC4", "RIC4", "BEV4"))){
#Set Marine Surv Params
if(inputs$TrndCycM == "Autoc"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
#RanMarine <- GammaSample(inputs$GammaMarA, inputs$GammaMarB)
RanMarine <- eval(parse(text = inputs$MarError.FUN))
RanMarine <- Autocorrel(inputs$TCM1, LastRanMarine, RanMarine)
}
}
if(inputs$TrndCycM == "Trend"){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
RanMarine <- inputs$MarAve
}else{
Mean <- Trend(inputs$TCM1, Year, inputs$MarAve, inputs$TCM2)
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- Mean
}else{
GammaMarTrA <- Mean * Mean / Var
GammaMarTrB <- Var / Mean
#RanMarine <- GammaSample(GammaMarTrA, GammaMarTrB)
RanMarine <- rgamma(1, GammaMarTrA, scale=GammaMarTrB)
}
}
}
if(inputs$TrndCycM == "Cycle"){
Mean <- Cycle(inputs$TCM1, inputs$TCM2, inputs$TCM3, Year)
# should use TCF4 as scalar in place of MarAve
Var <- (inputs$MarCV * Mean)^2
if(Var == 0){
RanMarine <- inputs$MarAve
}else{
GammaMarCyA <- Mean * Mean / Var
GammaMarCyB <- Var / Mean
#RanMarine <- GammaSample(GammaMarCyA, GammaMarCyB)
RanMarine <- rgamma(1, GammaMarCyA, scale=GammaMarCyB)
}
}
if(RanMarine < 0) RanMarine <- 0
LastRanMarine <- RanMarine
MS <- RanMarine^inputs$BSRc
if(MS < 0) MS <- 0
#Set the flow parameters
if(inputs$TrndCycF =="Autoc"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
#RanFlow <- GammaSample(inputs$GammaFlowA, inputs$GammaFlowB)
RanFlow <- rgamma(1, inputs$GammaFlowA, scale=inputs$GammaFlowB)
RanFlow <- Autocorrel(inputs$TCF1, LastRanFlow, RanFlow)
}
}
if(inputs$TrndCycF =="Trend"){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
RanFlow <- inputs$FlowAve
}else{
Mean <- Trend(inputs$TCF1, Year, inputs$FlowAve, inputs$TCF2)
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- Mean
}else{
GammaFlowTrA <- Mean * Mean / Var
GammaFlowTrB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowTrA, GammaFlowTrB)
RanFlow <- rgamma(1, GammaFlowTrA, scale=GammaFlowTrB)
}
}
}
if(inputs$TrndCycF =="Cycle"){
Mean <- Cycle(inputs$TCF1, inputs$TCF2, inputs$TCF3, Year)
# should use TCF4 as scalar in place of FlowAve
Var <- (inputs$FlowCV * Mean)^2
if(Var == 0){
RanFlow <- inputs$FlowAve
}else{
GammaFlowCyA <- Mean * Mean / Var
GammaFlowCyB <- Var / Mean
#RanFlow <- GammaSample(GammaFlowCyA, GammaFlowCyB)
RanFlow <- rgamma(1, GammaFlowCyA, scale=GammaFlowCyB)
}
}
if(RanFlow < 0) RanFlow <- 0
LastRanFlow <- RanFlow
FWS <- exp(inputs$BSRd * RanFlow)
if(FWS < 0) FWS <- 0
if(inputs$SRType == "HOC4"){
if(inputs$BSRa * Escpmnt > SRb){
r=SRb * MS * FWS
}else{ r <- inputs$BSRa * Escpmnt * MS * FWS }
}
if(inputs$SRType == "RIC4"){
r <- inputs$BSRa * Escpmnt * exp(-Escpmnt / SRb) * MS * FWS
}
if(inputs$SRType == "BEV4"){
r <- 0;
if(Escpmnt > 0) r <- MS * FWS / ((1 / SRb) + (1 / inputs$BSRa) * (1 / Escpmnt))
}
AEQRecruits <- r
if(inputs$SurvScale == "YES" & AEQRecruits > 0){
RanError <- inputs$ResCorParameter * LastRanError + (1 - inputs$ResCorParameter) * rnorm(1) * sqrt(inputs$SRErrorB)
LastRanError <- RanError
AEQRecruits <- AEQRecruits * exp(RanError)
}
} #end param 4 cases
if(inputs$depen == "YES" & Escpmnt < inputs$DL1 + 1){
if(Escpmnt < inputs$DL2 + 1){
fac <- Escpmnt / inputs$DL2 * inputs$DR
}else{
fac <- ((1 - (inputs$DL1 - Escpmnt) / (inputs$DL1 - inputs$DL2)) * (1 - inputs$DR)) + inputs$DR }
if(fac < 0) fac <- 0
AEQRecruits <- fac * AEQRecruits
}
# COMPUTE RECRUITS that are age 1;
# RecruitsAtAge1 is from Recruits() and is "total fraction of age 1 ind that eventually return"
# AEQRecruits/RecruitsAtAge1 <- Age 1 or Cohort[1]
CohortAge1 <- AEQRecruits / staticvars$RecruitsAtAge1
# ADD MARINE SURVIVAL IF STOCK RECRUIT FUNCTION IS SPAWNER TO SMOLT
if(inputs$MarSurv == "YES"){
BetaVariate <- CompBetaVariate(inputs$BetaMarA, inputs$BetaMarB)
CohortAge1 <- CohortAge1 * BetaVariate
}
repvars$Cohort[1] <- CohortAge1
repvars$LastRanFlow <- LastRanFlow
repvars$LastRanError <- LastRanError
repvars$LastRanMarine <- LastRanMarine
return(repvars)
}#END CompRecruits2
#' @title Title
#'
#' @param BufNum
#' @param inputs
#' @param BufSRb
#' @param YearStats
#' @param SummaryStats
#'
#' @return
#' @keywords internal
#'
#' @examples
CompStatsEEH <- function(BufNum, inputs, BufSRb, YearStats, SummaryStats){
#BufNum is BufMax, an integer
NYears <- inputs$NYears
NRuns <- inputs$NRuns
SRErrorA <- inputs$SRErrorA
SRErrorA <- inputs$SRErrorB
SRa <- inputs$BSRa
if(inputs$StepFunc=="POP"){
SRb <- BufSRb
}else{
SRB <- inputs$BSRb
}
SRc <- inputs$BSRc
SRd <- inputs$BSRd
SurvScale <- inputs$SurvScale
#In VB code Escapmnt (over ages) is redefined to be TotEscpmnt (sum over ages)
#kept terminology closer to what is being donw
TotEscpmnt <- YearStats$TotEscpmnt #a matrix of Escpmnt for each year
# COMPUTE STATISTICS
# Sum up over all reps and divide by reps
# COMPUTE AVERAGE MORTALITY OVER REPETITIONS
#divide by NRuns because we want the average over all runs
#Note in VB code, AEQMort is initially by age and then in SaveYearData, it is reassigned to be TotAEQMort.
#Here I keep AEQMort and TotAEQMort separate
SummaryStats$AvgAEQMort[BufNum] <- SummaryStats$AvgAEQMort[BufNum]+mean(YearStats$TotAEQMort)/inputs$NRuns
# COMPUTE AVERAGE CALENDAR YEAR HARVEST RATE
#CalendarHR <- YearStats$CalendarHR
#divide by NRuns because we want the average over all runs
SummaryStats$AvgCaHR[BufNum] <- SummaryStats$AvgCaHR[BufNum] + mean(YearStats$CalendarHR)/inputs$NRuns
# COMPUTE PROPORTION OF TIMES ESCAPEMENT WAS LESS THAN LOWER LEVEL
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$AvgECrit[BufNum] <- SummaryStats$AvgECrit[BufNum] + sum(TotEscpmnt < inputs$ECrit)/(inputs$NYears*inputs$NRuns)
# FIND MINIMUM, MAXIMUM, AND AVERAGE ESCAPEMENT
SummaryStats$MinEscpmnt[BufNum,] <- pmin(SummaryStats$MinEscpmnt[BufNum,],TotEscpmnt,na.rm=TRUE)
SummaryStats$MaxEscpmnt[BufNum,] <- pmax(SummaryStats$MaxEscpmnt[BufNum,],TotEscpmnt,na.rm=TRUE)
#divide these by NRuns because we want the average over all runs
SummaryStats$AvgEscpmnt[BufNum, ] <- SummaryStats$AvgEscpmnt[BufNum, ] + TotEscpmnt/inputs$NRuns
# COMPUTE PROPORTION OF TIMES ENDING ESCAPEMENT WENT BELOW QET (DL2) (depensation)
#uses last EndAv years
assessmentYears <- (inputs$NYears - (inputs$EndAv - 1)):inputs$NYears
# divide by NRuns at end to get frac of reps below
if(any(TotEscpmnt[assessmentYears] < (inputs$DL2 + 1))){
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$PropExt[BufNum] <- SummaryStats$PropExt[BufNum] + 1/inputs$NRuns
}
# COMPUTE PROPORTION OF TIMES ENDING ESCAPEMENT EXCEEDED UPPER LEVEL
EscpmntPositive <- TotEscpmnt
EscpmntPositive[EscpmntPositive<1] <- 1
#assessmentYears defined above; uses geometric mean over last EndAv years
GeomeanEscpmnt <- exp(mean(log(EscpmntPositive)[assessmentYears]))
if(GeomeanEscpmnt > inputs$ERecovery){
# add 1/NRuns because we want the fraction out of all runs
SummaryStats$PropRec[BufNum] <- SummaryStats$PropRec[BufNum] + 1/inputs$NRuns
}
# COMPUTE BROOD YEAR HARVEST RATES
BYrAEQMort <- rep(0, inputs$NYears - (inputs$MaxAge-inputs$MinAge)) # ReDim BYrAEQMort(-1 To NYears% - 5)
BYrEscpmnt <- rep(0, inputs$NYears - (inputs$MaxAge-inputs$MinAge)) #ReDim BYrEscpmnt(-1 To NYears% - 5)
# AGE 2
Age <- 2
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,2]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,2]
# AGE 3
Age<- 3
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,3]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,3]
# AGE 4
Age <- 4
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,4]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,4]
# AGE 5
Age <- 5
Years <- 1:(inputs$NYears - (inputs$MaxAge-inputs$MinAge)) + (Age-inputs$MinAge)
Byrs <- Years - Age
BYrAEQMort[Byrs+2] <- BYrAEQMort[Byrs+2] + YearStats$AEQMort[Years,5]
BYrEscpmnt[Byrs+2] <- BYrEscpmnt[Byrs+2] + YearStats$Escpmnt[Years,5]
TempBYrHR <- BYrAEQMort / (BYrAEQMort + BYrEscpmnt)
TempBYrHR[BYrAEQMort < 1E-16] <- 0
#divide these by NRuns because we want the average over all runs
SummaryStats$BufAvgBYrHR[BufNum] <- SummaryStats$BufAvgBYrHR[BufNum] + mean(TempBYrHR)/inputs$NRuns
SummaryStats$AvgBYrHR[BufNum, ] <- SummaryStats$AvgBYrHR[BufNum, ] + TempBYrHR/inputs$NRuns
#don't divide these by NRuns because Max/Min
SummaryStats$MaxBYrHR[BufNum, ] <- pmax(SummaryStats$MaxBYrHR[BufNum, ], TempBYrHR, na.rm=TRUE)
SummaryStats$MinBYrHR[BufNum, ] <- pmax(SummaryStats$MinBYrHR[BufNum, ], TempBYrHR, na.rm=TRUE)
SummaryStats$AveRanMarine[BufNum] <- SummaryStats$AveRanMarine[BufNum] + mean(YearStats$RanMarine)/inputs$NRuns
SummaryStats$AveRanFlow[BufNum] <- SummaryStats$AveRanFlow[BufNum] + mean(YearStats$RanFlow)/inputs$NRuns
return(SummaryStats)
}#END CompStatsEEH
#' @title Title
#'
#' @param TotAdultEscpmnt
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
CompStockStatus <- function(TotAdultEscpmnt, inputs){
# ***** CompStockStatus *****
# Sub CompStockStatus(TotAdultEscpmnt As Double, Status%)
# Dim Break%
# Dim EscpmntCheck$
#
# Break% <- 1
# EscpmntCheck$ <- "True"
# While Break% <= NumBreakPoints% And EscpmntCheck$ <- "True"
# If TotAdultEscpmnt > EscpmntBreakPoint(Break%) Then
# EscpmntCheck$ <- "True"
# Break% <- Break% + 1
# Else
# EscpmntCheck$ <- "False"
# End If
# Wend
#
# Status% <- Break%
#
# End Sub
Break <- 1
EscpmntCheck <- "True"
while(Break <= inputs$NumBreakPoints & EscpmntCheck == "True"){
if(TotAdultEscpmnt > inputs$EscpmntBreakPoint[Break]){
EscpmntCheck <- "True"
Break <- Break + 1
}else{
EscpmntCheck <- "False"
}
}
Status <- Break
return(Status)
}#END CompStockStatus
#' @title Title
#'
#' @param a
#' @param p
#' @param s
#' @param y
#'
#' @return
#' @keywords internal
#'
#' @examples
Cycle <- function(a, p, s, y){
# ***** FUNCTION Cycle *****
# Compute cyclic variable, a <- amplitude, p <- period, s <- start, y <- year, x <- mean value of variable
# NJS: created 7/9/02 corrected 9/16/03
#
# Function Cycle(a As Double, p As Double, s As Double, y, x As Double) As Double
# a is amplitude, p is period, s is starting point, y time period
# what is x doing here? It is average value and is not needed here.
# Dim cy As Double
# cy <- Sin(2# * 3.141592654 * (y + s - 1) / p)
# 'in good survival, cycle ranges from 1 to a (amplitude)
# in bad survival, cycle ranges from 1/a to 1 (this might be lower than expected)
# If cy >= 0 Then
# cy <- (cy * (a - 1)) + 1
# Else
# cy <- (cy * (1 - (1 / a))) + 1
# End If
# Cycle <- cy + x ' use if x is changed to scalar
# Cycle <- cy
#
# End Function
# a is amplitude, p is period, s is starting point, y time period
# EEH: changed to allow vectors of y
cy <- sin(2 * pi * (y + s - 1) / p)
# in good survival, cycle ranges from 1 to a (amplitude)
# in bad survival, cycle ranges from 1/a to 1 (this might be lower than expected)
newcy <- cy
newcy[cy >= 0]<- (cy[cy >= 0] * (a - 1)) + 1
newcy[cy < 0] <- (cy[cy < 0] * (1 - (1 / a))) + 1
return(newcy)
}#END Cycle
#' @title Title
#'
#' @param Alpha
#' @param Beta
#'
#' @return
#' @keywords internal
#'
#' @examples
GammaSample <- function(Alpha, Beta){
# ***** GammmaSample *****
# Function generates a random gamma deviate with shape
# paramater alpha and scale parameter beta
# ******************************************************************
# if(get("CONSTRUN", pkg_globals)){ #no random variables
# val <- Alpha * Beta # expected value
# }else{
val <- rgamma(1, Alpha, scale=Beta)
# }
return(val)
}#END GammaSample
#' @title Title
#'
#' @param InFile
#'
#' @return
#' @keywords internal
#'
#' @examples
GetInput <- function(InFile){
# ***** GetInput *****
# Read in a .rav file and assign all the variables
# Returns the list of all inputs
#
#The rav file has , as end of input/separator
inputs <- list()
inputs$InFile <- InFile
readit <- function(skip, n){ read.table(InFile,nrows=1,sep=",",stringsAsFactors=FALSE,skip=skip)[1,n] }
# GET TITLE FOR RUN
inputs$Title <- readit(0,1) #line 1
# --------------- RUN PARAMETERS SECTION
# INPUT RANDOM NUMBER SEED, NUMBER OF CYCLES AND REPETITIONS,
# MINIMUM AND MAXIMUM AGE
inputs$RanSeed <- readit(1,1) #line 2, RanSeed
inputs$NRuns <- readit(2,1) #line 2, NRuns
inputs$NYears <- readit(3,1) #line 3, NYears
inputs$MinAge <- readit(4,1) #line 4, NYears
inputs$MaxAge <- readit(4,2) #line 4, NYears
inputs$ConvergeCrit <- readit(5,1) #line 5, ConvergeCrit
inputs$Debugg <- readit(6,1) #line 6, Debugg
inputs$Debugg <- toupper(inputs$Debugg)
if(!(inputs$Debugg %in% c("YES", "NO"))) stop("Unknown debug selection (yes/no only)")
# ----- END OF RUN PARAMETERS SECTION -
# -- STOCK-RECRUIT SECTION ---
# 'GET FORM OF SPAWNER RECRUIT FUNCTION AND PARAMETERS
# njs updated these explanations
# BSRa <- productivity parameter
# BSRb <- density dependent prarameter
# BSRc <- marine survival paramater - M^c
# BSRd <- freshwater survival parameter - exp(dF)(d should be entered as negative)
# HOC2 - Hockey stick R=Min(aS,b) a <- producitvity b=MaxRecruits
# HOC3 - R <- Min(aS,) exp(dF) (freshwater index - may be used to predict smolts)
# HOC4 - R<- Min(aS,) M^c exp(dF)
# RIC2 - Ricker R=aS exp(-bS) a <- productvity b=1/capacity
# RIC3 - R<- aS exp(-bS+dF) (freshwater index - may be used to predict smolts)
# RIC4 - Ricker with marine survival and freshwater survival
# R=aS M^c exp(-bS+dF)
# BEV2 - Beverton-Holt R=1/(b+a/S) a <- 1/productivity b=1/MaxRecruits
# BEV3 - R=1/(b+a/S) exp(dF)(freshwater index - may be used to predict smolts)
# BEV4 - BH with marine survival and freshwater survival
# R=1/(a+b/S)M^c exp(dF)
inputs$SRType <- readit(7,1) #line 8, SRType
inputs$SRType <- toupper(inputs$SRType)
SRType=inputs$SRType
allowedSRType <- c("HOC2", "HOC3", "HOC4", "BEV2", "BEV3", "BEV4", "RIC2", "RIC3", "RIC4")
if(!(SRType %in% allowedSRType)) stop("Unknown stock recruit type")
if(SRType %in% c("HOC2", "RIC2", "BEV2")){
inputs$BSRa <- readit(8,1) #line 4, NYears
inputs$BSRb <- readit(8,2) #line 4, NYears
inputs$BSRc <- 0;
inputs$BSRd <- 0;
inputs$AveEnv <- 1
}
#Skip the next 6 lines and jump to line 16
if(SRType %in% c("HOC3", "RIC3", "BEV3")){
inputs$BSRa <- readit(8,1)
inputs$BSRb <- readit(8,2)
inputs$BSRc <- 0;
inputs$BSRd <- readit(8,3)
#skip next 3 lines
inputs$FlowAve <- readit(12,1)
inputs$FlowCV <- readit(12,2)
inputs$TrndCycF <- readit(13,1)
if(!(inputs$TrndCycF %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for Flow")
inputs$TCF1 <- readit(14,1)
inputs$TCF2 <- readit(14,2)
inputs$TCF3 <- readit(14,3)
inputs$BSRc <- 0
inputs$FlowSD <- inputs$FlowCV * inputs$FlowAve
if(inputs$FlowSD > 0){
inputs$GammaFlowA <- inputs$FlowAve * inputs$FlowAve / (inputs$FlowSD * inputs$FlowSD)
inputs$GammaFlowB <- (inputs$FlowSD * inputs$FlowSD) / inputs$FlowAve
}else{
inputs$GammaFlowA <- 0
inputs$GammaFlowB <- 0
}
inputs$AveEnv <- exp(inputs$BSRd * inputs$FlowAve)
}
if(SRType %in% c("HOC4", "RIC4", "BEV4")){
inputs$BSRa <- readit(8,1)
inputs$BSRb <- readit(8,2)
inputs$BSRc <- readit(8,3);
inputs$BSRd <- readit(8,4)
inputs$MarAve <- readit(9,1)
inputs$MarCV <- readit(9,2)
inputs$TrndCycM <- readit(10,1)
if(!(inputs$TrndCycM %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for marine survival")
inputs$TCM1 <- readit(11,1)
inputs$TCM2 <- readit(11,2)
inputs$TCM3 <- readit(11,3)
inputs$FlowAve <- readit(12,1)
inputs$FlowCV <- readit(12,2)
inputs$TrndCycF <- readit(13,1)
if(!(inputs$TrndCycF %in% c("Autoc", "Trend", "Cycle"))) stop("Unknown trend/cycle function for Flow")
inputs$TCF1 <- readit(14,1)
inputs$TCF2 <- readit(14,2)
inputs$TCF3 <- readit(14,3)
inputs$MarSD <- inputs$MarCV * inputs$MarAve
if(inputs$MarSD > 0){
inputs$GammaMarA <- inputs$MarAve * inputs$MarAve / (inputs$MarSD * inputs$MarSD)
inputs$GammaMarB <- (inputs$MarSD * inputs$MarSD) / inputs$MarAve
}else{
inputs$GammaMarA <- 0
inputs$GammaMarB <- 0
}
inputs$FlowSD <- inputs$FlowCV * inputs$FlowAve
if(inputs$FlowSD > 0){
inputs$GammaFlowA <- inputs$FlowAve * inputs$FlowAve / (inputs$FlowSD * inputs$FlowSD)
inputs$GammaFlowB <- (inputs$FlowSD * inputs$FlowSD) / inputs$FlowAve
}else{
inputs$GammaFlowA <- 0
inputs$GammaFlowB <- 0
}
inputs$AveEnv <- inputs$MarAve^inputs$BSRc * exp(inputs$BSRd * inputs$FlowAve)
}
# depensation with DL1=depensation esc; DL2 <- QET; DR % of predicted R realized at QET
inputs$depen <- readit(15,1)
inputs$depen <- toupper(inputs$depen)
if(!(inputs$depen %in% c("YES", "NO"))) stop("Unknown depensation selection (yes/no only)")
inputs$DL1 <- readit(16,1)
inputs$DL2 <- readit(16,2)
inputs$DR <- readit(16,3)
inputs$EscChoice <- readit(17,1)
inputs$EscChoice <- toupper(inputs$EscChoice)
if(!(inputs$EscChoice %in% c("YES", "NO"))) stop("Unknown escapement choice selection (yes/no only)")
# GET PARAMETERS TO ADD VARIABILITY TO STOCK RECRUIT FUNCTION
# SRErrorA and B are used for Hoc2, Ric2 and Bev2, SRErrorA, SRErrorB and ResCorPar (optional)
# are used for Hoc3, Ric3, Bev3, Hoc4, Ric4, and Bev4
inputs$SurvScale <- readit(18,1)
inputs$SurvScale <- toupper(inputs$SurvScale)
if(!(inputs$SurvScale %in% c("YES", "NO"))) stop("Unknown SR variability selection (yes/no only)")
if(inputs$SurvScale == "YES"){
inputs$SRErrorA <- readit(19,1)
inputs$SRErrorB <- readit(19,2)
inputs$ResCorParameter <- readit(19,3)
}
# GET PARAMETERS FOR SMOLT TO ADULT (MARINE) SURVIVAL IF STOCK RECRUIT FUNCTION
# IS FOR SMOLTS FROM SPAWNERS (FRESHWATER PRODUCTION)
inputs$MarSurv <- readit(20,1)
inputs$MarSurv <- toupper(inputs$MarSurv)
if(!(inputs$MarSurv %in% c("YES", "NO"))) stop("Unknown marine survival selection (yes/no only)")
if(inputs$MarSurv == "YES"){
inputs$BetaMarA <- readit(21,1)
inputs$BetaMarB <- readit(21,2)
inputs$CorrMS <- readit(21,3)
}
# - END OF STOCK RECRUIT SECTION ---
# FISHERY MANAGEMENT PARAMETERS --
# INPUT THE NUMBER OF BREAKPOINTS AND DIMENSION ARRAYS
inputs$NumBreakPoints <- readit(22,1)
# IDENTIFY WHICH LEVEL IS THE BASE REGIME
inputs$BaseRegime <- readit(23,1)
inputs$EscpmntBreakPoint <- rep(NA,inputs$NumBreakPoints+1)
inputs$TargetU <- rep(NA,inputs$NumBreakPoints+1)
#begin dynamically keeping track of line number
cline=24
# INPUT BREAKPOINTS AND TARGET EXPLOITATION RATES
for(BreakPoint in 1:(inputs$NumBreakPoints + 1)){
if(BreakPoint <= inputs$NumBreakPoints){
inputs$EscpmntBreakPoint[BreakPoint] <- readit(cline,1)
inputs$TargetU[BreakPoint] <- readit(cline,2)
cline=cline+1
}else{
inputs$TargetU[BreakPoint] <- readit(cline,1)
cline=cline+1
}
}
# INPUT PARAMETERS FOR MANAGEMENT ERROR
inputs$MgmtError <- readit(cline,1); cline=cline+1
inputs$MgmtError <- toupper(inputs$MgmtError)
if(!(inputs$MgmtError %in% c("YES", "NO"))) stop("Unknown manag error selection (yes/no only)")
# Read in dummy values if management error is No
if(inputs$MgmtError == "YES"){
inputs$GammaMgmtA <- readit(cline,1);
inputs$GammaMgmtB <- readit(cline,2); cline=cline+1
}else{
inputs$GammaMgmtA <- 0
inputs$GammaMgmtB <- 0
#still need to iterate by a step so the import isn't out of sync.
cline=cline+1
}
# ------- END OF FISHERY MANAGMENT INPUTS
inputs$ECrit <- readit(cline,1); cline=cline+1
inputs$ERecovery <- readit(cline,1);
inputs$EndAv <- readit(cline,2); cline=cline+1
# INPUT STEP SIZE AND RANGE FOR TARGET EXPLOITATION RATES OR STARTING ESCAPEMENT
# This program outputs information for a range of either exploitation rates or
# starting escapement levels.
# The BUFFER levels input below represent percentages of the base target
# exploitation rate or the starting escapement level.
# (e.g., a Buffer of .05 means that the ER will be 5% of the base target rate.)
# You may enter a number greater than 1.0, which means that the ER will be
# greater than the base target rate. Under normal runs with 0 breakpoints,
# the base target rate is not used other than to determine start, end, and step
# levels.
inputs$StepFunc <- readit(cline,1); cline=cline+1
inputs$StepFunc <- toupper(inputs$StepFunc)
if(!(inputs$StepFunc %in% c("POP", "ER"))) stop("Unknown step selection")
inputs$BufferStep <- readit(cline,1); cline=cline+1
inputs$BufferStart <- readit(cline,1);
inputs$BufferEnd <- readit(cline,2); cline=cline+1
#integer of the number of steps of ER or Pop Capacity to do # is 1:BufMax
inputs$BufMax <- round((inputs$BufferEnd - inputs$BufferStart) / inputs$BufferStep + 1)
# INPUT LOWER AND UPPER ESCAPEMENT TEST LEVELS
# The LOWER ESCAPEMENT LEVEL (ECrit) is the escapement level used by the
# program to test how often the observed escapements fall below this level.
# It may represent a "critical" level below which the spawner-recruit function
# destabilizes, and the stock increases risk of extinction, or it could be any
# level one just wanted to monitor frequency of achieving less than or equal
# to that level.
# The UPPER ECSAPEMENT LEVEL (ERecovery) is the escapement level used to compare
# against the geometric mean of the last n (EndAv) years of the run. It may be a
# management escapement goal, an interim recovery level, or some other target level.
# ------ BASE STOCK DATA SECTION
#
# INPUT INITIAL POPULATION SIZE BY AGE
# The initial population size by age, for all ages, is used to seed the model.
# In year 1, the management actions are applied to this population, and a portion of
# each size class escapes. By running the program over a range of starting population
# sizes, one can determine what minimum population size (or escapement) is needed
# to guarentee a probability of population viability.
# NJS Although the input allows for different MaxAge% from 5, the internal workings of the
# are fixed to 5 in many places.
inputs$CohortStart <- rep(0,inputs$MaxAge)
for(Age in (inputs$MinAge - 1):inputs$MaxAge){
inputs$CohortStart[Age] <- readit(cline,1); cline=cline+1
}
# INPUT NATURAL MORTALITY RATES
inputs$NatMort <- rep(0,inputs$MaxAge)
for(Age in (inputs$MinAge - 1):inputs$MaxAge){
inputs$NatMort[Age] <- readit(cline,1); cline=cline+1
}
# INPUT MATURATION RATES BY AGE (AEQ will be calculated)
inputs$MatRate <- rep(0,inputs$MaxAge)
for(Age in inputs$MinAge:inputs$MaxAge){
inputs$MatRate[Age] <- readit(cline,1); cline=cline+1
}
# INPUT PRETERMINAL AND MATURE EXPLOITATION RATES BY AGE
# These will be used to proportion target ER by age and fishery
# by def, PTU and MatU are 0 before MinAge
inputs$PTU <- rep(0,inputs$MaxAge)
inputs$MatU <- rep(0,inputs$MaxAge)
for(Age in inputs$MinAge:inputs$MaxAge){
inputs$PTU[Age] <- readit(cline,1);
inputs$MatU[Age] <- readit(cline,2); cline=cline+1
}
return(inputs)
}#END GetInput
#' @title Title
#'
#' @param InFile
#' @param OutFileBase
#' @param NRuns
#' @param silent
#' @param ...
#'
#' @return
#' @keywords internal
#'
#' @examples
Main <- function(InFile=NULL, OutFileBase=NULL, NRuns=-1, silent=FALSE,...){
#require(stringr)
#if not called with input file, then user is prompted to input one
if(is.null(InFile)) InFile <- file.choose()
if(!file.exists(InFile)) stop("Specified input file does not exist.")
if(is.null(OutFileBase)){
OutFileBase <- tools::file_path_sans_ext(basename(InFile))
#excluding this should allow the exclustion of stringr package:
# tmp <- strsplit(InFile,"\\\\")[[1]]
# InFileBase=tmp[length(tmp)]
# tmp <- strsplit(InFileBase,"/")[[1]]
# InFileBase <- tmp[length(tmp)]
# if(str_detect(InFileBase,"[.]")){
# OutFileBase <- strsplit(InFileBase,"[.]")[[1]][1];
# }else{
# OutFileBase <- InFileBase;
# }
}
#Two lists will be passed in and out of functions
# inputs <- list() #is everything from the .rav file
# staticvars <- list() #is anything computed from that; static
# READ INPUT DATA AND CALCULATE AEQs
# direct from .rav file or simple calculation from rav file inputs
#inputs <- GetInput(InFile)
time.start <- Sys.time()
cat("\n", InFile, time.start, "\n")
inputs.list <- read_rav2(InFile)
#input.vars has just the variable values and none of the metadata
inputs <- inputs.list$inputs.vars
if(NRuns>0) inputs$NRuns <- NRuns;
#add the output file names to the inputs
inputs <- SetOutFileNames(OutFileBase, inputs)
out <- RunSims(inputs, silent,...)
# SAVE SUMMARY RESULTS .sum
if(!silent) cat("\nSaving summary...\n")
data.summary <- SaveSummary(out$inputs, out$SummaryStats, out$staticvars)
# SAVE ESCAPEMENT DATA .esc
if(!silent) cat("Saving escapement data...\n")
data.escape.avg <- SaveEscpmntData(out$inputs, out$SummaryStats)
# SAVE BROOD YEAR EXPLOITATION RATE DATA .byr
if(!silent) cat("Saving BYr year data...\n")
data.byer.avg <- SaveBYrData(out$inputs, out$SummaryStats)
vrap.output <- list(list(inputs=inputs, output=out, data.summary=data.summary, data.escape.avg=data.escape.avg, data.byer.avg=data.byer.avg ))
names(vrap.output) <- tools::file_path_sans_ext(InFile)
saveRDS(vrap.output,file = paste0(tools::file_path_sans_ext(InFile), ".rds"))
# return(list(inputs=inputs, output=out, data.summary=data.summary, data.escape.avg=data.escape.avg, data.byer.avg=data.byer.avg ))
cat(Sys.time()-time.start,"\n")
}#END Main
#' @title Title
#'
#' @param prop
#' @param prev
#'
#' @return
#' @keywords internal
#'
#' @examples
progressBar <- function(prop=0, prev=0){
if (prev < 50) {
if (prop > 1) {
prop <- prop/100
}
bars <- round(prop * 50, digits <- 0) - prev
prev <- bars + prev
if (prop == 0) {
cat(" |2% |20% |40% |60% |80% |100%\n",
"Progress: ", sep <- "")
}
else {
while (bars > 0) {
cat("|")
bars <- bars - 1
}
}
if (prev == 50) {
cat("\n")
}
return(prev)
}
prev
}#END progressBar
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
Recruits <- function(inputs){
# ***** Recruits *****
# Function Recruits() As Double
# Compute factor to convert calculated spawner equivalent
# production to age cohort (source is PSC Chinook Model).
# #
# Tmp <- 0
# X9 <- 1 - NatMort(1)
# For Age% <- MinAge% To MaxAge%
# X9 <- X9 * (1 - NatMort(Age%))
# Tmp <- Tmp + X9 * MatRate(Age%)
# X9 <- X9 * (1 - MatRate(Age%))
# Next Age%
# Recruits <- Tmp
#
# End Function
# The SR function gets us the AEQRecruits from spawners in year t.
# We needs to then translate that to age 1 indiv in pop (Cohort[1])
# We know AEQRecruit. How many Age 1 individuals does that translate to? Age1 * (1- total fraction lost) <- AEQRecruits
# So Age1 <- AEQRecruits/(1-total fraction lost)
# Tmp here is total fraction of age 1 ind that eventually return
# AEQRecruits/Tmp <- Age 1 or Cohort[1]
Tmp <- 0
X9 <- 1 - inputs$NatMort[1] #survival ag 1
for(Age in inputs$MinAge:inputs$MaxAge){
X9 <- X9 * (1 - inputs$NatMort[Age])
Tmp <- Tmp + X9 * inputs$MatRate[Age]
X9 <- X9 * (1 - inputs$MatRate[Age])
}
#AEQRecruits / compvars$RecruitsAtAge1
return(Tmp) #Recruits at age 1
}#END Recruits
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
RepInit <- function(inputs){
rtn.list <- list()
#Set default values used when RanFlow and RanMarine are used (2 and 3 param SR models)
rtn.list$LastRanFlow <- 0
rtn.list$LastRanMarine <- 0
rtn.list$Cohort <- inputs$CohortStart
# SET INITIAL SEED FOR AUTOCORRELATED RESIDUALS
if(inputs$SurvScale == "YES"){
rtn.list$LastRanError <- rnorm(1, 0, sd=sqrt(inputs$SRErrorB))
}
# SET INITIAL VALUES FOR LAST RANMARINE AND RANFLOW
# EEH: VB code was misisng HOC3 and BEV3 from this
if(any(inputs$SRType==c("HOC3", "BEV3", "RIC3"))){
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
rtn.list$LastRanFlow <- inputs$FlowAve
}else{
if(inputs$TCF2 > 0){ rtn.list$LastRanFlow <- inputs$TCF2;
}else{ rtn.list$LastRanFlow <- rgamma(1,inputs$GammaFlowA, scale=inputs$GammaFlowB); }
}
}
if(any(inputs$SRType==c("HOC4", "BEV4", "RIC4"))){
if(inputs$GammaMarA + inputs$GammaMarB == 0){
rtn.list$LastRanMarine <- inputs$MarAve
}else{
if(inputs$TCM2 > 0){ rtn.list$LastRanMarine <- inputs$TCM2
}else{ rtn.list$LastRanMarine <- rgamma(1,inputs$GammaMarA, scale=inputs$GammaMarB) }
}
if(inputs$GammaFlowA + inputs$GammaFlowB == 0){
rtn.list$LastRanFlow <- inputs$FlowAve
}else{
if(inputs$TCF2 > 0){ rtn.list$LastRanFlow <- inputs$TCF2;
}else{ rtn.list$LastRanFlow <- rgamma(1,inputs$GammaFlowA, scale=inputs$GammaFlowB); }
}
}
return(rtn.list)
}#END RepInit
#' @title Title
#'
#' @param inputs
#' @param silent
#'
#' @return
#' @keywords internal
#'
#' @examples
RunSims <- function(inputs, silent){
# RunSims takes the input list and runs the VRAP functions
#Set up list that will hold static computed variables. These don't change with each rep or buffer
staticvars <- list()
# CALCULATE AEQs
# staticvars are computed variables; they are static
staticvars$AEQ <- AEQcalc(inputs)
# COMPUTE FACTOR TO TRANSLATE AEQ RECRUITMENT TO AGE 1
staticvars$RecruitsAtAge1 <- Recruits(inputs)
ptm <- proc.time()
# PROGRAM EXECUTION SECTION
Buffer <- inputs$BufferStart
#Set up the SummaryStats
SummaryStats <- SetupSummaryStats(inputs)
#Set up the YearStats
YearStats <- list()
YearStats$AEQMort <- matrix(0, inputs$NYears, inputs$MaxAge)
YearStats$Escpmnt <- matrix(0, inputs$NYears, inputs$MaxAge)
YearStats$TotAEQMort <- matrix(NA, inputs$NYears)
YearStats$TotAdultEscpmnt <- matrix(NA, inputs$NYears)
YearStats$TotEscpmnt <- matrix(NA, inputs$NYears)
YearStats$RanFlow <- matrix(NA, inputs$NYears)
YearStats$RanMarine <- matrix(NA, inputs$NYears)
#MF addition:
#inputs$output.raw is a string of the variable names to be saved for output
#if some variables are named in output.raw then build rawoutput to store them:
if(exists("output.raw", where=inputs)){
if(!is.na(inputs$output.raw)){ rawoutput <- vector(mode = "list", length = inputs$BufMax*inputs$NRuns)}
}
YearStats$Recruits <- rep(NA, inputs$NYears)
YearStats$AEQRecruits <- rep(NA, inputs$NYears)
YearStats$TargetUError <- rep(NA, inputs$NYears)
YearStats$ptu.scaled <- matrix(NA, inputs$NYears, 5)
YearStats$MatU.scaled <- matrix(NA, inputs$NYears, 5)
YearStats$Cohort <- matrix(NA, inputs$NYears, 5)
YearStats$tmpdat <- rep(NA, inputs$NYears)
if(!silent) prev <- progressBar()
# For each ER or Pop Cap level, go loop through NRuns, and for each NRun, loop through Year
for(BufNum in 1:inputs$BufMax){
# INITIALIZE BUFFER SPECIFIC PARAMETERS AND ARRAYS
out <- BufferInit(Buffer, inputs)
inputs$DR <- out$DR
#This is the ER or Pop Cap to use for this loop
#used in CompEscpmnt
BufTargetU <- out$BufTargetU
#This is the SRb to use when StepFunc <- "POP"; not used when StepFunc <- "ER"
#used in CompRecruits and CompStats
BufSRb <- out$BufSRb
# REPETITION LOOP
for(Rep in 1:inputs$NRuns){
if(!silent) prev <- progressBar((inputs$NRuns*(BufNum-1)+Rep)/(inputs$BufMax*inputs$NRuns),prev)
# INITIALIZE REPETITION SPECIFIC PARAMETERS AND ARRAYS
# repvar is a list of variables that change each year: Cohort, LastRanError, LastRanFlow, LastRanMarine
# These are updated with each year
# starts Cohort at the init # at each age in the input file
# Cohort seems to mean population #s at each age
# YearStats is a list of variables that are saved for each year
# SummaryStats passed in to fix first year Ran to be same across reps
repvars <- RepInit(inputs)
if(Rep == 1){
SummaryStats$FirstRanMarine <- repvars$LastRanMarine
SummaryStats$FirstRanFlow <- repvars$LastRanFlow
}
# BEGIN Year LOOP
for(Year in 1:inputs$NYears){
#Save this for Summary at end
YearStats$RanMarine[Year] <- repvars$LastRanMarine
YearStats$RanFlow[Year] <- repvars$LastRanFlow
# NATURAL MORTALITY PROCESS
# Update Cohort (pop size) at Age with Natural Mortality
repvars$Cohort <- CompNatMort(inputs, repvars$Cohort)
# COMPUTE ESCAPEMENT USING BASE LEVEL EXPLOITATION RATE
#Escapment is spawners
# cat("\n",Rep, class(staticvars$AEQ), Year, BufTargetU)
YearStats <- CompEscpmnt(inputs$BaseRegime, Year, inputs, BufTargetU, repvars$Cohort, staticvars$AEQ, YearStats)
YearStats$TempCohort <- round(YearStats$TempCohort,0)
repvars$TempCohort <- YearStats$TempCohort
#commented out since I want to test without the regime feature first
# # CHECK STOCK STATUS
# # Going to be a number between 1 and NumBreaks
# Status <- CompStockStatus(YearStats$TotAdultEscpmnt, inputs)
#
# # ADJUST REGIME IF WARRANTED
# if(Status !<- inputs$BaseRegime){
# YearStats=CompEscpmnt(Status, Year, inputs, BufTargetU, repvars$Cohort, staticvars$AEQ, YearStats)
# }
# AGE COHORT
# now update cohort (non-spawner pop size at age) for next year using the updated "cohort" from CompEscpmnt
repvars$Cohort <- CompAgeCohort(repvars$TempCohort, repvars$Cohort, inputs)
# COMPUTE RECRUITS FROM ESCAPEMENT
# Uses SR (Escpmnt to AEQRecruits) function with error added to that to get AEQRecruits from this years spawners (Escpmnt)
# Needs to then translate that to age 1 indiv in pop (Cohort[1])
# We know AEQRecruit. How many Age 1 individuals does that translate to? Age1 * (1- total fraction lost) <- AEQRecruits
# So Age1 <- AEQRecruits/(1-total fraction lost)
# Set Cohort[1], LastRanError, LastRanFlow, LastRanMarine
repvars <- CompRecruits(YearStats, Year, inputs, repvars, staticvars, BufSRb)
#repvars <- CompRecruits2(YearStats, Year, inputs, repvars, staticvars, BufSRb)
# SAVE YEAR DATA
# stores the year data that will be needed for the statistics at the end
YearStats <- SaveYearData(Year, YearStats)
#mf: I don't like this use of the term recruits. This is summing, the returns by age to get total returns for a return year:
YearStats$Recruits[Year] <- sum(repvars$Cohort, na.rm = TRUE)
YearStats$TargetU <- BufTargetU
#mf addition:
YearStats$AEQRecruits[Year] <- repvars$AEQRecruits
} #for loop for Year
#this cumulates some stats over each rep
SummaryStats<- CompStatsEEH(BufNum, inputs, BufSRb, YearStats, SummaryStats)
#mf:here is the output of all data:
if(exists("rawoutput")){
output.varnames <- unlist(strsplit(inputs$output.raw, split = " +"))
#convert appropriate data types to integer:
integer.data <- c("AEQMort", "Escpmnt", "TotAEQMort", "TotAdultEscpmnt", "TotEscpmnt", "Recruits", "AEQRecruits", "TempCohort")
for(i in names(YearStats)){
if(i %in% integer.data) storage.mode(YearStats[[i]]) <- "integer"}
YearStats$rep <- Rep
rawoutput[[(BufNum-1) *inputs$NRuns + Rep]] <- YearStats[output.varnames]
}#end if(exists(rawoutput)){
} #for loop for Rep
# next ER level or Pop Cap (SRb) level
Buffer <- Buffer + inputs$BufferStep
} #for loop for BufNum
#if not saving any detailed output
if(!exists("rawoutput"))rawoutput <- NULL
comp.time <- proc.time() - ptm
return(list(inputs=inputs, SummaryStats=SummaryStats, staticvars=staticvars, time=comp.time, rawoutput=rawoutput))
}#END RunSims
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveBYrData <- function(inputs, SummaryStats, writefile=FALSE){
file <- inputs$OutFileByr
output <- ""
BufNum <- 0
for(Buffer in seq(inputs$BufferStart,inputs$BufferEnd,inputs$BufferStep)){
BufNum <- BufNum + 1
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
for(Byr in -1:(inputs$NYears - inputs$MaxAge)){
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff,digits=2), nsmall=2,width=6));
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(Byr,nsmall=0,width=7));
#plus 2 because R indexing starts at 1
output <- c(output, format(round(SummaryStats$MinBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, format(round(SummaryStats$AvgBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, format(round(SummaryStats$MaxBYrHR[BufNum, Byr+2],digits=3), nsmall=3,width=8));
output <- c(output, "\r\n");
} #for loop for Byr
} #for loop for Buffer
if(writefile) cat(file=file, output)
invisible(output)
}#END SaveBYrData
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveEscpmntData <- function(inputs, SummaryStats, writefile=FALSE){
file <- inputs$OutFileEsc;
output <- ""
BufNum <- 0
for(Buffer in seq(inputs$BufferStart,inputs$BufferEnd,inputs$BufferStep)){
BufNum <- BufNum + 1
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
for(Year in 1:inputs$NYears){
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff,digits=2), nsmall=2,width=6));
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(Year,nsmall=0,width=7));
#plus 2 because R indexing starts at 1
output <- c(output, format(round(SummaryStats$MinEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, format(round(SummaryStats$MaxEscpmnt[BufNum, Year],digits=0), nsmall=0,width=8));
output <- c(output, "\r\n");
} #for loop for year
} #for loop for buffer
if(writefile) cat(file=file, output)
invisible(output)
}#END SaveEscpmntData
#' @title Title
#'
#' @param inputs
#' @param SummaryStats
#' @param staticvars
#' @param writefile
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveSummary <- function(inputs, SummaryStats, staticvars, writefile=FALSE){
file <- inputs$OutFileSum;
output=""
# PRINT HEADER INFORMATION
output <- c(output, format(paste("RAPVIABILITY (R) Version "),width=54), "Date:", format(Sys.Date()),"\r\n");
output <- c(output, format("Title:",width=14), inputs$Title, "\r\n");
output <- c(output, format("Input File:",width=19), inputs$InFile, "\r\n");
output <- c(output, "Copy of Input File:", inputs$InFile, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Basic Simulation Input Parameters:","\r\n");
output <- c(output, " ","# of Years=", inputs$NYears, " # of Reps=", inputs$NRuns, " HR Conv.Crit=", inputs$ConvergeCrit, " Seed", inputs$RanSeed, "\r\n");
output <- c(output, " Range start ", inputs$BufferStart, " end ", inputs$BufferEnd, " by ", inputs$BufferStep, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Stock Recruit Function Input Parameters:", "\r\n");
output <- c(output, " ", "Function Type: ", inputs$SRType, "\r\n");
if(inputs$SRType == "HOC2"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=maxrecruits=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "HOC3"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)* exp(d*FWindex)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, " b=maxrecruits=", inputs$BSRb, sep=""),width=34), " d=FW parameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "HOC4"){
output <- c(output, " ", "Recruits <- min(a*Spawners, b)* exp(d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits=", inputs$BSRb, " c=MSparameter=", inputs$BSRc, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "RIC2"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b)", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=capacity=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "RIC3"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b + d*FWindex)", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=capacity=", inputs$BSRb, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "RIC4"){
output <- c(output, " ", "Recruits <- a * Spawners * exp(-Spawners/b + d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=capacity=", inputs$BSRb, " c=MSparameter=", inputs$BSRc, " d=FW parameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "BEV2"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] ", "\r\n");
output <- c(output, format(paste(" ", "a=productivity=", inputs$BSRa, sep=""),width=24), " b=maxrecruits=", inputs$BSRb, "\r\n");
}
if(inputs$SRType == "BEV3"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] * exp(d*FWindex)", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits", inputs$BSRb, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
if(inputs$SRType == "BEV3"){
output <- c(output, " ", "Recruits <- 1/[(1/b) + 1/(a*Spawners)] * exp(d*FWindex) * MarineIndex^c", "\r\n");
output <- c(output, " ", "a=productivity=", inputs$BSRa, " b=maxrecruits", inputs$BSRb, " c=MSparameter", inputs$BSRc, " d=FWparameter<- ", inputs$BSRd, "\r\n");
}
output <- c(output, "\r\n");
if(inputs$SurvScale == "YES"){
if(inputs$SRType %in% c("HOC2", "RIC2", "BEV2")){
output <- c(output, "Stock-Recruit Error Parameters (gamma distr.) [R=f(s)*e]:", "\r\n");
output <- c(output, format(paste(" ", "A=", inputs$SRErrorA, sep=""), width=29), " B=", inputs$SRErrorB, "\r\n");
output <- c(output, "\r\n");
}
if(inputs$SRType %in% c("HOC3", "RIC3", "BEV3", "HOC4", "RIC4", "BEV4")){
output <- c(output, "Stock-Recruit Error Parameters [R=f(S)*exp(e)]:", "\r\n");
output <- c(output, format(paste(" ", "MSE=", inputs$SRErrorB, sep=""), width=29), " ResCor=", inputs$ResCorParameter, "\r\n");
output <- c(output, "\r\n");
}
}
if(inputs$depen == "YES"){
output <- c(output, "Depensation at escap:", inputs$DL1, " QET:", inputs$DL2, "fraction of depensation at QET", inputs$DR, "\r\n");
output <- c(output, "\r\n");
}
if(inputs$MarSurv == "YES"){
if( inputs$SRType %in% c("HOC2", "RIC2", "BEV2", "HOC3", "RIC3", "BEV3")){
output <- c(output, "Smolt to Adult Survival Rate Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Beta_A=", inputs$BetaMarA, sep=""), width=34), " Beta_B=", inputs$BetaMarB, "\r\n");
output <- c(output, "\r\n");
}else{
output <- c(output, "Should not be using this for ", inputs$SRType, "\r\n");
output <- c(output, "\r\n");
}
}
BufMax <- (inputs$BufferEnd - inputs$BufferStart) / inputs$BufferStep + 1
if( inputs$SRType %in% c("HOC3", "RIC3", "BEV3")){
output <- c(output, "Freshwater Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaFlowA, sep=""), width=34), " Gamma B=", inputs$GammaFlowB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$GammaFlowA * inputs$GammaFlowB, sep=""), width=34), " var=", inputs$GammaFlowB * inputs$GammaFlowB * inputs$GammaFlowA, "\r\n");
output <- c(output, " ", inputs$TrndCycM, inputs$TCF1, inputs$TCF2, inputs$TCF3, "\r\n");
output <- c(output, " ", "First Flow index generated <- ", SummaryStats$FirstRanFlow, "\r\n");
output <- c(output, " ", "Average Flow index generated <- ", SummaryStats$AveRanFlow / inputs$NYears / inputs$NRuns / inputs$BufMax, "\r\n");
output <- c(output, "\r\n");
}
if( inputs$SRType %in% c("HOC4", "RIC4", "BEV4")){
output <- c(output, "Marine Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaMarA, sep=""), width=34), " Gamma B=", inputs$GammaMarB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$MarAve, sep=""), width=34), " st dev=", inputs$MarSD, "\r\n");
output <- c(output, " ", inputs$TrndCycM, inputs$TCM1, inputs$TCM2, inputs$TCM3, "\r\n");
output <- c(output, " ", "First Marine index generated <- ", SummaryStats$FirstRanMarine, "\r\n");
output <- c(output, " ", "Average Marine index generated <- ", mean(SummaryStats$AveRanMarine), "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Freshwater Survival Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaFlowA, sep=""), width=34), " Gamma B=", inputs$GammaFlowB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$FlowAve, sep=""), width=34), " st dev=", inputs$FlowSD, "\r\n");
output <- c(output, " ", inputs$TrndCycF, inputs$TCF1, inputs$TCF2, inputs$TCF3, "\r\n");
output <- c(output, " ", "First Flow index generated <- ", SummaryStats$FirstRanFlow, "\r\n");
output <- c(output, " ", "Average Flow index generated <- ", mean(SummaryStats$AveRanFlow), "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, "Fishery Regime Parameters:", "\r\n");
if(inputs$NumBreakPoints >= 1){
for(Break in 1:inputs$NumBreakPoints){
output <- c(output, " ", "Breakpoint", format(Break, digits=2), "\r\n");
output <- c(output, " ", "HR Below Breakpoint=", inputs$TargetU[Break], "\r\n");
output <- c(output, " ", "Escapement=", inputs$EscpmntBreakPoint[Break], "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, " ", "HR Above BreakPoint=", inputs$TargetU[inputs$NumBreakPoints + 1], "\r\n");
output <- c(output, "\r\n");
}else{
output <- c(output, " ", "Base ER <- ", inputs$TargetU[inputs$BaseRegime], "\r\n");
}
if(inputs$MgmtError == "YES"){
output <- c(output, "Management Variability Parameters:", "\r\n");
output <- c(output, format(paste(" ", "Gamma A=", inputs$GammaMgmtA, sep=""), width=34), " Gamma B=", inputs$GammaMgmtB, "\r\n");
output <- c(output, format(paste(" ", "mean=", inputs$GammaMgmtA * inputs$GammaMgmtB, sep=""), width=34), " var=", inputs$GammaMgmtB * inputs$GammaMgmtB * inputs$GammaMgmtA, "\r\n");
output <- c(output, "\r\n");
}
output <- c(output, "AEQ for age class", "\r\n");
for(Age in inputs$MinAge:inputs$MaxAge){
output <- c(output, " ", "Age", Age, " AEQ =", staticvars$AEQ[Age], "\r\n");
}
output <- c(output, "Recruits At Age 1<- ", staticvars$RecruitsAtAge1, "\r\n");
output <- c(output, "\r\n");
output <- c(output, "Regime Evaluation Parameters: QET <- ", inputs$DL2, "\r\n");
output <- c(output, " ", "Lower Escapement Level (LEL)=", inputs$ECrit, "\r\n");
output <- c(output, " ", "Upper Escapement Level (UEL)=", inputs$ERecovery, "\r\n");
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output <- c(output, " ", "Max Return (under average variability) =", inputs$BSRa * inputs$AveEnv * inputs$BSRb / 2.71828, "\r\n");
}else{
output <- c(output, " ", "Max Return (under average variability) =", inputs$BSRb * inputs$AveEnv, "\r\n");
}
output <- c(output, "\r\n");
output <- c(output, "SUMMARY STATISTICS", "\r\n");
output <- c(output, "All statistics are averaged over repetitions", "\r\n");
output <- c(output, " . . . . __________Escapement___________", "\r\n");
output <- c(output, ". b Total-Exploit.-Rate . . #fish %runs %yrs %runs 1st LastYrs pop_size.", "\r\n");
output <- c(output, ". param. TgtER CYrER BYrER Mort. extnct <LEL end>UEL Year Ave. at_equil. ", "\r\n");
Buffer <- inputs$BufferStart
for(BufNum in 1:inputs$BufMax){
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
output <- c(output, format(round(inputs$BSRb * PBff),nsmall=0,width=7));
output <- c(output, format(round(inputs$TargetU[inputs$BaseRegime] * EBff, digits=2),nsmall=2,width=7));
output <- c(output, format(round(SummaryStats$AvgCaHR[BufNum], digits=3),nsmall=3,width=7));
output <- c(output, format(round(SummaryStats$BufAvgBYrHR[BufNum], digits=3),nsmall=3,width=7));
output <- c(output, format(round(SummaryStats$AvgAEQMort[BufNum]),nsmall=0,width=6));
output <- c(output, format(round(100 * SummaryStats$PropExt[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(100 * SummaryStats$AvgECrit[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(100 * SummaryStats$PropRec[BufNum],digits=1),nsmall=1,width=6));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, 1]),nsmall=0,width=7));
output <- c(output, format(round(SummaryStats$AvgEscpmnt[BufNum, inputs$NYears]),nsmall=0,width=6));
if(inputs$StepFunc == "POP"){
# njs MxR adjusted for average environmental conditions
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
output <- c(output, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime])),nsmall=0,width=8));
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output <- c(output, format(round(PBff * inputs$BSRb * log(inputs$AveEnv * inputs$BSRa * (1 - inputs$TargetU[inputs$BaseRegime]))),nsmall=0,width=8));
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
output <- c(output, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime] - ((1 / inputs$BSRa) * inputs$AveEnv))),nsmall=0,width=8));
}
}
output <- c(output, "\r\n");
Buffer <- Buffer + inputs$BufferStep
} #for(BufNum in 1:inputs$BufMax)
### this creates a table of just the summary columns:
output.tabulated <- data.frame(beta=numeric(), TgtER=numeric(), CYrER=numeric(), BYrER=numeric(), fishMort=numeric(), percent.runs.extinct=numeric(), percent.years.belowLEL=numeric(), percent.years.aboveUEL=numeric(), first.year=numeric(), lastyears.avg=numeric())
Buffer <- inputs$BufferStart
for(BufNum in 1:inputs$BufMax){
if(inputs$StepFunc == "POP"){
PBff <- Buffer
EBff <- 1
}
if(inputs$StepFunc == "ER"){
PBff <- 1
EBff <- Buffer
}
output.row <- cbind(format(round(inputs$BSRb * PBff),nsmall=0,width=7),
format(round(inputs$TargetU[inputs$BaseRegime] * EBff, digits=2),nsmall=2,width=7),
format(round(SummaryStats$AvgCaHR[BufNum], digits=3),nsmall=3,width=7),
format(round(SummaryStats$BufAvgBYrHR[BufNum], digits=3),nsmall=3,width=7),
format(round(SummaryStats$AvgAEQMort[BufNum]),nsmall=0,width=6),
format(round(100 * SummaryStats$PropExt[BufNum],digits=1),nsmall=1,width=6),
format(round(100 * SummaryStats$AvgECrit[BufNum],digits=1),nsmall=1,width=6),
format(round(100 * SummaryStats$PropRec[BufNum],digits=1),nsmall=1,width=6),
format(round(SummaryStats$AvgEscpmnt[BufNum, 1]),nsmall=0,width=7),
format(round(SummaryStats$AvgEscpmnt[BufNum, inputs$NYears]),nsmall=0,width=6))
if(inputs$StepFunc == "POP"){
# njs MxR adjusted for average environmental conditions
if(inputs$SRType %in% c("HOC2", "HOC3", "HOC4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime])),nsmall=0,width=8))
}
if(inputs$SRType %in% c("RIC2", "RIC3", "RIC4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * log(inputs$AveEnv * inputs$BSRa * (1 - inputs$TargetU[inputs$BaseRegime]))),nsmall=0,width=8))
}
if(inputs$SRType %in% c("BEV2", "BEV3", "BEV4")){
output.row <- c(output.row, format(round(PBff * inputs$BSRb * inputs$AveEnv * (1 - inputs$TargetU[inputs$BaseRegime] - ((1 / inputs$BSRa) * inputs$AveEnv))),nsmall=0,width=8))
}
}#if(inputs$StepFunc == "POP"){
output.row <- as.numeric(output.row)
output.tabulated <- rbind(output.tabulated, output.row)
Buffer <- Buffer + inputs$BufferStep
} #for(BufNum in 1:inputs$BufMax)
colnames(output.tabulated) <- c("beta", "TgtER", "CYrER", "BYrER", "fishMort", "percent.runs.extinct", "percent.years.belowLEL", "percent.years.aboveUEL", "first.year", "lastyears.avg")
### end table creation
if(writefile) cat(file=file, output)
invisible(list(output=output,output.tabulated=output.tabulated ))
}#END SaveSummary
#' @title Title
#'
#' @param Year
#' @param YearStats
#'
#' @return
#' @keywords internal
#'
#' @examples
SaveYearData <- function(Year, YearStats){
if(YearStats$TotAEQMort[Year] + YearStats$TotEscpmnt[Year] > 0.000000001){
HR <- YearStats$TotAEQMort[Year] / (YearStats$TotAEQMort[Year] + YearStats$TotEscpmnt[Year])
}else{
HR <- 0
}
#Update AEQMort to hold total
#in VB code Escpmnt (a vector for each age) is redefined to be TotEscpmnt (sum over ages)
#unclear why
# YearStats$AEQMort[Year] <- YearStats$TotAEQMort
# YearStats$Escpmnt[Year] <- YearStats$TotEscpmnt
YearStats$CalendarHR[Year] <- HR
# YearStats$Age2AEQMort[Year] <- YearStats$AEQMort[Year,2]
# YearStats$Age3AEQMort[Year] <- YearStats$AEQMort[Year,3]
# YearStats$Age4AEQMort[Year] <- YearStats$AEQMort[Year,4]
# YearStats$Age5AEQMort[Year] <- YearStats$AEQMort[Year,5]
# YearStats$Age2Escpmnt[Year] <- YearStats$Escpmnt[Year,2]
# YearStats$Age3Escpmnt[Year] <- YearStats$Escpmnt[Year,3]
# YearStats$Age4Escpmnt[Year] <- YearStats$Escpmnt[Year,4]
# YearStats$Age5Escpmnt[Year] <- YearStats$Escpmnt[Year,5]
return(YearStats)
}#END SaveYearData
#' @title Title
#'
#' @param BaseName
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
SetOutFileNames <- function(BaseName, inputs){
# ***** SetOutFileNames ******
# Used by GetOutFiles and also by GetCommandLine
# ***********************************************************************
# NoExtensionName 'stores path and name, but not extension
# Position 'position of final "\" in pathname
# PathName 'string containing the path (including last "\")
# InFileNameOnly 'name (excl path) of input file
# SET OUTPUT FILES
inputs$OutFileSum <- paste(BaseName,".sum",sep="")
inputs$OutFileEsc <- paste(BaseName,".esc",sep="")
inputs$OutFileByr <- paste(BaseName,".byr",sep="")
inputs$OutFileLog <- paste(BaseName,".log",sep="")
inputs$InFileCopy <- paste(BaseName,".rav",sep="")
return(inputs)
}#END SetOutFileNames
#' @title Title
#'
#' @param inputs
#'
#' @return
#' @keywords internal
#'
#' @examples
SetupSummaryStats <- function(inputs){
SummaryStats <- list()
#most get set to 0 because I add info/NRuns each rep
#mins and maxs get NA because I use min(a,b,na.rm=TRUE) to ignore initial setting
SummaryStats[["AvgEscpmnt"]] <- matrix(0, inputs$BufMax, inputs$NYears)
for(el in c("MinEscpmnt","MaxEscpmnt"))
SummaryStats[[el]] <- matrix(NA, inputs$BufMax, inputs$NYears)
SummaryStats[["AvgBYrHR"]] <- matrix(0, inputs$BufMax, inputs$NYears - (inputs$MaxAge - inputs$MinAge))
for(el in c("MinBYrHR","MaxBYrHR"))
SummaryStats[[el]] <- matrix(NA, inputs$BufMax, inputs$NYears - (inputs$MaxAge - inputs$MinAge))
for(el in c("AvgAEQMort","AvgECrit","AvgCaHR", "BufAvgBYrHR", "PropExt", "PropRec", "BufAvgBYrHR", "AveRanFlow", "AveRanMarine"))
SummaryStats[[el]] <- rep(0, inputs$BufMax)
return(SummaryStats)
}#END SetupSummaryStats
#' @title Title
#'
#' @param t
#' @param y
#' @param x
#' @param z
#'
#' @return
#' @keywords internal
#'
#' @examples
Trend <- function(t, y, x, z){
# ***** FUNCTION Trend *****
# Compute variable x with trend t for year y
# NJS: created 7/9/02
# Function Trend(t As Double, y, x As Double, z) As Double
# t is trend rate, y is time increment, x is first value, z is type of trend
#
# If z <- 0 Then
# Trend <- x * (1 + t) ^ y
# ElseIf z <- 1 Then
# Trend <- x + (y * t)
# Else
# If Trend < 0 Then Trend <- 0
#
# Print "Unknown trend/cycle function"
# Stop
# End If
#
# End Function
# t is trend rate, y is time increment, x is first value, z is type of trend 0 or 1
# EEH: changed to work with vectors of y
if(!(z == 0 | z==1)) stop("Unknown trend/cycle function")
if(z == 0){
val <- x * (1 + t) ^ y
}else{
val <- x + (y * t)
}
val[val < 0] <- 0
return(val)
}#END Trend
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.