#' Title
#'
#' @param vrap.martin
#'
#' @return
#' @export
#'
#' @examples
calcRERfrom_runSimulations <- function(vrap.martin){
totEsc <- vrap.martin$totEsc
#rows=er levels, cols=sims, 3rd dim=years
lel <- vrap.martin$input$ECrit
uel <- vrap.martin$input$ERecovery
NYears <- vrap.martin$input$NYears
NRuns <- vrap.martin$input$NRuns
EndAv <- vrap.martin$input$EndAv
targetER.vec <- vrap.martin$targetER.vec
#calc lel props when er=0
#making the assumption that row 1 will always represent when er=0 (risky)
prop.at.zeroER <- sum(totEsc[1,,]<lel)/(NYears*NRuns)
#calc lel props
prop.lel <- apply(totEsc,1, FUN = function(x) sum(x[,]<lel)) / (NYears*NRuns) + prop.at.zeroER
#calc uel props
uel.count <- apply(totEsc,c(1,2), FUN = function(esc.vec) {
year.index <- length(esc.vec):(length(esc.vec)-EndAv+1)
geo.mean <- exp(mean(log(esc.vec[year.index])))
uel.count <- ifelse(geo.mean>uel, 1, 0)
return(uel.count)
} )
prop.uel <- apply(uel.count, 1, FUN=function(x) sum(x)/NRuns)
er.df <- data.frame(er.target=targetER.vec, prop.lel=prop.lel, prop.uel=prop.uel)
er.qualify.df <- er.df[er.df$prop.lel<prop.at.zeroER+.05 & er.df$prop.uel > 0.8,]
rer <- max(er.qualify.df$er.target)
return(list(er.df=er.df, er.qualify.df=er.qualify.df, rer=rer, runtime=vrap.martin$time.run))
}#END calcRERfrom_runSimulations
#' Title
#'
#' @param rav2.filename
#'
#' @return
#' @export
#'
#' @examples
runMartinVRAP <- function(rav2.filename){
time.start <- Sys.time()
input.list <- VRAP2::read_rav2(rav2.filename)
input <- input.list$inputs.vars
#need to convert some variable names to satisfy his code
input$ERnumSteps <- input$BufMax
# parametere for setting exploitation rate start, stop and steps
input$baseER <- input$TargetU
input$stepSize <- input$BufferStep
input$ERmin <- input$BufferStart
input$ERmax <- input$BufferEnd
# derived ER parameters that make more sense
input$ERstepSize <- input$baseER*input$stepSize # stepsize in terms of the exploitation rate
input$ERnumSteps <- round(input$ERmax/input$stepSize) # number of steps
input$ERstart = input$ERmin
input$numSims <- input$NRuns
input$numYears <- input$NYears
input$mat <- input$MatRate
input$mort <- input$NatMort
input$initPop <- input$CohortStart
input$managementError <- ifelse(toupper(input$MgmtError)=="YES",TRUE,FALSE)
input$manageErrorA <- input$GammaMgmtA
input$manageErrorB <- input$GammaMgmtB
input$HRpt <- input$PTU
input$HRt <- input$MatU
input$prod <- input$BSRa
input$cap <- input$BSRb
input$SRerrorA <- input$SRErrorA
input$SRerrorB <- input$SRErrorB
esc.output <- runSimulations(input = input)
time.run <- Sys.time()-time.start
esc.output$time.run <- time.run
saveRDS(object = esc.output, file = paste0(tools::file_path_sans_ext(rav2.filename), ".rds"))
}
#' Title
#'
#' @param input
#' @param rngSeed
#'
#' @return
#' @export
#'
#' @examples
runSimulations <- function(input,rngSeed=NULL){
# initialize output array
# browser()
totEsc <- array(NA,dim=c(input$ERnumSteps+1,input$numSims,input$numYears))
HRscale <- 1 # multiply this times the harvest rates.
# calcualte the target exploitation rates based on start, step size, and steps
targetER <- input$ERstart + input$ERstepSize * (0:input$ERnumSteps)
# calculate AEQ and recruitsFromAgeOneFish
AEQ <- c(0,0,0,0,input$mat[5])
for(age in 4:2){
AEQ[age] <- input$mat[age] + (1-input$mort[age+1]) * (1 - input$mat[age]) * AEQ[age+1]
}
recruitsFromAgeOneFish <- (1-input$mort[1])*(1-input$mort[2])*AEQ[2]
for(ERind in 1:(input$ERnumSteps+1)){
print(paste("============= target ER =",targetER[ERind]))
for(sim in 1:input$numSims){ # loop through 1000 25 yr simulations
SRerror <- rnorm(1, 0, sd=sqrt(input$SRerrorB)) # not currently used
Cohort <- input$initPop # initialize population
for(year in 1:input$numYears){ # loop through 25 year simulation
# apply natural mortality
Cohort <- Cohort*(1-input$mort)
# generate management error
actualER <- targetER[ERind]
if(input$managementError){
repeat{
#mf modifed to resample:
ErrorScale <- rgamma(1, input$manageErrorA, scale=input$manageErrorB)
if(actualER * ErrorScale< 1 ) {break} }
actualER <- actualER * ErrorScale
#actualER <- min(actualER * rgamma(1, input$manageErrorA, scale=input$manageErrorB),1)
}
# loop to achieve target exploitation rate unless targetER=0
if(actualER==0){
HRptAdj <- 0
HRtAdj <- 0
}else{
numTrys <- 1
lastAEQmort <- 99
repeat{
# adjust preterminal and terminal fishing rates
HRptAdj <- input$HRpt*HRscale
HRtAdj <- input$HRt*HRscale
# can't be larger than 1
HRptAdj[HRptAdj>1] <- 1
HRtAdj[HRtAdj>1] <- 1
# calculate AEQ fishing mortality, escapement, and the exploitation rate
AEQmort <- Cohort*(HRptAdj*AEQ + (1-HRptAdj)*input$mat*HRtAdj)
Escpmnt <- Cohort*(1-HRptAdj)*(1-HRtAdj)*input$mat
totAEQmort <- sum(AEQmort)
totEscpmnt <- sum(Escpmnt)
ER <- totAEQmort/(totAEQmort+totEscpmnt)
# calculate the error rate (how far the actual ER is from the target)
ERerror <- abs(ER-actualER)/actualER
# exit loop if you are close enough OR other criteria are met. Otherwise adjust HRscale.
if(totAEQmort+totEscpmnt < 1 | totAEQmort==0 | numTrys > 100 | totAEQmort==lastAEQmort){
# cat(paste("Target ER = ",targetER[ERind]," Sim = ",sim," Year = ",year,
# " goal - actual = ",round(actualER,3)," - ",round(ER,3),
# " HRscale = ",round(HRscale,3)," numTrys = ",numTrys,
# " totEsc = ",round(totEscpmnt,1)," totAEQmort = ",round(totAEQmort,1),"\n",sep=""))
break
}else if(ERerror < 0.001) break
else HRscale <- HRscale*actualER/ER
numTrys <- numTrys+1
lastAEQmort <- totAEQmort
}
}
# calculate new cohort
newCohort <- Cohort*(1-HRptAdj)*(1-input$mat)
Escpmnt <- Cohort*(1-HRptAdj)*(1-HRtAdj)*input$mat
Escpmnt[Escpmnt < 1] <- 0
# calculate adult escapement
adultEscapement <- sum(Escpmnt[3:5])
# age the cohort
Cohort[2:5] <- newCohort[1:4]
# now fill in age 1 fish using the spawner-recruit function.
AEQrecruits <- input$prod * adultEscapement * exp(-adultEscapement / input$cap)
SRerror <- rgamma(1, input$SRerrorA, scale=input$SRerrorB)
Cohort[1] <- AEQrecruits*SRerror/recruitsFromAgeOneFish
totEsc[ERind,sim,year] <- sum(Escpmnt[2:5])
}
}
}
list(input=input, targetER.vec=targetER, totEsc=totEsc)
}#END runSimulations
# simFish <- function(){
# #include <Rcpp.h>
# using namespace Rcpp;
#
# // [[Rcpp::export]]
# NumericMatrix simFish(int numSims, int numYears,
# double targetER,
# bool managementError,
# double manageErrorA, double manageErrorB,
# String errorType,
# double SRerrorA, double SRerrorB,
# NumericVector initPop,
# double prod, double cap,
# NumericVector mat, NumericVector mort,
# NumericVector HRpt, NumericVector HRt,
# NumericVector AEQ){
#
# NumericMatrix totEsc(numSims, numYears);
# NumericVector Cohort(5);
# NumericVector newCohort(5);
# NumericVector HRptAdj(5);
# NumericVector HRtAdj(5);
# NumericVector AEQmort(5);
# NumericVector Escpmnt(5);
#
# double lastAEQmort, HRscale, totAEQmort, totEscpmnt, adultEscapement,
# AEQrecruits, ER, realizedER, ERerror, SRerror, logSRerror;
# bool converged;
# int numTrys;
#
# double recruitsFromAgeOneFish = (1-mort(0))*(1-mort(1))*AEQ(1);
#
# //Rcpp::Rcout << "targetER=" << targetER << "\n";
#
# for(int sim = 0; sim < numSims; sim++) {
# // loop through simulations
# Cohort = clone(initPop);
# // initialize population. Use clone to create a deep copy (i.e don't just copy a pointer)
# logSRerror = rnorm(1, 0, SRerrorA)[0];
# for(int year = 0; year < numYears; year++) { // loop through years
# Cohort = Cohort*(1-mort);
# if(managementError) ER = std::min(targetER * rgamma(1, manageErrorA, manageErrorB)[0],1.0); else ER = targetER;
# //Rcpp::Rcout << ER << ",";
# if(ER==0){
# std::fill(HRptAdj.begin(), HRptAdj.end(), 0); // set all elements to 0
# std::fill(HRtAdj.begin(), HRtAdj.end(), 0); // set all elements to 0
# }else{
# numTrys = 1;
# lastAEQmort = 99;
# converged = false;
# HRscale = 1;
# // Rcpp::Rcout <<"year=" << year << "Cohort=" << Cohort << "\n";
# while(!converged){
# // adjust preterminal and terminal fishing rates
# HRptAdj = pmin(HRpt*HRscale,1);
# HRtAdj = pmin(HRt*HRscale,1);
#
# // calculate AEQ fishing mortality, escapement, and the exploitation rate
# AEQmort = Cohort*(HRptAdj*AEQ + (1-HRptAdj)*mat*HRtAdj);
# Escpmnt = Cohort*(1-HRptAdj)*(1-HRtAdj)*mat;
# totAEQmort = sum(AEQmort);
# totEscpmnt = sum(Escpmnt);
# realizedER = totAEQmort/(totAEQmort+totEscpmnt);
# // calculate the error rate (how far the actual ER is from the target)
# // Rcpp::Rcout << "year=" << year << "ER=" << ER << "\n";
# ERerror = std::abs(ER-realizedER)/ER;
# // exit loop if you are close enough OR other criteria are met
#
# //Rcpp::Rcout << "HRptAdj,HRtAdj=" << HRptAdj << "," << HRtAdj << "\n";
# //Rcpp::Rcout << "numTrys=" << numTrys << " HRscale=" << HRscale << "\n";
# //Rcpp::Rcout << "totAEQmort=" << totAEQmort << " totEscpmnt=" << totEscpmnt << "\n";
# //Rcpp::Rcout << "ER=" << ER << " realizedER=" << realizedER << "\n";
# //Rcpp::Rcout << "ERerror=" << ERerror << "\n";
#
# if((totAEQmort+totEscpmnt < 1) || (totAEQmort==0) || (numTrys > 100) || (totAEQmort==lastAEQmort)){
# converged = true;
# }else if(ERerror < 0.001){
# converged = true;
# }else{
# HRscale = HRscale*ER/realizedER;
# }
# numTrys = numTrys+1;
# lastAEQmort = totAEQmort;
# }
# }
# // calculate new cohort
# newCohort = Cohort*(1-HRptAdj)*(1-mat);
# // Rcpp::Rcout << "newCohort=" << newCohort << "\n";
# Escpmnt = pmax(Cohort*(1-HRptAdj)*(1-HRtAdj)*mat,0);
# // calculate adult escapement
# adultEscapement = Escpmnt(2) + Escpmnt(3) + Escpmnt(4);
# // age the cohort
# for(int ageInd = 0; ageInd < 4; ageInd++) Cohort(ageInd+1) = newCohort(ageInd);
# // now fill in age 1 fish using the spawner-recruit function.
# AEQrecruits = prod * adultEscapement * exp(-adultEscapement / cap);
# if(errorType=="gamma"){
# SRerror = rgamma(1,SRerrorA,SRerrorB)[0];
# }else if(errorType=="logNormal"){
# // SRerrorA = lognormal sd, SRerrorB = autocorrelation
# logSRerror = SRerrorB*logSRerror + sqrt(1-pow(SRerrorB,2.0))*rnorm(1, 0, SRerrorA)[0];
# SRerror = exp(logSRerror);
# }
# Cohort(0) = AEQrecruits*SRerror/recruitsFromAgeOneFish;
# totEsc(sim,year) = Escpmnt(1) + Escpmnt(2) + Escpmnt(3) + Escpmnt(4);
# }
# }
# return totEsc;
# }
#
# }#END simFish
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.