R/VRAP_modified_functions.R

#' 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
MichaelFolkes/VRAP2 documentation built on June 9, 2019, 6:27 p.m.