Commit 3eb93d46 authored by Berges, Benoit's avatar Berges, Benoit
Browse files

commit changes to scripts and data

parent 207df5d0
## Preprocess data, write TAF data tables
## Part 1: Construct FLR objects
## Before: canum.txt, caton.txt, fleet.txt, fprop.txt, index.txt, lai.txt,
## matprop.txt, mprop.txt,
## Smoothed_span50_M_NotExtrapolated_NSASSMS2016.csv, weca.txt,
## west_raw.txt (bootstrap/data)
## canum
## After: data_sf.RData (data)
### ============================================================================
### init
### ============================================================================
rm(list=ls())
library(icesTAF)
library(tidyr)
library(FLCore)
library(methods)
data.source <- file.path(".","data")
# General
assessment_year <- "2021"
assessment_name <- paste0('NSAS_IBP', assessment_year)
### ============================================================================
### Prepare stock object for single fleet assessment
### ============================================================================
M2M1_raw <- read.csv(file.path(data.source,"SMS_NSAS_M_raw.csv"),header=TRUE,check.names = FALSE)
### ============================================================================
### process time series
### ============================================================================
keyRuns <- unique(M2M1_raw$Source)
for(iSMS in keyRuns){
M2M1_raw_filt <- subset(M2M1_raw,Source == iSMS & year >= 1974)
agesUnique <- unique(M2M1_raw_filt$age)
storeSmooth <- array(NA,dim=c(length(sort(unique(M2M1_raw_filt$age))),
length(sort(unique(M2M1_raw_filt$year))),
3),
dimnames=list(age=sort(unique(M2M1_raw_filt$age)),
year=sort(unique(M2M1_raw_filt$year)),
Fit=c("5%","50%","95")))
# loop on ages
for(iAge in agesUnique){
res <- predict(loess((M)~year,data=subset(M2M1_raw_filt, age == iAge),span=0.5),
newdata=expand.grid(year=sort(unique(subset(M2M1_raw_filt,age == iAge)$year))),
se=T)
yrs <- sort(unique(subset(M2M1_raw_filt, age == iAge)$year))
storeSmooth[ac(iAge),ac(yrs),] <- matrix(c(res$fit-1.96*res$se.fit,res$fit,res$fit+1.96*res$se.fit),nrow=length(yrs),ncol=3)
}
if(max(agesUnique) < 8){
res <- rbind(storeSmooth[,,"50%"],storeSmooth[8,,"50%"])
rownames(res) <- c(agesUnique,8)
}else{
res <- storeSmooth[,,"50%"]
}
write.csv(res,
file=file.path(data.source,
paste0("M_NSAS_smoothedSpan50_notExtrapolated_sms",ac(iSMS),".csv")))
}
## Preprocess data, write TAF data tables
## Part 1: Construct FLR objects
## Before: canum.txt, caton.txt, fleet.txt, fprop.txt, index.txt, lai.txt,
## matprop.txt, mprop.txt,
## Smoothed_span50_M_NotExtrapolated_NSASSMS2016.csv, weca.txt,
## west_raw.txt (bootstrap/data)
## canum
## After: data_sf.RData (data)
### ============================================================================
### init
### ============================================================================
rm(list=ls())
library(FLCore)
# sessionInfo()
# taf.session()
library(methods)
source(file.path("functions","utilities_data.R"))
mkdir("data")
data.source <- file.path(".","data")
# General
SMSkeyRunsMat <- c('2010','2013','2016','2019')
assessment_year <- "2021"
assessment_name <- paste0('NSAS_IBP', assessment_year)
for(SMSkeyRuns in SMSkeyRunsMat){
print(paste("SMS key run:", SMSkeyRuns ))
### ============================================================================
### Prepare stock object for single fleet assessment
### ============================================================================
## Load object
#NSH <- readFLStock("index.txt", no.discards=TRUE)
NSH <- readFLStock(file.path(data.source,"index.txt"),no.discards=TRUE,quiet=FALSE)
#Catch is calculated from: catch.wt * catch.n, however, the reported landings are
#normally different (due to SoP corrections). Hence we overwrite the calculate landings
# are we not using catches data then?
NSH@catch <- NSH@landings
units(NSH)[1:17] <- as.list(c(rep(c("tonnes","thousands","kg"),4),
rep("NA",2),"f",rep("NA",2)))
#Set object details
NSH@name <- assessment_name
range(NSH)[c("minfbar","maxfbar")] <- c(2,6)
NSH <- setPlusGroup(NSH,NSH@range["max"])
#Historical data is only provided for ages 0-8 prior to 1960 (rather than 0-9 after)
#We therefore need to fill in the last ages by applying the following assumptions
# - weight in the catch and weight in the stock at age 9 is the same as
# in period 1960-1983 (constant in both cases)
# - catch at age reported as 8+ is split evenly between age 8 and 9
# - natural mortality at age 9 is 0.1 (same as age 8)
# - proportion mature at age 9 is 1.0 (same as age 8)
# - harvest.spwn and m.spwn are the same as elsewhere
hist.yrs <- as.character(1947:1959)
NSH@catch.wt <- NSH@landings.wt #Automatic population of catch.wt introduces NAS
NSH@catch.wt["9",hist.yrs] <- 0.271
NSH@landings.wt["9",hist.yrs] <- 0.271
NSH@catch.n["9",hist.yrs] <- NSH@catch.n["8",hist.yrs]/2
NSH@catch.n["8",hist.yrs] <- NSH@catch.n["9",hist.yrs]
NSH@landings.n["9",hist.yrs] <- NSH@landings.n["8",hist.yrs]/2
NSH@landings.n["8",hist.yrs] <- NSH@landings.n["9",hist.yrs]
NSH@stock.wt["9",hist.yrs] <- 0.312
NSH@m["9",hist.yrs] <- 0.1
NSH@mat["9",hist.yrs] <- 1
#No catches of age 9 in 1977 so stock.wt does not get filled there.
#Hence, we copy the stock weight for that age from the previous year.
#Note that because we use a fixed stock.wt prior to 1983, there is no
#need to use averaging or anything fancier
NSH@stock.wt["9","1977"] <- NSH@stock.wt["9","1976"]
#Use a running mean(y-2,y-1,y) of input wests (i.e. west_raw) to calculate west
NSH@stock.wt[,3:dim(NSH@stock.wt)[2]] <- (NSH@stock.wt[,3:(dim(NSH@stock.wt)[2]-0)] +
NSH@stock.wt[,2:(dim(NSH@stock.wt)[2]-1)] +
NSH@stock.wt[,1:(dim(NSH@stock.wt)[2]-2)]) / 3
### ============================================================================
### Prepare natural mortality estimates
### ============================================================================
## Read in estimates from external file
M2 <- read.csv(file.path(data.source,
paste0("M_NSAS_smoothedSpan50_notExtrapolated_sms",ac(SMSkeyRuns),".csv")))
colnames(M2) <- sub("X","",colnames(M2))
rownames(M2) <- M2[,1]
M2 <- M2[,-1]# Trim off first column as it contains 'ages'
M2 <- M2[,apply(M2,2,function(x){all(is.na(x))==F})] # keep only years with data
#Extract key data from default assessment
NSHM2 <- NSH
NSHM2@m[] <- NA
yrs <- dimnames(NSHM2@m)$year
yrs <- yrs[which(yrs %in% colnames(M2))]
NSHM2@m[rownames(M2),yrs][] <- as.matrix(M2)
#- M extension to previous and future years (relative to vector given by WGSAM).
# One uses a 5 year running average for that
# year discrepency between assessment and M from WGSAM
extryrs <- dimnames(NSHM2@m)$year[which(!dimnames(NSHM2@m)$year %in% yrs)]
# year after those avaible for M from WGSAM
extryrsfw <- extryrs[which(extryrs > max(an(yrs)))]
# years prior to those available for M from WGSAM
extryrsbw <- extryrs[which(extryrs <= max(an(yrs)))]
ages <- dimnames(NSHM2@m)$age
extrags <- names(which(apply(M2,1,function(x){all(is.na(x))==T})==T))
yrAver <- 5
for(iYr in as.numeric(rev(extryrs))){
for(iAge in ages[!ages%in%extrags]){
if(iYr %in% extryrsbw){
NSHM2@m[ac(iAge),ac(iYr)] <- yearMeans(NSHM2@m[ac(iAge),ac((iYr+1):(iYr+yrAver)),],na.rm=T)
}
if(iYr %in% extryrsfw){
NSHM2@m[ac(iAge),ac(iYr)] <- yearMeans(NSHM2@m[ac(iAge),ac((iYr-1):(iYr-yrAver)),],na.rm=T)
}
}
}
if(length(extrags)>0){
for(iAge in extrags)
NSHM2@m[ac(iAge),] <- NSHM2@m[ac(as.numeric(min(sort(extrags)))-1),]
}
#Write new M values into the original stock object
#addM <- 0.11 #M profiling based on 2018 benchmark meeting
#addM <- 0
NSH@m <- NSHM2@m
### ============================================================================
### Prepare index object for assessment
### ============================================================================
## Load and modify all numbers at age data
NSH.tun <- readFLIndices(file.path(data.source,"fleet.txt"))
NSH.tun <- lapply(NSH.tun,function(x) {x@type <- "number"; return(x)})
NSH.tun[["IBTS0"]]@range["plusgroup"] <- NA
### ============================================================================
### Apply plus group to all data sets
### ============================================================================
pg <- max(an(rownames(M2)))
if(pg > 8) pg <- 8
#pg <- 8
#- This function already changes the stock and landings.wts correctly
NSH <- setPlusGroup(NSH,pg)
NSH.tun[["HERAS"]]@index[ac(pg),] <- quantSums(NSH.tun[["HERAS"]]@index[ac(pg:dims(NSH.tun[["HERAS"]]@index)$max),])
NSH.tun[["HERAS"]] <- trim(NSH.tun[["HERAS"]],age=dims(NSH.tun[["HERAS"]]@index)$min:pg)
NSH.tun[["HERAS"]]@range["plusgroup"] <- pg
NSH.tun[["IBTS-Q3"]] <- trim(NSH.tun[["IBTS-Q3"]],age=0:5)
NSH.tun[["IBTS-Q1"]] <- trim(NSH.tun[["IBTS-Q1"]],age=1)
NSH.tun[["HERAS"]] <- trim(NSH.tun[["HERAS"]],age=1:pg)
### ============================================================================
### Closure data deletion
### ============================================================================
## We don't believe the closure catch data, so put it to NA
NSH@catch.n[,ac(1978:1979)] <- NA
### ============================================================================
### Now the multifleet data
### ============================================================================
assessment_name <- paste0(assessment_name,'_mf')
# read the catch data
caaF <- read.table(file.path(data.source, "canum_mf.txt"),stringsAsFactors=F,header=T)
caA <- matrix(subset(caaF,Fleet=="A")[,"numbers"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))*1000
caB <- matrix(subset(caaF,Fleet=="B")[,"numbers"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))*1000
caC <- matrix(subset(caaF,Fleet=="C")[,"numbers"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))*1000
caD <- matrix(subset(caaF,Fleet=="D")[,"numbers"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))*1000
cwA <- matrix(subset(caaF,Fleet=="A")[,"weight"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))
cwB <- matrix(subset(caaF,Fleet=="B")[,"weight"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))
cwC <- matrix(subset(caaF,Fleet=="C")[,"weight"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))
cwD <- matrix(subset(caaF,Fleet=="D")[,"weight"],nrow=length(0:9),ncol=length(1997:range(NSH)["maxyear"]),dimnames=list(age=0:9,year=1997:range(NSH)["maxyear"]))
#- Functions to set the plusgroups
setMatrixPlusGroupN <- function(x,pg){
idxpg <- which(an(rownames(x)) == pg)
x[idxpg,] <- colSums(x[idxpg:nrow(x),],na.rm=T)
x <- x[1:idxpg,]
return(x)}
setMatrixPlusGroupWt <- function(x,y,pg){
idxpg <- which(an(rownames(x)) == pg)
y[idxpg,] <- colSums(x[idxpg:(nrow(x)),]*y[idxpg:nrow(y),],na.rm=T) / colSums(x[idxpg:nrow(x),],na.rm=T)
y[is.nan(y)] <- 0
y <- y[1:idxpg,]
return(y)}
NSH3f <- FLCore::expand(NSH,area=c("A","BD","C"))
NSH3f@catch.n[,ac(1997:range(NSH)["maxyear"]),,,"A"] <- setMatrixPlusGroupN(caA,pg)
NSH3f@catch.n[,ac(1997:range(NSH)["maxyear"]),,,"BD"] <- setMatrixPlusGroupN(caB+caD,pg)
NSH3f@catch.n[,ac(1997:range(NSH)["maxyear"]),,,"C"] <- setMatrixPlusGroupN(caC,pg)
NSH3f@catch.wt[,ac(1997:range(NSH)["maxyear"]),,,"A"] <- setMatrixPlusGroupWt(caA,cwA,pg)
NSH3f@catch.wt[,ac(1997:range(NSH)["maxyear"]),,,"BD"] <- setMatrixPlusGroupWt(caB+caD,(cwB*caB + cwD*caD)/(caB+caD),pg)
NSH3f@catch.wt[,ac(1997:range(NSH)["maxyear"]),,,"C"] <- setMatrixPlusGroupWt(caC,cwC,pg)
NSH3f@landings.wt <- NSH3f@catch.wt
NSH3f@catch.n[,ac(1947:1996),,,-1] <- NA
NSH3f@catch.wt[,ac(1947:1996),,,-1] <- NA
NSH3f@name <- assessment_name
NSHsum <- NSH[,ac(1947:1996)]
NSHres3 <- NSH3f[,ac(1997:range(NSH)["maxyear"])];
NSHres3@landings.n[] <- NSHres3@catch.n
NSHres3@catch <- computeCatch(NSHres3)
NSHs3 <- FLStocks(residual=NSHres3,sum=NSHsum)
### ============================================================================
### Save objects
### ============================================================================
save(NSH, NSH.tun,
file=file.path(data.source,paste0("data_sms",SMSkeyRuns,"_sf.RData")))
save(NSHs3, NSH3f, NSH.tun,
file=file.path(data.source,paste0("data_sms",SMSkeyRuns,"_mf.RData")))
}
## Run analysis, write model results
## Before: data_sf.RData and data_mf.RData (data)
## After: config_sf.RData, results_sf.RData,
## config_mf.RData and results_mf.RData (model)
rm(list=ls())
library(FLCore)
library(stockassessment)
library(FLSAM)
library(methods)
mkdir("model")
mkdir(file.path("model","config"))
mkdir(file.path("model","assessment"))
mkdir(file.path("model","packages"))
data.source <- file.path(".", "data")
model.save <- file.path(".",'model')
runName <- '01_baseline'
SMSkeyRuns <- 2019
source(file.path('functions','config_baseline_sf.R'))
source(file.path('functions','config_baseline_mf.R'))
### ============================================================================
### single fleet
### ============================================================================
load(file.path(data.source,paste0("data_sms",SMSkeyRuns,"_sf.RData"))) # NSH, NSH.tun
# create control
pg <- range(NSH)[2]
NSH.ctrl <- config_baseline_sf(NSH,NSH.tun,pg)
# run model
NSH.sam <- FLSAM(NSH, NSH.tun, NSH.ctrl)
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(NSH, NSH.tun, NSH.ctrl, NSH.sam,
file=file.path(model.save,
paste0("results_",runName,"_sf.RData")))
### ============================================================================
### multi fleet
### ============================================================================
load(file.path(data.source,paste0("data_sms",SMSkeyRuns,"_mf.RData")))
NSH3.ctrl <- config_baseline_mf(NSHs3,NSH.tun)
load(file.path(data.source,'results_mf_2020.Rdata'))
NSH3f.sam.stk0 <- NSH3f.sam
# run model
NSH3f.sam <- FLSAM(NSHs3,
NSH.tun,
NSH3.ctrl,
starting.values = NSH3f.sam.stk0)
save(NSHs3,
NSH.tun,
NSH3.ctrl,
NSH3f.sam,
file=file.path(model.save,
paste0("results_",runName,"_mf.RData")))
......@@ -3,7 +3,7 @@ library(TMB)
rm(list=ls())
path <- "J:/git/sanoba_nsas"
path <- "C:/git/sanoba_nsas"
try(setwd(path),silent=TRUE)
......@@ -32,7 +32,7 @@ figures.dir <- file.path(".","figures",runName)
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)
SST_pred <- seq(from=min(dat$sst),to=max(dat$sst),by=0.5)#min(dat$sst)#
################################################
# compile TMB function
......@@ -60,4 +60,4 @@ 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
lines(SSB_pred, logR-2*rep$sd, lty="dotted")
......@@ -2,7 +2,7 @@ library(FLCore)
rm(list=ls())
path <- "J:/git/sanoba_nsas"
path <- "C:/git/sanoba_nsas"
try(setwd(path),silent=TRUE)
......
#-------------------------------------------------------------------------------
# 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')))
age,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016
0,0.708213797,0.700420308,0.695702476,0.694734367,0.69649246,0.700220877,0.70721726,0.717491371,0.728355514,0.749058092,0.770926698,0.78341299,0.793727627,0.798940199,0.79555435,0.786964879,0.780578111,0.774778919,0.766525787,0.76108545,0.75534646,0.751198831,0.754413057,0.761415324,0.770958101,0.785369329,0.80382246,0.820338224,0.836227082,0.853640588,0.867861105,0.882440395,0.895153744,0.902616891,0.907615667,0.906022867,0.897765101,0.884524806,0.865860971,0.842214289,0.813413851,0.779131145,0.739540933
1,0.550199821,0.5733306,0.592828615,0.608218668,0.620141615,0.629109162,0.634325553,0.635818975,0.635337052,0.629126681,0.620425865,0.613087129,0.603888468,0.594452378,0.585667024,0.577216982,0.568919372,0.557081878,0.542547212,0.53193502,0.522277699,0.515338091,0.51669461,0.522309448,0.531473583,0.551205581,0.578789176,0.598913092,0.613712417,0.629189059,0.635831967,0.625124986,0.611273103,0.598987052,0.582702835,0.571579019,0.564276237,0.558272358,0.555707351,0.554592462,0.556119022,0.561635026,0.570009206
2,0.32851696,0.334817623,0.339680285,0.343027053,0.344954203,0.345531942,0.344443746,0.341743304,0.338056128,0.332887823,0.326609971,0.317536405,0.305990404,0.296910608,0.289645711,0.282487026,0.278787135,0.281021444,0.286585729,0.291213116,0.296292759,0.301901874,0.304883562,0.307545114,0.311848035,0.319140474,0.327679672,0.334790765,0.342878275,0.352360335,0.357582503,0.357093986,0.355199279,0.350617709,0.343419373,0.33860892,0.335584319,0.332726985,0.331007152,0.329510786,0.328891292,0.329889226,0.331901228
3,0.276135716,0.281502707,0.285333605,0.287428971,0.288000101,0.287260648,0.284841137,0.280761155,0.275817947,0.268207247,0.259420833,0.249299341,0.237202795,0.228151589,0.222774789,0.218739111,0.216621421,0.218960719,0.224640262,0.228753302,0.231205071,0.233476387,0.232821012,0.230763311,0.231165841,0.234720802,0.23941668,0.245060456,0.254144581,0.265680021,0.274251523,0.281723484,0.28820693,0.289276128,0.287923452,0.287230747,0.286962511,0.286031631,0.284932414,0.283192818,0.281167762,0.279318838,0.277321134
4,0.250468966,0.252641953,0.253829565,0.253913336,0.252969707,0.251132181,0.248328085,0.244541692,0.239991928,0.233817372,0.226628482,0.216615482,0.204199835,0.195155133,0.190110184,0.186426425,0.184976389,0.189017384,0.197034674,0.202742727,0.206967509,0.210561529,0.209307685,0.205897445,0.204447948,0.204006086,0.202951653,0.204401146,0.21062603,0.21941294,0.227151127,0.236981609,0.246731097,0.251487051,0.254867387,0.258085581,0.261064038,0.263416982,0.265354474,0.266678961,0.267530215,0.268123229,0.268327195
5,0.23447327,0.234464139,0.233861471,0.232578976,0.230672586,0.228241153,0.225275081,0.221754538,0.217763017,0.212802755,0.207195173,0.199647449,0.190494457,0.183490008,0.178962645,0.175292482,0.173310678,0.174870175,0.178922824,0.181948145,0.185006551,0.187800363,0.187011801,0.18494938,0.184494938,0.184968889,0.18521916,0.187418857,0.193271614,0.201153948,0.208248293,0.217045732,0.225845828,0.23089188,0.23497111,0.238794552,0.242298961,0.245269791,0.247835695,0.249876104,0.251501253,0.252876465,0.253900936
6,0.226462538,0.224001421,0.221304042,0.218318647,0.215064072,0.211604753,0.208007389,0.204243177,0.200233304,0.195739178,0.190870709,0.184522955,0.177070479,0.171176467,0.167144044,0.163817573,0.161672096,0.161910384,0.163852952,0.165151454,0.165969652,0.166805647,0.165960358,0.164463604,0.164561907,0.166078056,0.167968295,0.17117933,0.177158995,0.184987426,0.192059792,0.199895362,0.207695325,0.212850019,0.217261287,0.221614441,0.225830096,0.229732183,0.233475562,0.236915667,0.240143655,0.243292987,0.24627746
7,0.221422134,0.21794858,0.214359479,0.210619119,0.206735653,0.202754877,0.198736314,0.194656255,0.19043565,0.185868624,0.181063495,0.175046746,0.168154912,0.162745298,0.159091331,0.156145955,0.154285756,0.154660709,0.156657102,0.15800598,0.158745221,0.159565258,0.159242558,0.158461334,0.158793092,0.1599709,0.161311745,0.163763015,0.16833703,0.174321429,0.180076572,0.186471046,0.193014452,0.198221219,0.203211173,0.208214298,0.213195103,0.21810863,0.223017538,0.227863339,0.232680097,0.237519843,0.242350175
8,0.220713822,0.217241607,0.213626478,0.209812115,0.205831485,0.201753089,0.197606373,0.193367361,0.189024418,0.184166361,0.179081529,0.173343592,0.16701113,0.161849156,0.158164357,0.155137265,0.152843668,0.151961256,0.152238204,0.152356031,0.151754294,0.151438919,0.151223186,0.151034866,0.152076771,0.154901376,0.158718498,0.162716969,0.167417467,0.173021874,0.178240828,0.183792981,0.189268346,0.193115466,0.196438289,0.200237595,0.204448886,0.20873563,0.213270568,0.21788753,0.222670873,0.227739228,0.233012502
9,0.219419596,0.216781853,0.213809052,0.210470153,0.206770713,0.20275025,0.198442117,0.1938284,0.188871723,0.183251092,0.177197296,0.169730186,0.161164573,0.154233186,0.148995891,0.144262006,0.140888453,0.14021771,0.141407798,0.142038943,0.142086603,0.142379894,0.141312901,0.139652715,0.139819575,0.142165872,0.14540599,0.149503149,0.155907903,0.164184008,0.17146449,0.178806246,0.185929544,0.190503207,0.194201106,0.198158887,0.202280415,0.206195919,0.210111591,0.213830154,0.217471608,0.221197249,0.22489531
"","1960","1961","1962","1963","1964","1965","1966","1967","1968","1969","1970","1971","1972","1973","1974","1975","1976","1977","1978","1979","1980","1981","1982","1983","1984","1985","1986","1987","1988","1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018","2019"
"0",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.709921413147949,0.709826388591418,0.71208933596006,0.71756483140207,0.724962006684904,0.733556767852147,0.744629141708122,0.758121532017125,0.771303410262707,0.791362174597277,0.818320925889299,0.838702689696411,0.849313372339979,0.855909482295313,0.85843994835096,0.853127027940179,0.841641790874108,0.832106384508114,0.820299652113761,0.80326997114353,0.791043790224165,0.780340266694331,0.771977070835683,0.773394893887273,0.779433575178022,0.787366054883963,0.800281299543179,0.817951682616484,0.832675254664011,0.846038791703831,0.8615750626693,0.874501107493473,0.887030033664734,0.900364302400133,0.908204984225724,0.910436022236769,0.909866322615367,0.904584400517835,0.89474699405298,0.881178548774983,0.862969141008219,0.840023155844624,0.812798261794192,0.781214460308385,0.744980327441459,0.704250035055698
"1",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.471728637881525,0.492967856343673,0.511587999453601,0.527399180374252,0.540619947277215,0.551416559766628,0.559561676074361,0.565070731338135,0.568474150520857,0.568928434518853,0.566248145579157,0.561950800690747,0.553303541297621,0.5406285209331,0.530025774720142,0.521731455424586,0.513138412384935,0.506081725583136,0.49940196700942,0.492565658772043,0.488346806029996,0.48261449928481,0.479463244538165,0.485253653743783,0.494829729079409,0.506028684290285,0.526921974472893,0.555575038444979,0.574816065533686,0.584753538921638,0.593959079164122,0.597955581647271,0.591374845934748,0.577676627790935,0.565625756903042,0.554872410490515,0.54198534626784,0.531107264664795,0.521846662299363,0.512029784277938,0.503125953762908,0.495223016220327,0.487619734878851,0.480601809544772,0.474555806704829,0.469124361303812
"2",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.296392973965184,0.301815504094776,0.306256219497473,0.309649002315672,0.312056992998491,0.313533034423706,0.313942540923285,0.313305490989416,0.3119065737402,0.309423971072372,0.305750168613817,0.301491311079362,0.294036565513849,0.284073440123834,0.277155591366714,0.273986986285616,0.27180726018989,0.270973571936269,0.272841221012527,0.276749691086977,0.28003241728716,0.281975160393621,0.284812742322829,0.288838116459624,0.293387883156662,0.298784539697097,0.307454497272272,0.318197479899312,0.325895954100033,0.33182525155671,0.338328398344993,0.341889402576769,0.340661943093467,0.336771900712707,0.332661586510658,0.327327024781468,0.320277064447261,0.314723253526169,0.310483474520752,0.305831163344716,0.301661198132917,0.29800668657175,0.294452871780968,0.291222835045461,0.28863670894046,0.286431214983172
"3",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.269416624375446,0.272663475701948,0.274941357228651,0.276107228710793,0.276324319441853,0.275724096144664,0.274166333452517,0.271650709762875,0.268518547444832,0.264164026290329,0.258459511458872,0.2524636033709,0.243710356632656,0.232694418972307,0.224915280493173,0.221608253695254,0.219929319505116,0.219258970665935,0.221059485467143,0.225050757501953,0.227999702469222,0.228441672311039,0.229494107481884,0.231965396296555,0.234847938457495,0.239057799513604,0.246397364609169,0.255519123046274,0.262614296610513,0.269859926874748,0.27862138452767,0.283872385028623,0.283825847851519,0.281448661926847,0.278845496898809,0.274677363329874,0.268700052753795,0.264660406271702,0.262329560511742,0.259682531112042,0.257771977688701,0.256626164885734,0.255774927061536,0.255523278403773,0.256320252948006,0.25781184285594
"4",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.2548493763822,0.253006003926819,0.250780155034322,0.2480122953896,0.244883956728261,0.241546878579024,0.237894442917284,0.233914600463324,0.229900943863203,0.225238660958449,0.219804254619566,0.214628290036073,0.208453994210787,0.201345990803881,0.196316939990913,0.195178625125575,0.19610530768606,0.196653365034677,0.197008278700776,0.198195509815172,0.199031405142373,0.197278756131041,0.196013337233711,0.196601959358774,0.197239188196585,0.199952477190409,0.206937377751956,0.216413275910719,0.224443377036632,0.233779455520721,0.245472999555492,0.253036785213104,0.254731427142719,0.254218001944305,0.253109978616231,0.249977150321792,0.244841021801767,0.241492629420751,0.239730498028851,0.237503333964369,0.235827445824194,0.23472236749636,0.233743435498734,0.23321066237295,0.233614832583837,0.234586595626584
"5",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.235236570783192,0.233242646083143,0.230956937476693,0.228277529366476,0.225308749969553,0.222148944907754,0.218735405574499,0.215059426124899,0.211300356323323,0.207063225013875,0.202260858604742,0.197534303206416,0.191475214001953,0.184360267414968,0.179358041416786,0.177983906678266,0.178333673726225,0.178446580315929,0.178883438110543,0.180385337049957,0.181261932304527,0.179946366087793,0.17898071711658,0.178533556900502,0.177573298821532,0.17877276561348,0.183510432856329,0.189998301961351,0.196159484409293,0.204785186338715,0.216045140181864,0.223930246860015,0.227474414546427,0.229852641211785,0.231267462473509,0.230455554360649,0.227929400023327,0.226555245760563,0.226176759414282,0.225254014634103,0.224577866430582,0.224166394461238,0.22367154444172,0.22333497046978,0.223537069748537,0.223996512772068
"6",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.226790676232422,0.223092332030349,0.219379736862161,0.21562173514262,0.211828976988414,0.208042402789847,0.20430907393585,0.20060917332407,0.196883323607885,0.193150793087087,0.18937202357216,0.185399397982519,0.180120308358457,0.174004434054065,0.169277588959446,0.166550650689558,0.164647369906335,0.163112870070341,0.162170641855727,0.161914934116723,0.161697007584014,0.160518301020245,0.159916986477407,0.160254594163257,0.160757418483391,0.162878966399282,0.167597114141485,0.173832335072889,0.179952031853466,0.188448577742565,0.199307954185902,0.207137976122993,0.211336184591924,0.214710160581286,0.216982778956353,0.216988575748974,0.215393393788248,0.21466706813533,0.214684729955431,0.214101955122025,0.213559245758666,0.213063186573613,0.21234187094656,0.211618064008208,0.211236142842971,0.210941388097937
"7",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.222042827715914,0.217218480455241,0.212526699527718,0.207939496750029,0.203470906915947,0.199158004026946,0.195046452985936,0.19111704055914,0.187309523543002,0.18363369560465,0.180062156874596,0.17647887384436,0.172337692786757,0.167882210286677,0.164181018303751,0.161505212495675,0.159355632850881,0.157645487825109,0.156484076035187,0.155812093738129,0.15528984584995,0.154149846447417,0.153456139571118,0.153418662674261,0.153528512880303,0.155115029556374,0.158772804963683,0.16364334842475,0.168876374898667,0.17653795324499,0.186291493822537,0.193682723839292,0.198726182220381,0.203647772045644,0.207310266755129,0.208677328564758,0.208706587175766,0.209251767256171,0.210230432705222,0.210611765185002,0.210826916007775,0.210868016106532,0.21056434665471,0.210102784325237,0.209764822302907,0.209343893466986
"8",NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,0.222878498269777,0.217608170016783,0.212393632108218,0.207159875795553,0.201989065928003,0.196959189235232,0.192065590389752,0.187289867959511,0.182689389517322,0.178007580783677,0.173201469868345,0.168648827928049,0.163808701984887,0.158743935386049,0.154737993775898,0.152402306512936,0.151054391416269,0.15003684682117,0.149499095695815,0.149608142011315,0.149732955607346,0.149338071132228,0.149298705685537,0.149677580209853,0.150219000483615,0.151930337277169,0.155251814288771,0.159504883971799,0.164018760229285,0.17035221859375,0.178270869481044,0.184414179119196,0.188805835969803,0.193139363125494,0.196558978760842,0.198326081503052,0.199133907165352,0.200274429701353,0.201675670683879,0.202608970975081,0.20338128346257,0.203983289379358,0.204298276573988,0.204472107857364,0.20471572244753,0.204870013779705
This diff is collapsed.
......@@ -15,7 +15,14 @@ Type objective_function<Type>::operator() ()
PARAMETER(logSigma);
vector<Type> pred_model=logA+log(SSB)-exp(logB)*SSB+exp(logC)*SST;
Type nl=-sum(dnorm(logR,pred_model,exp(logSigma),true));
vector<Type> pred=logA+log(SSB_pred)-exp(logB)*SSB_pred+exp(logC)*SST_pred;
matrix<Type> pred(SSB_pred.size(),SST_pred.size());
/*pred=logA+log(SSB_pred)-exp(logB)*SSB_pred+exp(logC)*SST_pred;*/
for(int iSST=0; iSST<SST_pred.size(); ++iSST){
for(int iSSB=0; iSSB<SSB_pred.size(); ++iSSB){
pred(iSSB,iSST)=logA+log(SSB_pred(iSSB))-exp(logB)*SSB_pred(iSSB)+exp(logC)*SST_pred(iSST);
}
}
ADREPORT(pred)
/*ADREPORT(pred_model)*/
return nl;
......
No preview for this file type
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_VECTOR(SSB_pred)
DATA_VECTOR(SST_pred)
DATA_VECTOR(SSB)
DATA_VECTOR(SST)
DATA_VECTOR(logR);
PARAMETER(logA);
PARAMETER(logB);
PARAMETER(logC);
PARAMETER(logSigma);
vector<Type> pred_model=logA+log(SSB)-exp(logB)*SSB+exp(logC)*SST;
Type nl=-sum(dnorm(logR,pred_model,exp(logSigma),true));
vector<Type> pred=logA+log(SSB_pred)-exp(logB)*SSB_pred+exp(logC)*SST_pred;
ADREPORT(pred)
/*ADREPORT(pred_model)*/
return nl;
}
config_baseline_mf <- function( NSHs3,
NSH.tun){
NSH3.ctrl <- FLSAM.control(NSHs3,NSH.tun,sumFleets=dimnames(NSHs3[["residual"]]@catch)$area)
catchRow <- grep("catch",rownames(NSH3.ctrl@f.vars))
NSH3.ctrl@states["catch A",] <- c(-1,0:6,6)
NSH3.ctrl@states["catch BD",] <- c(7:9,10,10,10,rep(-1,3))
NSH3.ctrl@states["catch C",] <- c(-1,11:13,14,14,14,rep(-1,2))
NSH3.ctrl@catchabilities["HERAS",ac(1:8)] <- c(101,