pollock and codGB problems

[mle2718]

I am working in the EconOnly2 branch, which was recently synched with master. edit: The codGOM model runs just fine by setting

stockExclude <- c('haddockGB', 'codGB', 'pollock', 'yellowtailflounderGB')

When I try to run a standardFisheries model with any of these values of stockExclude:

stockExclude <- c('haddockGB', 'pollock', 'yellowtailflounderGB')
stockExclude <- c('haddockGB', 'codGB', 'yellowtailflounderGB')

I get the following message:

Error in eval(substitute(expr), data, enclos = parent.frame()) : 
  object 'Rhat' not found

The traceback looks like:

16: eval(substitute(expr), data, enclos = parent.frame())
15: eval(substitute(expr), data, enclos = parent.frame())
14: with.default(stockEnv, {
        if ("rho" %in% names(par)) {
            if (par["rho"] < -1 || par["rho"] > 1) {
                stop("Rfun: par[rho] must be between 0 and 1")
            }
            if (is.null(R_ym1) || is.null(Rhat_ym1)) {
                stop("Rfun: if rho is given R_ym1 and Rhat_ym1 must be provided")
            }
            if (is.na(R_ym1) || is.na(Rhat_ym1)) {
                warning("Rfun: NA encountered in R_ym1 or Rhat_y-1 -- autocorrelation\n               error not included")
                R_ym1 <- 1
                Rhat_ym1 <- 1
            }
        }
        else {
            par["rho"] <- 0
            R_ym1 <- 1
            Rhat_ym1 <- 1
        }
        if (type == "BH") {
            Rhat <- (par["a"] * SSB)/(1 + (par["b"] * SSB)) * exp(par["g"] * 
                TAnom_y)
            Rhat <- Rhat * 1000
        }
        else if (type == "BHSteep") {
            if (is.na(par["beta1"])) {
                par["beta1"] <- 0
            }
            if (is.na(par["beta2"])) {
                par["beta2"] <- 0
            }
            if (is.na(par["beta3"])) {
                par["beta3"] <- 0
            }
            Rhat <- with(as.list(par), {
                gamma <- -0.5 * log((1 - 0.2)/(h - 0.2) - 1) + beta1 * 
                    TAnom_y
                hPrime <- 0.2 + (1 - 0.2)/(1 + exp(-2 * gamma))
                R0Prime <- R0 * exp(beta2 * TAnom_y)
                num <- 4 * hPrime * (SSB/(SSBRF0))
                den <- ((1 - hPrime) + (5 * hPrime - 1) * (SSB/(R0Prime * 
                    SSBRF0)))
                z <- num/den * exp(beta3 * TAnom_y)
                return(z)
            })
        }
        else if (type == "HS") {
            if (stock[[i]]$stockName == "codGOM") {
                Rhat <- with(as.list(par), {
                    SSBhinge <- SSB_star
                    if (type2 == "True") {
                      assess_vals <- get_HistAssess(stock = stock[[i]])
                      if (SSB >= SSBhinge) {
                        pred <- cR * remp(1, tail(as.numeric(assess_vals$assessdat$R), 
                          Rnyr))
                      }
                      else if (SSB < SSBhinge) {
                        pred <- cR * (SSB/SSBhinge) * remp(1, tail(as.numeric(assess_vals$assessdat$R), 
                          Rnyr))
                      }
                    }
                    else {
                      if (SSB >= SSBhinge) {
                        pred <- cR * remp(1, tail(as.numeric(R_est), 
                          Rnyr))
                      }
                      else if (SSB < SSBhinge) {
                        pred <- cR * (SSB/SSBhinge) * remp(1, tail(as.numeric(R_est), 
                          Rnyr))
                      }
                    }
                    return(pred)
                })
            }
            if (stock[[i]]$stockName == "haddockGB") {
                Rhat <- with(as.list(par), {
                    if (type2 == "True") {
                      assess_vals <- get_HistAssess(stock = stock[[i]])
                      pred <- remp(1, tail(as.numeric(assess_vals$assessdat$R, 
                        20)))
                    }
                    else {
                      pred <- remp(1, as.numeric(R_est))
                    }
                    return(pred)
                })
            }
        }
        ac <- par["rho"] * log(R_ym1/Rhat_ym1)
        rc <- rnorm(1, mean = 0, sd = pe_R)
        R <- Rhat * exp(ac + rc)
        out <- list(Rhat = unname(Rhat), R = unname(R))
        return(out)
    })
