Commit 207df5d0 authored by bberges's avatar bberges
Browse files

initial commit

parents
#-------------------------------------------------------------------------------
# WKNSMSE
#
# Author: Benoit Berges
# WMR, The Netherland
# email: benoit.berges@wur.nl
#
# script to test the effect of dropping LAI indices
#
# Date: 2018/12/11
#
# Build for R3.5.1, 64bits
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Running assessment
#-------------------------------------------------------------------------------
rm(list=ls())
library(FLSAM)
library(FLEDA)
#library(R.utils)
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- '/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/'
#path <- 'E:/wk_WKNSMSE_her.27.3a47d/R'
assessment_name_noLAI <- "NSAS_WKNSMSE2018"
try(setwd(path),silent=TRUE)
# paths to different subfolders
dataPath <- file.path(".","data/")
outPath <- file.path(".","results/")
scriptPath <- file.path(".","side_scripts/")
# running assessment with no LAI
source(file.path(scriptPath,"setupAssessmentObjects_sf_noLAI.R"))
source(file.path(scriptPath,"setupControlObject_sf_noLAI.R"))
NSHnoLAI.ctrl <- NSH.ctrl
NSHnoLAI.tun <- NSH.tun
NSHnoLAI <- NSH
NSHnoLAI.sam <- FLSAM(NSHnoLAI,
NSHnoLAI.tun,
NSHnoLAI.ctrl)
NSHnoLAI@stock.n <- NSHnoLAI.sam@stock.n[,ac(range(NSH)["minyear"]:range(NSH)["maxyear"])]
NSHnoLAI@harvest <- NSHnoLAI.sam@harvest[,ac(range(NSH)["minyear"]:range(NSH)["maxyear"])]
NSHnoLAI.tun@residuals <- F
n.retro.years <- 7
NSHnoLAI.retro <- retro( NSHnoLAI,
NSHnoLAI.tun,
NSHnoLAI.ctrl,
retro=n.retro.years)
#-------------------------------------------------------------------
# saving
#-------------------------------------------------------------------
save(NSHnoLAI,
NSHnoLAI.tun,
NSHnoLAI.ctrl,
NSHnoLAI.sam,
NSHnoLAI.retro,
file=file.path(outPath,paste0(assessment_name_noLAI,'_sf_noLAI_retro.Rdata')))
# modify object names so it fits those in the other scripts
NSH <- NSHnoLAI
NSH.tun <- NSHnoLAI.tun
NSH.ctrl <- NSHnoLAI.ctrl
NSH.sam <- NSHnoLAI.sam
save(NSH,
NSHnoLAI.tun,
NSHnoLAI.ctrl,
NSHnoLAI.sam,
file=file.path(outPath,paste0(assessment_name_noLAI,'_sf.Rdata')))
#-------------------------------------------------------------------
# loading assessments with and without LAI index
#-------------------------------------------------------------------
outputName <- 'LAI_comparison'
PDF <- FALSE
PNG <- ifelse(PDF,F,T)
if(PDF) pdf(file.path(outPath,'plots',paste0(outputName,".pdf")))
if(PNG) png(file.path(outPath,'plots',paste0(outputName,"_%02d.png")),
units = "px",
height=800,
width=672,
bg = "white")
load(file.path(outPath,paste0('NSAS_WKNSMSE2018_sf_retro.Rdata')))
load(file.path(outPath,paste0('NSAS_WKNSMSE2018_sf_noLAI_retro.Rdata')))
#-------------------------------------------------------------------
# Compute mohn's rhos
#-------------------------------------------------------------------
# SSB
SSB_mrNoLAI <- mohns.rho(NSHnoLAI.retro,span=n.retro.years,ref.year=2017,type="ssb")
mean(SSB_mrNoLAI$rho)
SSB_mrLAI <- mohns.rho(NSH.retro,span=n.retro.years,ref.year=2017,type="ssb")
mean(SSB_mrLAI$rho)
# fbar
fbar_mrNoLAI <- mohns.rho(NSHnoLAI.retro,span=n.retro.years,ref.year=2017,type="fbar")
mean(fbar_mrNoLAI$rho)
fbar_mrLAI <- mohns.rho(NSH.retro,span=n.retro.years,ref.year=2017,type="fbar")
mean(fbar_mrLAI$rho)
# recruitment
rec_mrNoLAI <- mohns.rho(NSHnoLAI.retro,span=n.retro.years,ref.year=2017,type="rec")
mean(rec_mrNoLAI$rho)
rec_mrLAI <- mohns.rho(NSH.retro,span=n.retro.years,ref.year=2017,type="rec")
mean(rec_mrLAI$rho)
#-------------------------------------------------------------------
# plotting
#-------------------------------------------------------------------
ssbNoLAI <- ssb(NSHnoLAI.sam)
recruitNoLAI <- rec(NSHnoLAI.sam)
fbarNoLAI <- fbar(NSHnoLAI.sam)
ssbLAI <- ssb(NSH.sam)
recruitLAI <- rec(NSH.sam)
fbarLAI <- fbar(NSH.sam)
yrs <- ssbLAI$year
par(mfrow=c(3,1))
# SSB
# with LAI
plot(yrs,
ssbLAI$value,
xlab='year',ylab='SSB',main='SSB',type='l',col=rgb(0,0,1,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(ssbLAI$lbnd,rev(ssbLAI$ubnd)),
col=rgb(0,0,1,0.1),lty=0)
# without LAI
lines(yrs,
ssbNoLAI$value,
col=rgb(1,0,0,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(ssbNoLAI$lbnd,rev(ssbNoLAI$ubnd)),
col=rgb(1,0,0,0.1),lty=0)
legend("topright",
c("with LAI","without LAI"),
col=c(rgb(0,0,1,1),rgb(1,0,0,1)),lty=1,lwd=3)
# Fbar
# with LAI
plot(yrs,
fbarLAI$value,
xlab='year',ylab='Fbar',main='Fbar',type='l',col=rgb(0,0,1,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(fbarLAI$lbnd,rev(fbarLAI$ubnd)),
col=rgb(0,0,1,0.1),lty=0)
# without LAI
lines(yrs,
fbarNoLAI$value,
col=rgb(1,0,0,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(fbarNoLAI$lbnd,rev(fbarNoLAI$ubnd)),
col=rgb(1,0,0,0.1),lty=0)
# Recruitment
# with LAI
plot(yrs,
recruitLAI$value,
xlab='year',ylab='rec',main='rec',type='l',col=rgb(0,0,1,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(recruitLAI$lbnd,rev(recruitLAI$ubnd)),
col=rgb(0,0,1,0.1),lty=0)
# without LAI
lines(yrs,
recruitNoLAI$value,
col=rgb(1,0,0,1),lwd=3)
polygon(c(yrs,rev(yrs)),
c(recruitNoLAI$lbnd,rev(recruitNoLAI$ubnd)),
col=rgb(1,0,0,0.1),lty=0)
# retrospective comparison
yrs_retro <- SSB_mrNoLAI$year
par(mfrow=c(3,1))
# SSB
plot(yrs_retro,
SSB_mrLAI$rho,
xlab='year',ylab='mohn rho SSB',main='mohn rho SSB',type='l',col=rgb(0,0,1,1),lwd=3)
lines(yrs_retro,
SSB_mrNoLAI$rho,
col=rgb(1,0,0,1),lwd=3)
legend("topright",
c("with LAI","without LAI"),
col=c(rgb(0,0,1,1),rgb(1,0,0,1)),lty=1,lwd=3)
# fbar
plot(yrs_retro,
fbar_mrLAI$rho,
xlab='year',ylab='mohn rho fbar',main='mohn rho fbar',type='l',col=rgb(0,0,1,1),lwd=3)
lines(yrs_retro,
fbar_mrNoLAI$rho,
col=rgb(1,0,0,1),lwd=3)
# recruitment
plot(yrs_retro,
rec_mrLAI$rho,
xlab='year',ylab='mohn rho rec',main='mohn rho rec',type='l',col=rgb(0,0,1,1),lwd=3)
lines(yrs_retro,
rec_mrNoLAI$rho,
col=rgb(1,0,0,1),lwd=3)
dev.off()
library(FLCore)
rm(list=ls())
path <- "J:/git/sanoba_nsas"
try(setwd(path),silent=TRUE)
################################################
# setup paths
################################################
data.dir <- file.path(".","data")
script.dir <- file.path(".","side_scripts")
function.dir <- file.path(".","functions")
results.dir <- file.path(".","results")
figures.dir <- file.path(".","figures")
runName <- '01a_growth_VB_fit'
dir.create(file.path(results.dir,runName),showWarnings = FALSE)
dir.create(file.path(figures.dir,runName),showWarnings = FALSE)
results.dir <- file.path(".","results",runName)
figures.dir <- file.path(".","figures",runName)
################################################
# read data
################################################
NSH <- readFLStock(file.path(data.dir,"index.txt"),no.discards=TRUE,quiet=FALSE)
stk.wt <- window(NSH@stock.wt, start=1984)
window(NSH)
base<-(read.csv("baseW.csv"))
base<-base[-1470,]
names(base)[3]<-"w"
SST<-read.csv("SST.csv")
surf<-read.csv("area.csv")
#par(mfrow=c(3,4))
################################################
# compute growth curve per cohort
################################################
stocknonbalt<-c(2,4:6,8:14,16:18,20)
results<-data.frame()
for (i in stocknonbalt)
{
###########
#calcule la courbe de croissance du stock
print(i)
stock<-levels(base$stock)[i]
growth.fit<- nls(w~winf*(1-exp(-k*(age-t0)))^3,data=base[base$stock==stock,],start=list(winf=0.5,k=0.3,t0=-0.5))
winfest<-summary(growth.fit)$coefficients[1,1]
kest<-summary(growth.fit)$coefficients[2,1]
t0est<-summary(growth.fit)$coefficients[3,1]
YC<-as.double(levels(as.factor(base$YC[base$stock==stock])))
t<-seq(from=0,to=max(base$age[base$stock==stock]),by=0.1)
#windows()
plot(t,winfest*(1-exp(-kest*(t-t0est)))^3,type="l",ylim=c(0,max(base$w[base$stock==stock])*1.1),main= paste(stock," nls"),xlab="age",ylab="weight")
points(base$age[base$stock==stock],base$w[base$stock==stock])
#############
#########################""""
#cree les tableau pour stocker les rsultats
paraestnls<-data.frame(stock,YC,rep(NA,length(YC)),rep(NA,length(YC)),rep(NA,length(YC)),rep(NA,length(YC)),rep(NA,length(YC)),rep(NA,length(YC)),rep(NA,length(YC)))
SS<-rep(0,length(YC))
names(paraestnls)<-c("stock","YC","winf","k","t0","SST","dens","nb_years_W","nb_years_SST")
#################################
# ajuste le model pour chaque cohorte
for (j in 1:length(YC))
{
years<-base$year[base$stock==stock & base$YC==YC[j] ]
paraestnls$nb_years_W[j]<-length(years)
paraestnls$nb_years_SST[j]<-length(SST$Temp[SST$stock==stock & is.element(SST$year,years)])
paraestnls$SST[j]<- mean(SST$Temp[SST$stock==stock & is.element(SST$year,years)])
try(paraestnls$dens[j]<-mean (na.omit(base$SSB[base$stock==stock & is.element(base$year,years)]))/ surf$area[surf$stock==stock],silent=T)
paraestnls$nb_years_SST[j]<-length(base$SSB[base$YC==YC[j] & base$stock==stock & is.element(base$year,years)])
if (length(years)>(max( base$age[base$stock==stock])-min( base$age[base$stock==stock])-1))
{
try(
{
growth.fitI<- nls(w~winf*(1-exp(-k*(age-t0)))^3,data=base[base$stock==stock & base$YC==YC[j] ,],start=list(winf=0.4,k=0.3,t0=-0.5))
#growth.fitI<- nls(w~winf*(1-exp(-k*(age+exp(lnt0))))^3,data=base[base$stock==stock & base$YC==YC[j] ,])
SS[j]<-sum((residuals(growth.fitI))^2)
winfestI<-summary(growth.fitI)$coefficients[1,1]
kestI<-summary(growth.fitI)$coefficients[2,1]
t0estI<-(summary(growth.fitI)$coefficients[3,1])
lines(t,winfestI*(1-exp(-kestI*(t-t0estI)))^3)
paraestnls$winf[j]<-winfestI
paraestnls$k[j]<-kestI
paraestnls$t0[j]<-t0estI
rm(winfestI,kestI,t0estI)
} )
}
}
SSnls <-sum(SS)
results<-rbind(results,paraestnls)
resultssave<-results
}
\ No newline at end of file
library(FLCore)
library(TMB)
rm(list=ls())
path <- "J:/git/sanoba_nsas"
try(setwd(path),silent=TRUE)
################################################
# setup paths
################################################
data.dir <- file.path(".","data")
script.dir <- file.path(".","side_scripts")
function.dir <- file.path(".","functions")
results.dir <- file.path(".","results")
figures.dir <- file.path(".","figures")
runName <- '01b_SRR_fit'
dir.create(file.path(results.dir,runName),showWarnings = FALSE)
dir.create(file.path(figures.dir,runName),showWarnings = FALSE)
results.dir <- file.path(".","results",runName)
figures.dir <- file.path(".","figures",runName)
################################################
# read data
################################################
dat <- read.csv(file.path(results.dir,'..','00_NSAS_SST','SST_NSAS.csv'))
SSB_pred <- seq(from=min(dat$ssb),to=max(dat$ssb),by=1000)
SST_pred <- min(dat$sst)#seq(from=min(dat$sst),to=max(dat$sst),by=0.5)
################################################
# compile TMB function
################################################
compile(file.path(function.dir,"C++","ricker_T.cpp"))
dyn.load(dynlib(file.path(function.dir,"C++","ricker_T")))
data <- list(SSB=dat$ssb,
logR=log(dat$rec),
SST=dat$sst,
SSB_pred=SSB_pred,
SST_pred=SST_pred)
parameters <- list(logA=0,
logB=0,
logC=0,
logSigma=0)
obj <- MakeADFun(data,parameters,DLL="ricker_T")
opt <- nlminb(obj$par,obj$fn,obj$gr)
rep <- sdreport(obj)
logR <- rep$par.fixed[1]+log(SSB_pred)-exp(rep$par.fixed[2])*SSB_pred+exp(rep$par.fixed[3])*SST_pred
plot(dat$ssb,log(dat$rec), main = "Ricker", cex.main = 2,cex.lab = 1.5,
xlab = "SSB",ylab = "R")
lines(SSB_pred, logR, lwd=3)
lines(SSB_pred, logR+2*rep$sd, lty="dotted")
lines(SSB_pred, logR-2*rep$sd, lty="dotted")
\ No newline at end of file
#-------------------------------------------------------------------------------
# WKNSMSE
#
# Author: Niels Hintzen
# IMARES, The Netherland
#
# MSE of North Sea Herring
#
# Date: 2018/11/18
#
# Build for R3.4.3, 64bits
#-------------------------------------------------------------------------------
rm(list=ls())
library(FLSAM)
library(FLEDA)
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
path <- 'E:/wk_WKNSMSE_her.27.3a47d/R'
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- "/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/"
assessment_name <- "NSAS_WKNSMSE2018"
try(setwd(path),silent=TRUE)
# paths to different subfolders
dataPath <- file.path(".","data/")
outPath <- file.path(".","results/")
scriptPath <- file.path(".","side_scripts/")
#-------------------------------------------------------------------------------
# 0): multi-fleet assessment
#-------------------------------------------------------------------------------
source(file.path(scriptPath,"setupAssessmentObjects_mf_noLAI.r"))
source(file.path(scriptPath,"setupControlObject_mf_noLAI.r"))
NSH3f.sam <- FLSAM(NSHs3,
NSH.tun,
NSH3.ctrl)
# save assessment object
save(NSHs3,
NSH.tun,
NSH3.ctrl,
NSH3f.sam,
file=file.path(outPath,paste0(assessment_name,'_mf_noLAI.Rdata')))
#-------------------------------------------------------------------------------
# 0): Single fleet assessment
#-------------------------------------------------------------------------------
source(file.path(scriptPath,"setupAssessmentObjects_sf_noLAI.r"))
source(file.path(scriptPath,"setupControlObject_sf_noLAI.r"))
#- Perform the assessment
NSH.sam <- FLSAM(NSH,
NSH.tun,
NSH.ctrl)
name(NSH.sam) <- assessment_name
NSH@stock.n <- NSH.sam@stock.n[,ac(range(NSH)["minyear"]:range(NSH)["maxyear"])]
NSH@harvest <- NSH.sam@harvest[,ac(range(NSH)["minyear"]:range(NSH)["maxyear"])]
# save assessment object
save(NSH,
NSH.tun,
NSH.ctrl,
NSH.sam,
file=file.path(outPath,paste0(assessment_name,'_sf_noLAI.Rdata')))
This diff is collapsed.
#-------------------------------------------------------------------------------
# WKNSMSE
#
# Author: Benoit Berges
# WMR, The Netherland
# email: benoit.berges@wur.nl
#
# MSE of North Sea Herring
#
# Date: 2018/11/18
#
# Build for R3.5.1, 64bits
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 1) load packages
# setup paths
# load functions
#-------------------------------------------------------------------------------
#rm(list=ls())
args=(commandArgs(TRUE))
#args <- 'ftar=0.261_btrig=1.4e6_HCR=1_IAV=0_BB=0'
args <- strsplit(args,"_")
#ftarget <- as.numeric(substr(args[[1]][1],6,9))
#btrigger<- as.numeric(substr(args[[1]][2],7,11))
ftarget <- as.numeric(strsplit(args[[1]][1],'=')[[1]][2])
btrigger <- as.numeric(strsplit(args[[1]][2],'=')[[1]][2])
# HCR
HCR <- ifelse(as.numeric(substr(args[[1]][3],5,nchar(args[[1]][3])))==1,"A","B")
# IAV
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 0){
IAV <- 0
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 1){
IAV <- 'A'
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 2){
IAV <- 'B'
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 3){
IAV <- 'E'
}
# BB
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][4]))) == 0){
BB <- 0
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 1){
BB <- 'A'
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 2){
BB <- 'B'
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 3){
BB <- 'E'
}
if(IAV == 0){
IAV <- NULL
}else if(IAV == "B"){
IAV <- c("A","B")
}
if(BB == 0){
BB <- NULL
}else if(BB == "B"){
BB <- c("A","B")
}
cat(ftarget,"\t",btrigger,"\t",HCR,"\t",IAV,"\t",BB,"\n")
library(FLSAM)
library(FLEDA)
library(FLFleet)
library(minpack.lm) # install.packages("minpack.lm")
library(stats)
# define path to directory
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- 'E:/wk_WKNSMSE_her.27.3a47d/R'
#path <- 'D:/Repository/NSAS_MSE/wk_WKNSMSE_her.27.3a47d/R/'
#path <- "/home/hintz001/wk_WKNSMSE_her.27.3a47d/R"
path <- '/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/'
#path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
assessment_name <- "NSAS_WKNSMSE2018"
try(setwd(path),silent=TRUE)
# paths to different subfolders
dataPath <- file.path(".","data/")
outPath <- file.path(".","results/")
scriptPath <- file.path(".","side_scripts/")
functionPath <- file.path(".","functions/")
source(file.path(functionPath,"MSE_assessment.R"))
source(file.path(functionPath,"forecastScenarios.r"))
source(file.path(functionPath,"forecastFunctions.r"))
#-------------------------------------------------------------------------------
# 2) Initialize
#
# Define MSE parameters,
# load objects initialized previously
# * Biology
# - biol
# - surveys
# * Fisheries
# - FCProp
# - catchD
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 1000
# load object
load(file.path(outPath,paste0(assessment_name,'_init_MSE_',ac(nits),'.RData')))
stkAssessment.ctrl <- NSH.ctrl
# load initialization for SAM