Commit 6a9164b2 authored by Benoit Berges's avatar Benoit Berges
Browse files

perf metrics

parent a7d9fd1b
rm(list=ls())
load('model/NSAS_sanoba_200.RData')
metricsPeriod <- projPeriod[(length(projPeriod)-10):(length(projPeriod)-1)]
f26 <- ac(2:6)
f01 <- ac(0:1)
biol@catch <- computeCatch(biol)
biol@stock <- computeStock(biol)
biol@landings <- computeLandings(biol)
# risk of SSB < Blim
SSB <- ssb(biol[,metricsPeriod])
SSB <- drop(SSB)
SSB_riskMat <- array(FALSE,dim=dim(SSB))
SSB_bool <- array(FALSE,dim=c(1,nits))
for(idxIter in 1:nits){
# store value per year
SSB_riskMat[which(SSB[,idxIter] < referencePoints$Blim),idxIter] <- TRUE
# TRUE/FALSE for each iteration
if(length(which(SSB[,idxIter] < referencePoints$Blim)!=0))
SSB_bool[idxIter] <- TRUE
}
SSB_prob <- array(NA,dim=c(1,length(metricsPeriod)))
for(idxProb in 1:length(metricsPeriod)){
SSB_prob[idxProb] <- length(which(SSB_riskMat[idxProb,] == TRUE))/nits
}
SSBTail <- apply(SSB, 1, quantile, probs=c(0.05, 0.5, 0.95), na.rm=TRUE)
lmFrame <- as.data.frame(t(rbind(an(metricsPeriod),SSB_prob,SSBTail[1,])))
colnames(lmFrame) <- c('years','risk','SSB5per')
riskLR <- lm(risk~years,data=lmFrame)
SSBTailLR <- lm(SSB5per~years,data=lmFrame)
# plot
par(mfrow=c(2,1))
plot(metricsPeriod,SSBTail[2,],type='l',
ylim=c(0,max(SSBTail[3,])),
ylab='SSB',
xlab='years',main=paste0('ftar_',ftarget,'_btrig_',btrigger))
lines(metricsPeriod,SSBTail[1,],type='l',col='blue',lty=2)
lines(metricsPeriod,SSBTail[3,],type='l',col='blue',lty=2)
plot(metricsPeriod,SSB_prob,type='l',
ylim=c(0,0.2),
ylab='Annual risk',
xlab='years',
xlim=c(an(min(metricsPeriod)),an(max(metricsPeriod))))
lines(c(2010,2040),c(0.05,0.05),lty=2,col='green')
# Fbar
FHCR26 <- FHCR[f26,ac(2021:2041)]
FHCR01 <- FHCR[f01,ac(2021:2041)]
f01Mat <- apply(drop(FHCR01[,metricsPeriod]),c(2,3),'mean')
f26Mat <- apply(drop(FHCR26[,metricsPeriod]),c(2,3),'mean')
f01Quant <- apply(f01Mat, 1, quantile, probs=c(0.05, 0.5, 0.95), na.rm=TRUE)
f26Quant <- apply(f26Mat, 1, quantile, probs=c(0.05, 0.5, 0.95), na.rm=TRUE)
# LTY
catchQuant <- apply(drop(biol[,metricsPeriod]@catch), 1, quantile, probs=c(0.05, 0.5, 0.95), na.rm=TRUE)
LTR1 <- mean(SSB_prob) #length(which(SSB_bool))/nits
LTR3 <- max(SSB_prob) #length(which(SSB_bool))/nits
LTY <- mean(catchQuant['50%',])
LTF01 <- mean(f01Quant['50%',])
LTF26 <- mean(f26Quant['50%',])
LT_RiskTrend <- summary(riskLR)$coefficients[2,1]
LT_RiskTrendSig <- ifelse(summary(riskLR)$coefficients[2,4] > 0.05,0,1) # 0 is not significant. 1 is significant
LT_SSBTailTrend <- summary(SSBTailLR)$coefficients[2,1]
LT_SSBTailTrendSig <- ifelse(summary(SSBTailLR)$coefficients[2,4] > 0.05,0,1) # 0 is not significant. 1 is significant
\ No newline at end of file
......@@ -107,7 +107,7 @@ SSBHCR <- FLQuant(NA, dimnames=list(age="all",year=ac(an(proj
start.time <- Sys.time()
#an(projPeriod)
for (iYr in 2021:2042){
for (iYr in 2032:2042){
cat(iYr,"\n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment