search

lkerr / groundfish-MSE

Development of Robust Management Strategies for Northeast Groundfish Fisheries in a Changing Climate
MIT License
4 stars 2 forks source link

Error in get_fillRepArrays() #227

Closed mle2718 closed 2 years ago

mle2718 commented 2 years ago

When I run the simulation using the Fix_pollock_and_codGB_problems branch, the following lines of code:

https://github.com/lkerr/groundfish-MSE/blob/6549b3dcbb0b149aad98c61c0d8e47c6a929ff94/processes/runSim.R#L133-L138 produce an error message, when I try to run just a CodGOM model.

This occurs when nrep<-2, but not when nrep<-1

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]])

Originally posted by @mle2718 in https://github.com/lkerr/groundfish-MSE/issues/218#issuecomment-1154257836

mle2718 commented 2 years ago

A couple of thoughts.

  1. I'm running in Linux and calling ASAP3 natively (not with WINE). I'm running inside Rstudio, usually by sourcing the "runSim.R" file. But sometimes, by just running bit of code.

  2. It looks like the final ASAP model run may not be converging.
    a. I get red colored text on the final model run about the covariance matrix not being positive definite. b. I don't love this explanation for things going wrong because I also get this message in earlier years.

  3. runSim.R seems to run just fine when I set nrep<-1. However, when I set nrep<-2, this error message appears.

mle2718 commented 2 years ago

Okay, this is definitely due to the me resetting the nrep after the containers are setup.

samtruesdell commented 2 years ago

So convergence issues resolved?

mle2718 commented 2 years ago

@samtruesdell : I don't know. They were not the cause of the problem (it was my bad code).