13: with(stockEnv, {
        if ("rho" %in% names(par)) {
            if (par["rho"] < -1 || par["rho"] > 1) {
                stop("Rfun: par[rho] must be between 0 and 1")
            }
            if (is.null(R_ym1) || is.null(Rhat_ym1)) {
                stop("Rfun: if rho is given R_ym1 and Rhat_ym1 must be provided")
            }
            if (is.na(R_ym1) || is.na(Rhat_ym1)) {
                warning("Rfun: NA encountered in R_ym1 or Rhat_y-1 -- autocorrelation\n               error not included")
                R_ym1 <- 1
                Rhat_ym1 <- 1
            }
        }
        else {
            par["rho"] <- 0
            R_ym1 <- 1
            Rhat_ym1 <- 1
        }
        if (type == "BH") {
            Rhat <- (par["a"] * SSB)/(1 + (par["b"] * SSB)) * exp(par["g"] * 
                TAnom_y)
            Rhat <- Rhat * 1000
        }
        else if (type == "BHSteep") {
            if (is.na(par["beta1"])) {
                par["beta1"] <- 0
            }
            if (is.na(par["beta2"])) {
                par["beta2"] <- 0
            }
            if (is.na(par["beta3"])) {
                par["beta3"] <- 0
            }
            Rhat <- with(as.list(par), {
                gamma <- -0.5 * log((1 - 0.2)/(h - 0.2) - 1) + beta1 * 
                    TAnom_y
                hPrime <- 0.2 + (1 - 0.2)/(1 + exp(-2 * gamma))
                R0Prime <- R0 * exp(beta2 * TAnom_y)
                num <- 4 * hPrime * (SSB/(SSBRF0))
                den <- ((1 - hPrime) + (5 * hPrime - 1) * (SSB/(R0Prime * 
                    SSBRF0)))
                z <- num/den * exp(beta3 * TAnom_y)
                return(z)
            })
        }
        else if (type == "HS") {
            if (stock[[i]]$stockName == "codGOM") {
                Rhat <- with(as.list(par), {
                    SSBhinge <- SSB_star
                    if (type2 == "True") {
                      assess_vals <- get_HistAssess(stock = stock[[i]])
                      if (SSB >= SSBhinge) {
                        pred <- cR * remp(1, tail(as.numeric(assess_vals$assessdat$R), 
                          Rnyr))
                      }
                      else if (SSB < SSBhinge) {
                        pred <- cR * (SSB/SSBhinge) * remp(1, tail(as.numeric(assess_vals$assessdat$R), 
                          Rnyr))
                      }
                    }
                    else {
                      if (SSB >= SSBhinge) {
                        pred <- cR * remp(1, tail(as.numeric(R_est), 
                          Rnyr))
                      }
                      else if (SSB < SSBhinge) {
                        pred <- cR * (SSB/SSBhinge) * remp(1, tail(as.numeric(R_est), 
                          Rnyr))
                      }
                    }
                    return(pred)
                })
            }
            if (stock[[i]]$stockName == "haddockGB") {
                Rhat <- with(as.list(par), {
                    if (type2 == "True") {
                      assess_vals <- get_HistAssess(stock = stock[[i]])
                      pred <- remp(1, tail(as.numeric(assess_vals$assessdat$R, 
                        20)))
                    }
                    else {
                      pred <- remp(1, as.numeric(R_est))
                    }
                    return(pred)
                })
            }
        }
        ac <- par["rho"] * log(R_ym1/Rhat_ym1)
        rc <- rnorm(1, mean = 0, sd = pe_R)
        R <- Rhat * exp(ac + rc)
        out <- list(Rhat = unname(Rhat), R = unname(R))
        return(out)
    }) at get_recruits.R#51
12: get_recruits(type = type, type2 = type2, par = Rpar, SSB = SSB, 
        TAnom = TAnom, pe_R = sdR, R_ym1 = 1, block = "late", Rhat_ym1 = 1, 
        R_est = parpop$R) at Rfun_BmsySim.R#25
