Commit d6f616c3 authored by Iago Mosqueira's avatar Iago Mosqueira
Browse files

Added eqsim refpts calculation

parent a555af2a
No preview for this file type
# model.R - DESC
# /model.R
# Copyright Iago MOSQUEIRA (WMR), 2019
......@@ -57,46 +58,114 @@ years <- setNames(as.list(seq(2018, by=-1, length=6)),
retro <- FLStocks(parallel::mclapply(years, function(x) {
aap(window(stock, end=x), lapply(indices[c("GAM", "SNS")], window, end=x),
control=control, stdfile=btsgam@stdfile, wkdir=tempfile()) + window(stock, end=x)
control=control, stdfile=btsgam@stdfile, wkdir=tempfile()) +
window(stock, end=x)
}, mc.cores=length(years)))
# --- STOCKS & RESIDUALS
stocks <- FLStocks(lapply(runs, "+", stock))
residuals <- lapply(runs[!names(runs) %in% "McMC"], residuals, stock)
# --- SRR
srrs <- FLSRs(setNames(lapply(c("bevholt", "ricker", "segreg"), function(x)
fmle(as.FLSR(stocks[["BTS-GAM"]], model=x, name=x))), c("bevholt", "ricker", "segreg")))
srrs <- setNames(lapply(c("bevholt", "ricker", "segreg"), function(x)
as.FLSR(stocks[["BTS-GAM"]], model=x, name=x)),
c("bevholt", "ricker", "segreg"))
# HACK bevholt
bhlogl <- function(a, b, rec, ssb) {
if (b<=0)
return(-1e100)
else
return(loglAR1(log(rec), log(a*ssb/(b+ssb))))
}
logl(srrs[[1]]) <- bhlogl
srrs <- FLSRs(lapply(srrs, fmle))
# --- REFERENCE points
srfit <- eqsr_fit(stocks[["BTS-GAM"]],
nsamp = 1000,
models = c("Ricker", "Bevholt"))
eqsr_plot(srfit, n=2e4, Scale=1e3)
eqsr_plot(srfit, n=2e4, ggPlot=TRUE, Scale=1e3)
srsim <- eqsim_run(srfit,
bio.years = c(2004, 2018),
sel.years = c(2010, 2018),
Fcv = 0.24,
Fphi = 0.42,
# Blim, Bpa
Blim = 20000,
Bpa = 22000,
library(msy)
Fs <- seq(0, 1.5, length=51)
nsamp <- 1000
# FIT segreg to obtain BLIM & BPA
srfit0 <- eqsr_fit(stocks[["BTS-GAM"]],
nsamp = nsamp, models = c("Segreg"))
Blim <- srfit0[["sr.det"]][,"b"]
Bpa <- 1.4 * Blim
# FIT all models
srfit1 <- eqsr_fit(stocks[["BTS-GAM"]],
nsamp = nsamp,
models = c("Ricker", "Bevholt", "Segreg"))
# SIMULATE all models w/10 y, Fcv=Fphi=0, Btrigger=0
srsim1 <- eqsim_run(srfit1,
bio.years = c(2009, 2018), sel.years = c(2009, 2018),
Fcv = 0, Fphi = 0,
Btrigger=0, Blim = Blim, Bpa = Bpa,
Fscan = Fs,
verbose = FALSE)
# EXTRACT Flim and Fpa
Flim <- srsim1$Refs2["catF", "F50"]
Fpa <- Flim / 1.4
# FIT all models, REMOVE last 3 years
srfit2 <- eqsr_fit(stocks[["BTS-GAM"]],
nsamp = nsamp,
models = c("Ricker", "Bevholt", "Segreg"),
remove.years=ac(2016:2018))
# SIMULATE w/ 5 yrs, Fcv=0.212, Fphi=0.423 (WKMSYREF4)
srsim2 <- eqsim_run(srfit2,
bio.years = c(2013, 2018), sel.years = c(2013, 2018),
bio.const = FALSE, sel.const = FALSE,
Fcv=0.212, Fphi=0.423,
Btrigger=0, Blim = Blim, Bpa = Bpa,
Fscan = Fs,
verbose = FALSE)
Fmsy <- srsim$Refs2["lanF", "medianMSY"]
F05 <- srsim$Refs2["catF", "F05"]
# SIMULATE for Btrigger
srsim3 <- eqsim_run(srfit2,
bio.years = c(2009, 2018), sel.years = c(2009, 2018),
bio.const = FALSE, sel.const = FALSE,
Fcv = 0, Fphi = 0,
Btrigger=0, Blim = Blim, Bpa = Bpa,
Fscan = Fs,
verbose = FALSE)
x <- srsim3$rbp[srsim3$rbp$variable=="Spawning stock biomass", ]
Btrigger <- x[which(abs(x$Ftarget - Fmsy) == min(abs(x$Ftarget - Fmsy))), "p05"]
# SIMULATE
srsim4 <- eqsim_run(srfit,
bio.years = c(2009, 2018), sel.years = c(2009, 2018),
bio.const = FALSE, sel.const = FALSE,
Fcv=0.212, Fphi=0.423,
Btrigger=Btrigger, Blim = Blim, Bpa = Bpa,
Fscan = seq(0, 1.2, len = 40),
verbose = FALSE)
eqsim_plot(srsim, catch=TRUE)
Fmsy <- srsim4$Refs2["lanF", "medianMSY"]
F05 <- srsim4$Refs2["catF", "F05"]
eqsim_plot_range(srsim, type="median")
eqsim_plot_range(srsim, type="ssb")
rps <- FLPar(Fmsy=Fmsy, F05=F05, Flim=Flim, Fpa=Fpa,
Btrigger=Btrigger, Blim=Blim, Bpa=Bpa, units=c(rep("", 3), rep("tonnes", 3)))
# --- RUN short term forecast
# --- SAVE
save(runs, stocks, residuals, retro, runmcmc, srrs, srfit, srsim,
save(runs, stocks, residuals, retro, runmcmc, srrs, rps, srfit2, srsim4,
file="model/aap.RData", compress="xz")
No preview for this file type
Supports Markdown
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