gadget_projections.RdSetup model and parameter files for forward simulations and deterministic projections.
gadget_project_time(
path = ".",
num_years = 100,
variant_dir = getwd() %>% stringr::str_count("/") %>% rep("../", .) %>%
paste(collapse = "") %>% paste(tempdir(), sep = "")
)
gadget_project_stocks(
path,
imm.file,
mat.file,
spawnfunction = "hockeystick",
proportionfunction = "constant 1",
mortalityfunction = "constant 0",
weightlossfunction = "constant 0"
)
gadget_project_fleet(
path,
pre_fleet = "comm",
post_fix = "pre",
fleet_type = "linearfleet",
common_mult = NULL,
pre_proportion = NULL,
type = "standard",
...
)
gadget_project_prognosis_likelihood(
path,
stocks,
pre_fleets = "comm",
post_fix = "pre",
mainfile_overwrite = TRUE,
firsttacyear = NULL,
assessmentstep = 2,
weightoflastyearstac = 0,
...
)
gadget_project_recruitment(
path,
stock,
recruitment = NULL,
n_replicates = 100,
params.file = "PRE/params.pre",
method = "AR",
...
)
gadget_project_advice(
path,
params.file = "PRE/params.pre",
harvest_rate = 0.2,
advice_cv = 0.2,
advice_rho = 0.6,
pre_fleet = "comm",
post_fix = "pre",
n_replicates = 100
)
gadget_project_ref_points(path, ref_points, params.file = "PRE/params.pre")
gadget_project_output(
path,
imm.file,
mat.file,
pre_fleets = "comm",
post_fix = "pre",
f_age_range = NULL,
rec_age = NULL,
print_block = "1"
)gadget variant director
number of years
location of the
location of the immature stock file
location of the mature stock file
what spawn function to use (only hockeystick atm)
proportion suitability
mortataliy suitability
weightloss suitability
name of the fleet projections are based on (original fleet). Only a single fleet should be implemented here at a time. Multiple fleets require multiple gadget_project_advice calls.
A single label (vector of length 1) used to distinguish the original fleet and the projection fleet (with same parameterisations). Should be the same for all projection fleets.
type of gadget fleet for the projections
a string with a base name for a vector of harvest/catch rate multipliers (shared between fleets)
proportion of the effort taken
type of projection, either "standard" or "prognosis"
additional arguments for the
names of of the stocks
vector of fleets on which the projections are based
should the likelihood be overwritten, defaults to TRUE, but should be set as FALSE for multiple prognosis components (different groups of fleets)
when does the TAC scheme start
when is the assessment
running average TAC..
name of immature stock
recruitment time series
number of simulations
name of the parameter files
prediction method, built in options are 'AR', 'bootstrap' or constant
median harvest rate
assessment error cv
assessment error correlation
tibble with reference points
F age range, specified in the format a1:a2
age of recruitment (as reported)
file suffix for the printfile
gadget variant directory
gadget project recruitment: Sets up the process error for the stock recruitment relationship. The user has a choice of three types, constant, AR and block bootstrap based on a time-series of recruitment
gadget project stocks: Here the user can define the SSB - Rec relationship used in the projections. It defines link between two stocks, immature and mature stock and in the case of minage (immstock) > 0 an interim bookeeping stock from age 0 to minage.
Limitiations: At present theres is no way to define recruitment into multiple stocks (ie. via multiple interim stocks). This function also assumes that there is no SSB - Rec relationship, since if it is already present the simulation should use that. Multi-annual recruitment is currently not defined.
if (FALSE) {
fit <- gadget.fit()
res <-
gadget_project_time() %>%
gadget_project_stocks(imm.file = 'Modelfiles/cod.imm',mat.file = 'Modelfiles/cod.mat') %>%
gadget_project_fleet(pre_fleet = 'comm') %>%
gadget_evaluate(params.out = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
params.in = 'WGTS/params.final') %>%
gadget_project_recruitment(stock = 'codimm',
recruitment = fit$stock.recruitment %>%
filter(stock == 'codimm',
year > 1980),
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
gadget_project_ref_points(ref_points = tibble(codmat.blim = 207727665),
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
gadget_project_advice(pre_fleet = 'comm',
harvest_rate = 1:100/100,
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
gadget_project_output(imm.file = 'Modelfiles/cod.imm',mat.file = 'Modelfiles/cod.mat',
pre_fleet = 'comm') %>%
gadget_evaluate(params.in = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
{read.printfiles(paste(attr(.,'variant_dir'),'out',sep='/'))} %>%
map(mutate, trial=cut(1:length(year),c(0,which(diff(year)<0),1e9),labels = FALSE)) %>%
set_names(c("catch.F","catch.lw",'codimm.rec','codmat.ssb')) %>%
map(left_join,tibble(trial=1:10000,harvest_rate = rep(1:100/100,100)))
yield_curve <-
res$catch.lw %>%
filter(year>2050) %>%
group_by(trial,harvest_rate,year) %>%
summarise(c=sum(biomass_consumed)/1e6) %>%
group_by(harvest_rate) %>%
summarise(m=median(c),u=quantile(c,0.95),l=quantile(c,0.05))
ssb_curve <-
res$codmat.ssb %>%
filter(year>2050) %>%
group_by(trial,harvest_rate,year) %>%
summarise(c=sum(number*mean_weight)/1e6) %>%
group_by(harvest_rate) %>%
summarise(m=median(c),u=quantile(c,0.95),l=quantile(c,0.05))
f.curve <-
res$catch.F %>%
filter(year>2050) %>%
group_by(trial,harvest_rate,year) %>%
summarise(c=median(mortality)) %>%
group_by(harvest_rate) %>%
summarise(m=median(c),u=quantile(c,0.95),l=quantile(c,0.05))
blim <-
fit$res.by.year %>%
filter(grepl('mat',stock)) %>%
summarise(b=min(total.biomass)/1e6) %>%
.$b
bpa <- 1.4*blim
hr_msy <-
yield_curve %>% filter(m==max(m)) %>% .$harvest_rate
hr_lim <-
ssb_curve %>%
filter(m>blim) %>%
filter(harvest_rate == max(harvest_rate)) %>%
.$harvest_rate
f.msy <-
f.curve %>%
filter(harvest_rate == hr_msy) %>%
.$m
f.lim <-
f.curve %>%
filter(harvest_rate == hr_lim) %>%
.$m
f.pa <-
f.lim/1.4
hr_pa <-
f.curve %>%
filter(m < f.pa) %>%
summarise(hr = max(harvest_rate)) %>%
.$hr
library(patchwork)
yield_curve %>%
left_join(f.curve %>%
select(harvest_rate,F=m)) %>%
ggplot(aes(F,m)) +
geom_ribbon(aes(ymin=l,ymax=u),fill = 'gold') +
geom_line() +
geom_vline(xintercept = min(c(f.msy,f.pa))) +
geom_vline(xintercept = f.pa,lty=2,col='red') +
geom_vline(xintercept = f.lim,lwd=1.1,col='red') +
ssb_curve %>%
left_join(f.curve %>%
select(harvest_rate,F=m)) %>%
ggplot(aes(F,m)) +
geom_ribbon(aes(ymin=l,ymax=u),fill = 'gold') +
geom_line() +
geom_vline(xintercept = min(c(f.msy,f.pa))) +
geom_vline(xintercept = f.pa,lty=2,col='red') +
geom_vline(xintercept = f.lim,lwd=1.1,col='red') +
geom_hline(yintercept = blim, col = 'red', lwd = 1.1) +
geom_hline(yintercept = bpa,col = 'red',lty = 2)
## run prognosis
res <-
gadget_project_time(num_years = 5, variant_dir = 'PRG') %>%
gadget_project_stocks(imm.file = 'Modelfiles/cod.imm',mat.file = 'Modelfiles/cod.mat') %>%
gadget_project_fleet(pre_fleet = 'comm') %>%
gadget_evaluate(params.out = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
params.in = 'WGTS/params.final') %>%
gadget_project_recruitment(stock = 'codimm',
recruitment = fit$stock.recruitment %>%
filter(stock == 'codimm',
year > 1980),
method = 'constant',
n_replicates = 1,
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
gadget_project_ref_points(ref_points = tibble(codmat.blim = 207727665),
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
gadget_project_advice(pre_fleet = 'comm',
harvest_rate = hr_msy,
params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
n_replicates = 1,
advice_cv = 0) %>%
gadget_project_output(imm.file = 'Modelfiles/cod.imm',mat.file = 'Modelfiles/cod.mat',
pre_fleet = 'comm') %>%
gadget_evaluate(params.in = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
{read.printfiles(paste(attr(.,'variant_dir'),'out',sep='/'))} %>%
set_names(c("catch.F","catch",'rec','ssb'))
}