11: Rfun(type = stockEnv$R_typ, parpop = parpop, parenv = parenv, 
        parmgt = parmgt, SSB = c(N[y - 1, ]) %*% c(parpop$waa * parpop$mat), 
        sdR = stockEnv$pe_R, TAnom = Tanom[y], Rest = Rest) at get_proj.R#156
10: get_proj(type = "FREF", parmgt = parmgt, parpop = parpop, parenv = parenvTemp, 
        Rfun = Rfun, F_val = candF[x], ny = 200, stReportYr = 2, 
        stockEnv = stockEnv) at get_FBRP.R#93
9: FUN(X[[i]], ...)
8: lapply(1:length(candF), function(x) {
       get_proj(type = "FREF", parmgt = parmgt, parpop = parpop, 
           parenv = parenvTemp, Rfun = Rfun, F_val = candF[x], ny = 200, 
           stReportYr = 2, stockEnv = stockEnv)
   }) at get_FBRP.R#92
7: get_FBRP(parmgt = parmgt, parpop = parpop, parenv = parenv, Rfun_lst = Rfun_lst, 
       stockEnv = stock) at get_burnF.R#27
6: get_burnF(stockT) at get_popInit.R#88
5: eval(substitute(expr), e)
4: eval(substitute(expr), e)
3: within.list(stock, {
       stockT <- stock
       stockT$R_mis <- FALSE
       burnFmsy <- get_burnF(stockT)
       burnFmean <- burnFmsyScalar * burnFmsy
       F_full[(fyear + 1):fmyearIdx] <- rlnorm(fmyearIdx - (fyear + 
           1) + 1, log(burnFmean), burnFsd)
   })
2: within(stock, {
       stockT <- stock
       stockT$R_mis <- FALSE
       burnFmsy <- get_burnF(stockT)
       burnFmean <- burnFmsyScalar * burnFmsy
       F_full[(fyear + 1):fmyearIdx] <- rlnorm(fmyearIdx - (fyear + 
           1) + 1, log(burnFmean), burnFsd)
   }) at get_popInit.R#85
1: get_popInit(stock[[i]])

I believe this problem is coming at this line in runSim.R

 # Initialize stocks and determine burn-in F
    for(i in 1:nstock){
      stock[[i]] <- get_popInit(stock[[i]])
    }

[samtruesdell]

So is this problem related to codGOM then? Will that run by itself?

[mle2718]

@samtruesdell sorry, I should have clarified that the codGOM works. I will edit

[lkerr]

Is this still an issue after you implemented: histAssess <- FALSE in set_om_parameters_global.R to use parameters burnFmsyScalar and burnFsd

This fix addressed the burn-in for F.

[mle2718]

@lkerr: unfortunately, not fixed. No change in the error message.

[lkerr]

Jerelle is looking into the different stock parameter files to see what might be happening here-work in progress.

[jerellejesse]

To me it looks like the HS recruitment was never set up for pollock in get_recruits., only codGOM and haddock are listed. Also the SSB hinge point wasn't located in the assessment materials for the stock parameter file. I can work on this if it is the problem.

[mle2718]

@jerellejesse: That would be super helpful to me. My guess is that there will be similar issues for yellowtailGB.

My 'wrinkle' is that I'd like to simultanously run as many of the 5 stocks as I can. Which necessitates passing in the same options from the mproc.csv file and the set_om_parameters_global.R.

Which may be both tricky and unrealistic.

[mdmazur]

Yes, HS recruitment was not set up for any stock aside from codGOM and haddock. I provided a fix for this in the new branch Fix_pollock_and_codGB_problems that is the same as EconOnly2 but has a change in get_recruits so that the hockey-stick option will run for any stock. Before it was only running for GOM cod and GB haddock. Haddock's hockey-stick recruitment is different since there is not a hinge point. So, I created a if statement for haddock but then all other stocks follow the hockey-stick recruitment form used for GOM cod. All stocks need SSB_star (which is the hinge point), cR (a conversion factor) and Rnyr (the number of years of recruitment used in the empirical distribution function). Those should be specified in the scripts in stockParameters. One other issue I ran into while running this was that if the minimum catch constraint is on, all stocks will need a csv with discard data: data/data_raw/AssessmentHistory/'stockname'_Discard.csv. The minimum catch constraint is based on the last ten years of discards. Right now only GOM cod and GB haddock have these csvs. We will also need an if statement for those stocks in get_nextF around line 220. Or those lines could be changed so that the stock name is inserted into the csv name to be loaded.

