gadget_projections: Gadget projection function

In gadget-framework/rgadget: Tools to Create and Maintain Gadget Models

Gadget projection function

Description

Setup model and parameter files for forward simulations and deterministic projections.

Usage

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"
)

Arguments

path	gadget variant director
num_years	number of years
variant_dir	location of the
imm.file	location of the immature stock file
mat.file	location of the mature stock file
spawnfunction	what spawn function to use (only hockeystick atm)
proportionfunction	proportion suitability
mortalityfunction	mortataliy suitability
weightlossfunction	weightloss suitability
pre_fleet	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.
post_fix	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.
fleet_type	type of gadget fleet for the projections
common_mult	a string with a base name for a vector of harvest/catch rate multipliers (shared between fleets)
pre_proportion	proportion of the effort taken
type	type of projection, either "standard" or "prognosis"
...	additional arguments for the
stocks	names of of the stocks
pre_fleets	vector of fleets on which the projections are based
mainfile_overwrite	should the likelihood be overwritten, defaults to TRUE, but should be set as FALSE for multiple prognosis components (different groups of fleets)
firsttacyear	when does the TAC scheme start
assessmentstep	when is the assessment
weightoflastyearstac	running average TAC..
stock	name of immature stock
recruitment	recruitment time series
n_replicates	number of simulations
params.file	name of the parameter files
method	prediction method, built in options are 'AR', 'bootstrap' or constant
harvest_rate	median harvest rate
advice_cv	assessment error cv
advice_rho	assessment error correlation
ref_points	tibble with reference points
f_age_range	F age range, specified in the format a1:a2
rec_age	age of recruitment (as reported)
print_block	file suffix for the printfile

Details

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.

Value

gadget variant directory

Examples

## Not run: 

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')) 
 
 
 
## End(Not run)
 
 

gadget-framework/rgadget documentation built on July 20, 2022, 12:16 p.m.