Commit 480ff2c2 authored by Berges, Benoit's avatar Berges, Benoit
Browse files

moving IBP stuff to the attic

parent de31e9c5
......@@ -9,9 +9,11 @@ rm(list=ls())
library(FLSAM)
library(FLEDA)
library(FLFleet)
library(tidyverse)
# paths to different subfolders
dataPath <- file.path(".","data/")
resultPath <- file.path(".","results/")
modelPath <- file.path(".","model/")
functionPath <- file.path(".","functions/")
......@@ -26,7 +28,7 @@ source(file.path(functionPath,"make.arma.resid.lst.R"))
#-------------------------------------------------------------------------------
# 2) load assessment objects (single and multi fleet)
# load assessment objects (single and multi fleet)
# define MSE parameters
# load raw M
#
......@@ -37,6 +39,9 @@ source(file.path(functionPath,"make.arma.resid.lst.R"))
# complicated to implement
#-------------------------------------------------------------------------------
SST <- read.csv(file.path(resultPath,'00_NSAS_SST','SST_NSAS.csv'))
SST$sstResi <- (SST$sst-mean(SST$sst))/sd(SST$sst)
load(file.path(modelPath,
paste0(runName,"_assessment_results_sf.RData")))
load(file.path(modelPath,
......@@ -72,7 +77,7 @@ colnames(raw_M) <- histMinYr:histMaxYr
raw_M <- raw_M+0.06
#-------------------------------------------------------------------------------
# 3) create random samples using variance/covariance matrix
# create random samples using variance/covariance matrix
# initialize biol object
#
# Note: this takes forever, therefore saving the object after this process
......@@ -120,7 +125,7 @@ biol <- window(window(biol,end=histMaxYr+1),start=histMinYr,end=futureM
biol@m.spwn[,ac((histMaxYr+1):futureMaxYr)] <- 0.67
#-------------------------------------------------------------------------------
# 4) create FLStocks object using random samples (with future years as NA)
# create FLStocks object using random samples (with future years as NA)
#-------------------------------------------------------------------------------
biol@harvest.spwn[,projPeriod] <- biol@harvest.spwn[,ac(histMaxYr)] # propagate Fprop before spawning
......@@ -128,7 +133,7 @@ biol@m.spwn[,projPeriod] <- biol@m.spwn[,ac(histMaxYr)] # propagate Fprop
biol@stock <- computeStock(biol)
#-------------------------------------------------------------------------------
# 5) allocating future maturity, stock weight and M at age
# allocating future maturity, stock weight and M at age
#
# w@a and mat are on the same randomization
# M gets an independent randomization
......@@ -231,7 +236,7 @@ biol@landings.wt <- biol@catch.wt
biol@m[,histPeriod] <- as.matrix(raw_M)
#-------------------------------------------------------------------------------
# 6) creating survey indices
# creating survey indices
#
# Filling matrix of catchabilities and residuals for the different surveys
#
......@@ -303,7 +308,7 @@ for(idxIter in 1:nits){
}
#-------------------------------------------------------------------------------
# 7) creating catches from each random samples
# creating catches from each random samples
#
# creating residuals and Fc proportions for future years
#
......@@ -348,7 +353,7 @@ for(idxIter in 1:nits){
biol@catch <- computeCatch(biol)
#-------------------------------------------------------------------------------
# 8) C fleet: proportion of F of the C fleet for the future years
# C fleet: proportion of F of the C fleet for the future years
# The proportion (all ages combined) is obtained from the multi-fleet assessment
# using a random draw similar to M and weight at age
# For the given quantity, we randomize the following:
......@@ -407,7 +412,7 @@ for(idxIter in 1:nits){
}
#-------------------------------------------------------------------------------
# 9) create selection patterns for the different fleets
# create selection patterns for the different fleets
# Selectivity of fleet projected forward using a random walk, using results from
# the multi-fleet assessment
#-------------------------------------------------------------------------------
......@@ -490,50 +495,41 @@ for(idxFleet in 1:length(fleets)){
}
#-------------------------------------------------------------------------------
# 10) Future recruitment
# surface temperature residuals
#-------------------------------------------------------------------------------
recPeriod <- ac(2002:2021)
biol.sr <- fmle(as.FLSR(biol,model='ricker')) # just to populate the structure
params(biol.sr) <- params(fmle(FLSR(rec = rec(biol)[,ac(an(recPeriod)[2]:max(an(recPeriod)))], # rec from 1948 to 2017
ssb = ssb(biol)[,ac(an(recPeriod)[1]:(max(an(recPeriod))-1))], # ssb from 1947 to 2016
model='ricker')))
SST_resi <- as.data.frame(matrix(rlnorm( nFutureyrs*nits,
0,
sd(SST$sstResi)),
nrow=nits,byrow=TRUE))
colnames(SST_resi) <- projPeriod
rownames(SST_resi) <- 1:nits
# THIS IS MODIFIED SO THAT AN ARIMA MODEL IS FITTED FOR THE RESIDUALS OF EACH ITERATION
# AND USE TO PRODUCE THE FUTURE DEVIATIONS FOR THE CORRESPONDING ITERATION
# I want to add something that takes autocorrelation in SR relationship into account
# to do this I use an arima model
### S/R residuals - with autocorrelation
rec.res <- residuals(biol.sr)[,ac(an(recPeriod)[2]:(max(an(recPeriod)))-1)]
#-------------------------------------------------------------------------------
# Future recruitment
#-------------------------------------------------------------------------------
# autoregressive model order 1
set.seed(108)
biol.sr <- fmle(as.FLSR(biol[,histPeriod],model='ricker')) # just to populate the structure
# a list with one model per iteration
#params(biol.sr) <- params(fmle(FLSR(rec = rec(biol)[,ac(an(recPeriod)[2]:max(an(recPeriod)))], # rec from 1948 to 2017
# ssb = ssb(biol)[,ac(an(recPeriod)[1]:(max(an(recPeriod))-1))], # ssb from 1947 to 2016
# model='ricker')))
arima.fit.lst <- list()
for(its in 1:dims(biol)$iter)
arima.fit.lst[[its]] <- try(arima(an(iter(rec.res,its)),order=c(1,0,0)))
idx <- which(unlist(lapply(arima.fit.lst,function(x){class(x)=="try-error"}))==T)
for(its in idx)
arima.fit.lst[[its]] <- try(arima(an(iter(rec.res,its))))
table(unlist(lapply(arima.fit.lst,class)))
#recPeriod <- ac(2002:2021)
sr.res <- quantSums(FLQuant(NA,dimnames=dmns))
#ny <- 20 # number of years to project - Usually 20
#dy <- range(stkMC)["maxyear"] # data year
#ay <- dy # assessment year
#iy <- ay+1 # initial projections year (also intermediate)
#fy <- iy + ny -1 # final year
# generate random blocks for residuals
yrChain <- randBlocks(an(fecYears),an(projPeriod),nits)
# create autocorrelation in residuals and propagate throughout stock into the future
# from initial year of projections (iy) to last of projections (ny-1)
sr.res <- make.arma.resid.lst(arima.fit.lst, age = 0, years = an(projPeriod[1]):max(an(projPeriod)), rec.res)
for(idxIter in 1:nits){
print(paste('init step 11 bio variables - iter=',idxIter))
iter(sr.res[,projPeriod],idxIter) <- iter(residuals(rec(biol.sr), predict(biol.sr))[,ac(yrChain[[1]])],idxIter)
}
#-------------------------------------------------------------------------------
# 11) process error
# process error
#-------------------------------------------------------------------------------
projYearsCohort <- (an(projPeriod)[1]-8):(max(an(projPeriod)))
......@@ -565,7 +561,7 @@ procError <- surv[ac(1:8),histPeriod[1:(length(histPeriod)-1)]]/ # survivors ag
for(idxIter in 1:nits){
print(paste('init step 11 process error - iter=',idxIter))
print(paste('init step 12 process error - iter=',idxIter))
# covariance accross the ages using a 10 year period of full cohorts
covMat <- cov(t(FLCohort(log(procError))[,ac(1999:2008),,,,idxIter,drop=T]))
......@@ -576,7 +572,7 @@ for(idxIter in 1:nits){
}
#------------------------------------------------------------------------------#
# 12) Define TACs for A, B and D fleets.
# Define TACs for A, B and D fleets.
# TACs for A and B fleets are taken out of HAWG2018. This needs updating
#
# Note 1: The Cfleet is defined as a proportion of F.
......@@ -644,7 +640,7 @@ TAC_var[,,'Csplit'] <- t(Csplit)
#-------------------------------------------------------------------------------
# 13) tidying up and saving objects for next step
# tidying up and saving objects for next step
#-------------------------------------------------------------------------------
# prepare stock object
......@@ -654,7 +650,7 @@ units(biol) <- units(NSH)
recFuture <- array(NA,dim=c(nits,1))
for(idxIter in 1:nits){
# recruitment out of SAM for 2018
recFuture[idxIter] <- as.array(subset(rec(NSH.sim[[idxIter]]),year==2018)$value)
recFuture[idxIter] <- as.array(subset(rec(NSH.sim[[idxIter]]),year==2021)$value)
}
......@@ -662,11 +658,9 @@ for(idxIter in 1:nits){
save( biol, # biology object
varProccError, # future process error (across cohorts)
catchVar, # future observation variance for the catches
recFuture, # recruitment for 2018 for stf from SAM object
recFuture, # recruitment for 2021 for stf from SAM object
sr.res, # residuals for stock-recruitment
biol.sr, # stock-recruitment fits
itersSR, # iterations for segmented regression
itersRI, # iterations for Ricket
fishery, # fishery object (fleet wise), contains selection patterns + catch.wt
surveys, # survey object
surveyVars, # future catchabilities and residuals for the surveys
......
......@@ -54,7 +54,7 @@ source(file.path(functionPath,"find.FCAR.r"))
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 200
nits <- 10
#load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
......@@ -135,24 +135,9 @@ for (iYr in 2021:2042){
biol@harvest[,ac(iYr-1)] <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,])
#- Update recruitment
recruitBio <- array( 0, dim=c(1,nits)) # initialize array
biol.sr@ssb[,ac(iYr-1)] <- ssb(biol[,ac(iYr-1)])
recruitBio <- predict(biol.sr[,ac(iYr-1)]) * exp(sr.res[,ac(iYr),drop=T])
paramRec <- params(biol.sr)
ssbRec <- drop(ssb(biol[,ac(iYr-1)]))
# Ricker
recruitBio[itersRI] <- paramRec['a',itersRI]*ssbRec[itersRI]*exp(-paramRec['b',itersRI]*ssbRec[itersRI])
# Segmented regression
idxSSR1 <- which(ssbRec[itersSR] <= paramRec['b',itersSR])
# if loop to cover the case of idxSSR1 being empty
if(length(idxSSR1)!=0){
recruitBio[itersSR[idxSSR1]] <- paramRec['a',itersSR[idxSSR1]]*ssbRec[itersSR[idxSSR1]] # SSB < b (slope)
recruitBio[itersSR[-idxSSR1]] <- paramRec['a',itersSR[-idxSSR1]]*paramRec['b',itersSR[-idxSSR1]] # SSB > b (plateau)
}else{
recruitBio[itersSR[idxSSR1]] <- paramRec['a',itersSR]*paramRec['b',itersSR]
}
recruitBio <- recruitBio * exp(sr.res[,ac(iYr),drop=T])
stock.n(biol)[1,ac(iYr)] <- recruitBio
#- Plusgroup
......@@ -214,7 +199,7 @@ for (iYr in 2021:2042){
idx0 <- stkAssessment.tun
sam0.ctrl <- stkAssessment.ctrl
escapeRuns <- escapeRuns
resInit <- stkAssessment.init
#resInit <- stkAssessment.init
# Run stock assessment
ret <- MSE_assessment( stkAssessment,
......
......@@ -19,10 +19,10 @@ biol@stock <- computeStock(biol)
biol@landings <- computeLandings(biol)
compObj[['no_transfer']] <- biol
plot(window(compObj,start=2000))
#plot(window(compObj,start=2000))
#metricsPeriod <- projPeriod[(length(projPeriod)-10):(length(projPeriod)-1)]
metricsPeriod <- ac(2031:2039)
metricsPeriod <- ac(2031:2042)
f26 <- ac(2:6)
f01 <- ac(0:1)
......
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