[jerellejesse]

@mdmazur That's great, thank you! I was looking into this, but had some questions, so it will be great to see what you did. The minimum catch constraint was on my list also, so I can help with that.

[mle2718]

@mdmazur Thanks!! I've merge this into ``EconOnly2'' (#225) and I'll see how it works.

[mle2718]

@mdmazur and @jerellejesse : After merging, I can run GOM Cod and GB haddock.

# Available stocks: haddockGB, codGOM, codGB, pollock, yellowtailflounderGB
stockExclude <- c('codGB', 'pollock', 'yellowtailflounderGB')

I have a few more bit to work on for the Econ model, I'd say I've got about a week of work on it, but I don't know when I can fit it in.

It looks like the other stocks just needs a handful of parameters, which seems like an easy fix, but might also be a painful fix.

[mle2718]

ah -- spoke too soon. I get this message:

Error in `[<-`(`*tmp*`, r, m, , value = c(80886432.4078435, 80886432.4078435,  : 
  subscript out of bounds
Called from: eval(substitute(expr), e)

and the traceback is:

5: eval(substitute(expr), e)
4: eval(substitute(expr), e)
3: within.list(stock, {
       omval$N[r, m, ] <- rowSums(J1N)
       omval$SSB[r, m, ] <- SSB
       omval$R[r, m, ] <- R
       omval$F_full[r, m, ] <- F_full
       omval$sumCW[r, m, ] <- sumCW
       omval$OFdStatus[r, m, ] <- OFdStatus
       omval$mxGradCAA[r, m, ] <- mxGradCAA
       omval$F_fullAdvice[r, m, ] <- F_fullAdvice
       omval$ACL[r, m, ] <- ACL
       omval$OFgStatus[r, m, ] <- OFgStatus
       omval$SSB_cur[r, m, ] <- SSB_cur
       omval$natM[r, m, ] <- natM
       omval$annPercentChange[r, m, ] <- NA
       meanSizeCN <- sapply(1:nrow(CN), function(x) laa[x, ] %*% 
           paaCN[x, ])
       omval$meanSizeCN[r, m, ] <- meanSizeCN
       meanSizeIN <- sapply(1:nrow(CN), function(x) laa[x, ] %*% 
           paaIN[x, ])
       omval$meanSizeIN[r, m, ] <- meanSizeIN
       omval$FPROXY[r, m, ] <- RPmat[, 1]
       omval$SSBPROXY[r, m, ] <- RPmat[, 2]
       omval$FPROXYT[r, m, ] <- RPmat[, 3]
       omval$SSBPROXYT[r, m, ] <- RPmat[, 4]
       omval$FRATIO[r, m, y] <- stock$res$F.report[length(stock$res$F.report)]/RPmat[, 
           1][y]
       omval$SSBRATIO[r, m, y] <- stock$res$SSB[length(stock$res$SSB)]/RPmat[, 
           2][y]
       if (mproc[m, "rhoadjust"] == "TRUE" & y > fmyearIdx & Mohns_Rho_SSB[y] > 
           0.15) {
           omval$FRATIO[r, m, y] <- (stock$res$F.report[length(stock$res$F.report)]/(Mohns_Rho_F[y] + 
               1))/RPmat[, 1][y]
           omval$SSBRATIO[r, m, y] <- (stock$res$SSB[length(stock$res$SSB)]/(Mohns_Rho_SSB[y] + 
               1))/RPmat[, 2][y]
       }
       omval$FRATIOT[r, m, y] <- stock$F_full[y]/RPmat[, 3][y]
       omval$SSBRATIOT[r, m, y] <- stock$SSB[y]/RPmat[, 4][y]
       if (y == nyear) {
           if (nyear > length(cmip_dwn$YEAR)) {
               nprologueY <- nyear - length(cmip_dwn$YEAR)
               prologueY <- (cmip_dwn$YEAR[1] - nprologueY):(cmip_dwn$YEAR[1] - 
                   1)
               yrs <- c(prologueY, cmip_dwn$YEAR)
           }
           else {
               yrs <- rev(rev(cmip_dwn$YEAR)[1:nyear])
           }
           omval$YEAR <- yrs
       }
       omval$relE_qI[r, m, ] <- relE_qI
       omval$relE_qC[r, m, ] <- relE_qC
       omval$relE_selCs0[r, m, ] <- relE_selCs0
       omval$relE_selCs1[r, m, ] <- relE_selCs1
       omval$relE_ipop_mean[r, m, ] <- relE_ipop_mean
       omval$relE_ipop_dev[r, m, ] <- relE_ipop_dev
       omval$relE_R_dev[r, m, ] <- relE_R_dev
       omval$relE_SSB[r, m, ] <- relE_SSB
       omval$relE_N[r, m, ] <- relE_N
       omval$relE_IN[r, m, ] <- relE_IN
       omval$relE_R[r, m, ] <- relE_R
       omval$relE_F[r, m, ] <- relE_F
       omval$conv_rate[r, m, ] <- conv_rate
       omval$Mohns_Rho_SSB[r, m, ] <- Mohns_Rho_SSB
       omval$Mohns_Rho_N[r, m, ] <- Mohns_Rho_N
       omval$Mohns_Rho_F[r, m, ] <- Mohns_Rho_F
       omval$Mohns_Rho_R[r, m, ] <- Mohns_Rho_R
       omval$mincatchcon[r, m, ] <- mincatchcon
       omval$SSBest[y, 1:length(stock$res$SSB)] <- stock$res$SSB
       omval$Fest[y, 1:length(stock$res$SSB)] <- stock$res$F.report
       omval$Catchest[y, 1:length(stock$res$SSB)] <- stock$res$catch.pred
       omval$Rest[y, 1:length(stock$res$SSB)] <- stock$res$N.age[, 
           1]
       if (y == nyear) {
           omval$relTermE_SSB[r, m, ] <- relTermE_SSB
           omval$relTermE_CW[r, m, ] <- relTermE_CW
           omval$relTermE_CW[r, m, ] <- relTermE_CW
           omval$relTermE_IN[r, m, ] <- relTermE_IN
           omval$relTermE_qI[r, m, ] <- relTermE_qI
           omval$relTermE_R[r, m, ] <- relTermE_R
           omval$relTermE_F[r, m, ] <- relTermE_F
       }
   })
2: within(stock, {
       omval$N[r, m, ] <- rowSums(J1N)
       omval$SSB[r, m, ] <- SSB
       omval$R[r, m, ] <- R
       omval$F_full[r, m, ] <- F_full
       omval$sumCW[r, m, ] <- sumCW
       omval$OFdStatus[r, m, ] <- OFdStatus
       omval$mxGradCAA[r, m, ] <- mxGradCAA
       omval$F_fullAdvice[r, m, ] <- F_fullAdvice
       omval$ACL[r, m, ] <- ACL
       omval$OFgStatus[r, m, ] <- OFgStatus
       omval$SSB_cur[r, m, ] <- SSB_cur
       omval$natM[r, m, ] <- natM
       omval$annPercentChange[r, m, ] <- NA
       meanSizeCN <- sapply(1:nrow(CN), function(x) laa[x, ] %*% 
           paaCN[x, ])
       omval$meanSizeCN[r, m, ] <- meanSizeCN
       meanSizeIN <- sapply(1:nrow(CN), function(x) laa[x, ] %*% 
           paaIN[x, ])
       omval$meanSizeIN[r, m, ] <- meanSizeIN
       omval$FPROXY[r, m, ] <- RPmat[, 1]
       omval$SSBPROXY[r, m, ] <- RPmat[, 2]
       omval$FPROXYT[r, m, ] <- RPmat[, 3]
       omval$SSBPROXYT[r, m, ] <- RPmat[, 4]
       omval$FRATIO[r, m, y] <- stock$res$F.report[length(stock$res$F.report)]/RPmat[, 
           1][y]
       omval$SSBRATIO[r, m, y] <- stock$res$SSB[length(stock$res$SSB)]/RPmat[, 
           2][y]
       if (mproc[m, "rhoadjust"] == "TRUE" & y > fmyearIdx & Mohns_Rho_SSB[y] > 
           0.15) {
           omval$FRATIO[r, m, y] <- (stock$res$F.report[length(stock$res$F.report)]/(Mohns_Rho_F[y] + 
               1))/RPmat[, 1][y]
           omval$SSBRATIO[r, m, y] <- (stock$res$SSB[length(stock$res$SSB)]/(Mohns_Rho_SSB[y] + 
               1))/RPmat[, 2][y]
       }
       omval$FRATIOT[r, m, y] <- stock$F_full[y]/RPmat[, 3][y]
       omval$SSBRATIOT[r, m, y] <- stock$SSB[y]/RPmat[, 4][y]
       if (y == nyear) {
           if (nyear > length(cmip_dwn$YEAR)) {
               nprologueY <- nyear - length(cmip_dwn$YEAR)
               prologueY <- (cmip_dwn$YEAR[1] - nprologueY):(cmip_dwn$YEAR[1] - 
                   1)
               yrs <- c(prologueY, cmip_dwn$YEAR)
           }
           else {
               yrs <- rev(rev(cmip_dwn$YEAR)[1:nyear])
           }
           omval$YEAR <- yrs
       }
       omval$relE_qI[r, m, ] <- relE_qI
       omval$relE_qC[r, m, ] <- relE_qC
       omval$relE_selCs0[r, m, ] <- relE_selCs0
       omval$relE_selCs1[r, m, ] <- relE_selCs1
       omval$relE_ipop_mean[r, m, ] <- relE_ipop_mean
       omval$relE_ipop_dev[r, m, ] <- relE_ipop_dev
       omval$relE_R_dev[r, m, ] <- relE_R_dev
       omval$relE_SSB[r, m, ] <- relE_SSB
       omval$relE_N[r, m, ] <- relE_N
       omval$relE_IN[r, m, ] <- relE_IN
       omval$relE_R[r, m, ] <- relE_R
       omval$relE_F[r, m, ] <- relE_F
       omval$conv_rate[r, m, ] <- conv_rate
       omval$Mohns_Rho_SSB[r, m, ] <- Mohns_Rho_SSB
       omval$Mohns_Rho_N[r, m, ] <- Mohns_Rho_N
       omval$Mohns_Rho_F[r, m, ] <- Mohns_Rho_F
       omval$Mohns_Rho_R[r, m, ] <- Mohns_Rho_R
       omval$mincatchcon[r, m, ] <- mincatchcon
       omval$SSBest[y, 1:length(stock$res$SSB)] <- stock$res$SSB
       omval$Fest[y, 1:length(stock$res$SSB)] <- stock$res$F.report
       omval$Catchest[y, 1:length(stock$res$SSB)] <- stock$res$catch.pred
       omval$Rest[y, 1:length(stock$res$SSB)] <- stock$res$N.age[, 
           1]
       if (y == nyear) {
           omval$relTermE_SSB[r, m, ] <- relTermE_SSB
           omval$relTermE_CW[r, m, ] <- relTermE_CW
           omval$relTermE_CW[r, m, ] <- relTermE_CW
           omval$relTermE_IN[r, m, ] <- relTermE_IN
           omval$relTermE_qI[r, m, ] <- relTermE_qI
           omval$relTermE_R[r, m, ] <- relTermE_R
           omval$relTermE_F[r, m, ] <- relTermE_F
       }
   }) at get_fillRepArrays.R#14
1: get_fillRepArrays(stock = stock[[i]])

It's a little hard for me to figure this one out, it is happening when y=fmyearIdx maybe not.

[mdmazur]

I am not getting any errors when I run this. Is this on the EconOnly2 branch?

[mle2718]

@mdmazur. Yes, it's the EconOnly2 branch, which now contains the changes from Fix_pollock_and_codGB_problems. I will try running the Fix_pollock_and_codGB_problems branch with

stockExclude <- c('codGB', 'pollock', 'yellowtailflounderGB')

to figure out if it's a problem with the EconOnly2 branch or something else. I have a bad feeling it's a something really obscure with R.

[mle2718]

All stocks need SSB_star (which is the hinge point), cR (a conversion factor) and Rnyr (the number of years of recruitment used in the empirical distribution function). Those should be specified in the scripts in stockParameters.

@jerellejesse --this is something you can address, right?

[jerellejesse]

Yes, I will double check that they all have that set up correctly!