Commit fd722666 authored by Benoit Berges's avatar Benoit Berges
Browse files

script update correcting bug

parent 6d8ce18d
......@@ -10,6 +10,7 @@ library(FLCore)
library(stockassessment)
library(FLSAM)
library(methods)
library(icesTAF)
mkdir("model")
mkdir(file.path("model","config"))
......
......@@ -9,6 +9,7 @@ rm(list=ls())
library(FLSAM)
library(FLEDA)
library(FLFleet)
library(tidyverse)
# paths to different subfolders
dataPath <- file.path(".","data/")
......@@ -56,20 +57,22 @@ selPeriod <- ac(2010:2020)
fecYears <- ac(2010:2020)
nits <- 200 # number of random samples
load(file.path(dataPath,paste0("data_sms",SMSkeyRuns,"_sf.RData")))
# reading the raw M and applying plus group
#raw_M <- read.csv(file.path(dataPath,"Smoothed_span50_M_NotExtrapolated_NSASSMS2016.csv"),header=TRUE)
raw_M <- read.csv(file.path(dataPath,
paste0("M_NSAS_smoothedSpan50_notExtrapolated_sms",ac(SMSkeyRuns),".csv")),
header=TRUE)
colnames(raw_M) <- sub("X","",colnames(raw_M))
rownames(raw_M) <- raw_M[,1]
raw_M <- raw_M[,-1]# Trim off first column as it contains 'ages'
raw_M <- cbind(replicate(as.numeric(colnames(raw_M)[1])-histMinYr,raw_M[,1]), raw_M)
raw_M <- cbind(raw_M,raw_M[,dim(raw_M)[2]])
colnames(raw_M) <- histMinYr:histMaxYr
raw_M <- raw_M+0.06
#raw_M <- read.csv(file.path(dataPath,
# paste0("M_NSAS_smoothedSpan50_notExtrapolated_sms",ac(SMSkeyRuns),".csv")),
# header=TRUE)
M2M1_raw <- read.csv(file.path(dataPath,"SMS_NSAS_M_raw.csv"),header=TRUE,check.names = FALSE)
M2M1_raw <- M2M1_raw[M2M1_raw$Source == '2019',]
M2M1_raw <- M2M1_raw %>% select(-Source) %>% pivot_wider(names_from = year, values_from = M)
M2M1_raw <- M2M1_raw[,-1]# Trim off first column as it contains 'ages'
M2M1_raw <- cbind(replicate(as.numeric(colnames(M2M1_raw)[1])-histMinYr,M2M1_raw[,1]), M2M1_raw)
M2M1_raw <- cbind(M2M1_raw,M2M1_raw[,dim(M2M1_raw)[2]])
colnames(M2M1_raw) <- histMinYr:histMaxYr
raw_M <- M2M1_raw+0.06
#plot(colnames(raw_M),raw_M[1,],type="l")
#-------------------------------------------------------------------------------
# 3) create random samples using variance/covariance matrix
......@@ -109,6 +112,9 @@ biol <- FLStock( catch.wt=FLQuant(dim=C1,dimnames=dmns),
harvest.spwn=FLQuant(dim=C1,dimnames=dmns),
m.spwn=FLQuant(dim=C1,dimnames=dmns))
range(biol)['minfbar'] <- 2
range(biol)['maxfbar'] <- 6
units(biol) <- units(NSH)
for(idxIter in 1:nits){
......
......@@ -133,6 +133,8 @@ for (iYr in 2021:2042){
biol@harvest[,ac(iYr-1)] <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,])
#- Update recruitment
#biol.sr@rec <- rec(biol[,ac(iYr-1)])
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])
stock.n(biol)[1,ac(iYr)] <- recruitBio
......@@ -160,11 +162,11 @@ for (iYr in 2021:2042){
# filter stock object up to intermediate year to include biological variables
stkAssessment <- window(biol,end=an(TaY))
stkAssessment@catch.n <- stkAssessment@catch.n * catchVar[,ac(histMinYr:TaY),,1]
stkAssessment@catch.n <- stkAssessment@catch.n * catchVar[,ac(histMinYr:TaY),,'residuals']
stkAssessment@landings.n <- stkAssessment@catch.n
stkAssessment@landings <- computeLandings(stkAssessment)
stkAssessment <- window(biol,end=an(TaY))
#stkAssessment <- window(biol,end=an(TaY))
# smooth M prior to running the assessment, median filter of order 5
require(doParallel); ncores <- detectCores()-1; ncores <- ifelse(nits<ncores,nits,ncores);cl <- makeCluster(ncores); registerDoParallel(cl)
......@@ -191,10 +193,10 @@ for (iYr in 2021:2042){
stkAssessment.ctrl@range[5] <- an(TaY)+1
stkAssessment.ctrl@residuals <- F
stk0 <- stkAssessment
idx0 <- stkAssessment.tun
sam0.ctrl <- stkAssessment.ctrl
escapeRuns <- escapeRuns
#stk0 <- stkAssessment
#idx0 <- stkAssessment.tun
#sam0.ctrl <- stkAssessment.ctrl
#escapeRuns <- escapeRuns
#resInit <- stkAssessment.init
# Run stock assessment
......@@ -207,6 +209,14 @@ for (iYr in 2021:2042){
escapeRuns <- ret$escapeRuns
stkAssessment <- ret$stk
#compObj <- new('FLStocks')
#compObj[['assessment']] <- stkAssessment
#compObj[['biol']] <- biol
#windows()
#plot(compObj)
#print(escapeRuns)
#print(stkAssessment@stock.n[,ac(TaY)] / biol@stock.n[,ac(TaY)])
......@@ -260,6 +270,8 @@ for (iYr in 2021:2042){
lower=rep(1e-8,4))$par)
stopCluster(cl); detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE)
#catchVar[1,ac(iYr),,"FCprop",,]
#Check for very high F
idx <- which(quantMeans(sweep(landings.sel(fishery)[,ac(ImY)],3:6,t(multso),"*")[ac(2:6),,,,"A"]) > 5)
if(length(idx)>0){
......
......@@ -7,8 +7,8 @@
rm(list=ls())
HCR <- 'FMSYAR'
ftarget <- 0.33
btrigger <- 1190149
ftarget <- 0.31
btrigger <- 1232828
cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
......@@ -56,12 +56,14 @@ source(file.path(functionPath,"find.FCAR.r"))
nits <- 200
load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
#load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
# load object
load(file.path(modelPath,paste0(runName,'_MSE_init_',ac(nits),'.RData')))
stkAssessment.ctrl <- NSH.ctrl
# load MSE parameters
load(file.path(modelPath,paste0(runName,'_MSE_params_',ac(nits),'.RData')))
fishery@landings.sel[,projPeriod] <- sweep(fishery@landings.sel[,projPeriod],c(2:6),quantMeans(fishery@landings.sel[,projPeriod]),"/")
......@@ -80,13 +82,13 @@ runNameSave <- paste0(runName,"_MSE.init_",nits)
#
#-------------------------------------------------------------------------------
referencePoints <- list(Fmsy = 0.33,
referencePoints <- list(Fmsy = 0.31,
Fsq = NA,
Flim = 0.42,
Fpa = 0.33,
Blim = 833125,
Bpa = 1157692,
MSYBtrigger = 1190149,
Flim = 0.39,
Fpa = 0.31,
Blim = 874198,
Bpa = 956483,
MSYBtrigger = 1232828,
Ftarget = ftarget,
F01 = 0.05,
Btrigger = btrigger)
......
......@@ -10,7 +10,7 @@ library(FLSAM)
library(FLEDA)
library(FLFleet)
# paths to different subfolders
# paths to different subfoldersa
dataPath <- file.path(".","data/")
modelPath <- file.path(".","model/")
functionPath <- file.path(".","functions/")
......
......@@ -11,6 +11,12 @@ TAC2sel <- function(mult,iYr,iBiol,iFishery,iTAC,catchVar,TAC_var,iTer){
Ctarget <- mean(rowSums(Fs)*catchVar[1,iYr,,"FCprop",,iTer,drop=T])
Dtarget <- iTAC[,iYr,,,"D",iTer,drop=T] * TAC_var[iYr,iTer,"Dsplit"] * TAC_var[iYr,iTer,"Duptake"]
print(iBiol@stock.n[,iYr,,,,iTer,drop=T])
#print(Atarget)
#print(Btarget)
#print(Ctarget)
#print(Dtarget)
ret <- sqrt(c(Atarget - Cs[1],Btarget - Cs[2],Ctarget - mean(Fs[,3],na.rm=T),Dtarget - Cs[4])^2)
......
......@@ -5,7 +5,6 @@ config_baseline_sf <- function( NSH,
NSH.ctrl <- FLSAM.control(NSH,NSH.tun)
catchRow <- grep("catch",rownames(NSH.ctrl@f.vars))
laiRow <- grep("LAI",rownames(NSH.ctrl@obs.vars))
#Variance in N random walks (set 1st one free is usually best)
NSH.ctrl@logN.vars[] <- c(1,rep(2,length(1:pg)))
......@@ -21,13 +20,11 @@ config_baseline_sf <- function( NSH,
NSH.ctrl@obs.vars["IBTS-Q1",ac(1)] <- 201
NSH.ctrl@obs.vars["IBTS0",ac(0)] <- 301
NSH.ctrl@obs.vars["IBTS-Q3",ac(0:5)] <- c(rep(400,2),rep(401,4))#c(400,401,rep(402,4))
NSH.ctrl@obs.vars[laiRow,1] <- 501
#Catchabilities of the surveys. Set LAI all to 1 value, rest can be varied
NSH.ctrl@catchabilities["HERAS",ac(1:pg)] <- c(rep(100,2),rep(101,length(3:pg)))
NSH.ctrl@catchabilities["IBTS-Q1",ac(1)] <- c(201)
NSH.ctrl@catchabilities["IBTS-Q3",ac(0:5)] <- 300:305#c(300,rep(301,3),302,303)
NSH.ctrl@catchabilities[laiRow,1] <- 401
#Add correlation correction for Q3 survey, not for HERAS
idx <- which(rownames(NSH.ctrl@cor.obs)=="IBTS-Q3")
......@@ -36,7 +33,7 @@ config_baseline_sf <- function( NSH,
#Variance of F random walks
NSH.ctrl@f.vars[1,] <- c(101,101,rep(102,4),rep(103,length(6:pg)))
#No correlation structure among ages in F random walks
#correlation structure among ages in F random walks
NSH.ctrl@cor.F <- 2
#Finalise
NSH.ctrl@name <- "North Sea Herring"
......
......@@ -65,7 +65,7 @@ fmsyAR_fun <- function( stf,
stf@harvest[,FcY,c('A','B','C','D')] <- sweep(stf@harvest[,FcY,c('A','B','C','D')],
c(3,6),res,"*")
print(fbar(stf[,FcY,'A']))
#print(fbar(stf[,FcY,'A']))
#stf@harvest[,FcY,c('A','C','D')] <- sweep(stf@harvest[,FcY,c('A','C','D')],
# c(3,6),res,"*")
......@@ -90,6 +90,11 @@ fmsyAR_fun <- function( stf,
iYr=FcY,
CATCH=CATCH[,FcY])
#print(stf@harvest[,ImY])
#print(stf@stock.n[,ImY,,,,1,drop=T])
#print(ssb(stf[,ImY]))
# update catch and landings for each fleet in Forecast year
for(i in dms$unit){
stf@catch.n[,FcY,i] <- stf@stock.n[,FcY,i]*(1-exp(-unitSums(stf@harvest[,FcY])-stf@m[,FcY,i]))*(stf@harvest[,FcY,i]/(unitSums(stf@harvest[,FcY])+stf@m[,FcY,i]))
......
......@@ -88,6 +88,9 @@ projectNSH_FMSYAR <- function(iStocks,
stf@landings[,ImY,i] <- computeLandings(stf[,ImY,i])
}
#print(fbar(stf[,,'A']))
#print(iCATCH[,ImY,,,"A"])
#===============================================================================
# Forecast year
#===============================================================================
......@@ -111,7 +114,7 @@ projectNSH_FMSYAR <- function(iStocks,
f26)
#print(res$stf@catch[,FcY,'A'])
#print(fbar(res$stf[,FcY,'A']))
print(fbar(res$stf[,FcY,'A']))
iCATCH[,FcY,,,c("A","B","C","D")] <- res$stf@catch[,FcY]
totF <- unitSums(res$stf@harvest[,FcY])
......
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