##' \code{gadget.fit} calls Gadget and compiles fit statistics, relevant biomass trajectories and consumption.
##' It returns an object of class 'gadget.fit' which is essentially a list containing the following objects: \describe{
##' \item{sidat}{data.frame containing the model fit to the surveyindices}
##' \item{resTable}{Table of likelihood component scores from the different stages of the iterative reweighting run}
##' \item{nesTable}{Same as resTable but normalised with the minimum score for each component}
##' \item{suitability}{The model estimated fleet predator suitability}
##' \item{stock.recruitment}{Recruitment by stock}
##' \item{res.by.year}{Results by year}
##' \item{stomachconten}{Fit combined with data from the stomachcontent likelihood}
##' \item{likelihoodsummary}{Likelihood component scores by year and step}
##' \item{catchdist.fleets}{Data compared with model output for the catchdistribution components}
##' \item{stockdist}{Data compared with model output for the stockdistribution components}
##' \item{SS}{The likelihood output from the model}
##' \item{stock.full}{Modeled abundance and mean weight by year,step, length and stock}
##' \item{stock.std}{Modeled abundance, mean weight, number consumed by age, stock and year}
##' \item{stock.prey}{Consumption of prey (or fleet catches) by prey, length, year and step}
##' \item{fleet.info}{Information on catches, harvest rate and harvestable biomass by fleet, year and step}
##' \item{predator.prey}{consumption of preys by predators, length, year and step}
##' \item{params}{parameter values used for the fit}
##' \item{catchstatistics}{Data compared with model output for the catchdstatistics components}
##' }
##' @title Gadget fit
##' @rdname gadget.fit
##' @param wgts Location of the iterative reweighting folder, if null gadget.fit requires a parameter file
##' @param main.file if the main file is different from the 'main' or 'wgts/main.final'
##' @param fleet.predict data.frame containing the basis fleets used to calculate the harvestable biomass. Fleet names should be specified in the fleet column.
##' @param mat.par parameters for the maturity ogive
##' @param params.file parameter file used in the fit (defaults to "WGTS/params.final")
##' @param fit.folder location of the output
##' @param printatstart should the stock standard output be printed at the beginning or the end of the timestep
##' @param steps what steps should be printed
##' @param f.age.range data.frame describing the desired age range where the F's are calculated, if null this defaults to the apical F fro all stocks.
##' Input columns should include stock, age.min and age.max
##' @param recruitment_step_age data frame defining the recruitment age by stock and time step. Default is the minimum age at step 1. Expects columns stock, age and step
##' @param gd gadget_directory object to optionally set the location of the gadget directory
##' @return list containing the output from Gadget.
##' @importFrom rlang .data
##' @author Bjarki Thor Elvarsson
##' @export
gadget.fit <- function(wgts = 'WGTS',
main.file = NULL,
fleet.predict = NULL,
mat.par=NULL,
params.file=NULL,
f.age.range=NULL,
fit.folder = 'FIT',
printatstart = 1,
steps = 1,
recruitment_step_age = NULL,
gd = NULL){
warning('gadget.fit is deprecated, please use gadget_fit')
old.dir <- getwd()
on.exit(setwd(old.dir))
if(!is.null(gd)){
setwd(gd)
} else {
gd <- gadget.variant.dir('.')
}
if(!is.null(f.age.range) & !('data.frame' %in% class(f.age.range))){
stop('F age range should be specified as a data.frame with columns stock, age.min and age.max')
}
if(is.null(main.file)) {
if(is.null(wgts)){
main.file <- 'main'
} else {
main.file <- sprintf('%s/main.final',wgts)
}
}
main <- read.gadget.main(main.file)
if(!is.null(wgts)){
resTable <- read.gadget.results(wgts=wgts)
nesTable <- read.gadget.results(wgts=wgts,normalize = TRUE)
if(is.null(params.file)){
params.file <- sprintf('%s/params.final',wgts)
}
params <- read.gadget.parameters(params.file)
lik <- read.gadget.likelihood(sprintf('%s/likelihood.final',wgts))
} else {
resTable <- list()
nesTable <- list()
wgts <- fit.folder
dir.create(fit.folder,showWarnings = FALSE)
params <- if(!is.null(params.file)) read.gadget.parameters(params.file) else NULL
lik <- read.gadget.likelihood(main$likelihoodfiles)
}
print('Reading input data')
lik.dat <- read.gadget.data(lik)
stocks <-
main$stockfiles %>%
purrr::map(~read.gadget.file(path='.',file_name = .,file_type = 'stock',recursive = FALSE))
names(stocks) <- stocks %>% purrr::map(1) %>% purrr::map('stockname') %>% unlist()
## model output, i.e printfiles
make.gadget.printfile(main.file = main.file,
file = sprintf('%s/printfile.fit',wgts),
gd = list(dir='.',output = sprintf('%s/out.fit',wgts),
aggfiles = sprintf('%s/print.aggfiles',wgts)),
recruitment_step_age = recruitment_step_age,
printatstart = printatstart,
steps = steps)
on.exit(unlink(sprintf('%s/out.fit',wgts),recursive = TRUE))
main$printfiles <- sprintf('%s/printfile.fit',wgts)
write.gadget.main(main,file = sprintf('%s/main.print',wgts))
print("Running Gadget")
callGadget(s=1,
i = params.file,
main = sprintf('%s/main.print',wgts),
o = sprintf('%s/SS.print',wgts))
print("Reading output files")
out <- read.printfiles(sprintf('%s/out.fit',wgts)) %>%
setNames(gsub("^[^a-zA-Z]*","",names(.)))
SS <- tryCatch(read.gadget.lik.out(sprintf('%s/SS.print',wgts)),
error = function(e) 'SS could not be read')
#stocks <- read.gadget.stockfiles(main$stockfiles)
print('Gathering results')
stock.recruitment <-
out[sprintf('%s.recruitment',names(stocks))] %>%
purrr::set_names(.,names(stocks)) %>%
dplyr::bind_rows(.id='stock') %>%
dplyr::select(.data$stock,.data$year,.data$area,.data$step, recruitment=.data$number)
stock.full <-
out[sprintf('%s.full',names(stocks))] %>%
purrr::set_names(.,names(stocks)) %>%
dplyr::bind_rows(.id='stock') %>%
dplyr::mutate(length=as.numeric(gsub('len','',.data$length)))
stock.std <-
out[sprintf('%s.std',names(stocks))] %>%
purrr::set_names(.,names(stocks)) %>%
dplyr::bind_rows(.id='stock')
if(sum(grepl('prey$',names(out)))>0){
stock.prey <-
out[sprintf('%s.prey',names(stocks))] %>%
purrr::set_names(.,names(stocks)) %>%
dplyr::bind_rows(.id='stock')
} else {
stock.prey <-
tibble::tibble(year=NA_real_, step=NA_real_, area=NA_character_, stock = NA_character_,
age=NA_real_,biomass_consumed=NA_real_,number_consumed=NA_real_, mortality = NA_real_)
}
if(sum(grepl('prey',names(out)))>0){
predator.prey <-
out[grepl('.+\\.prey\\..+',names(out))] %>%
purrr::set_names(.,names(.)) %>%
purrr::keep(~'number_consumed' %in% names(.)) %>%
dplyr::bind_rows(.id='stock') %>%
tidyr::separate(.data$stock,c('prey','predator'),sep='\\.prey\\.') %>%
dplyr::group_by(.data$year,.data$step,.data$prey,.data$predator) %>%
dplyr::mutate(suit = .data$mortality/max(.data$mortality),
suit = ifelse(is.finite(.data$suit),.data$suit,0),
length = gsub('len', '', .data$length) %>%
as.numeric())
} else {
predator.prey <-
tibble::tibble(year=NA_real_, step=NA_real_, area=NA_character_,
predator = NA_character_,
prey = NA_character_,
length = NA_real_,
suit = NA_real_,
biomass_consumed = NA_real_, mortality = NA_real_)
}
fleet.catches <-
predator.prey %>%
dplyr::group_by(.data$year,.data$area,.data$predator,.data$prey) %>%
dplyr::summarise(amount = sum(.data$biomass_consumed)) %>%
dplyr::rename(fleet = .data$predator, stock = .data$prey)
fleet.info <-
stock.full %>%
dplyr::left_join(predator.prey %>%
dplyr::select(.data$year,
.data$step,
.data$area,
fleet = .data$predator,
stock = .data$prey,
.data$length,
.data$suit),
by = c("stock", "year", "step", "area", "length")) %>%
dplyr::group_by(.data$year,.data$step,.data$area,.data$fleet) %>%
dplyr::summarise(harv.bio = sum(.data$suit*.data$number*.data$mean_weight)) %>%
dplyr::left_join(fleet.catches %>%
dplyr::group_by(.data$year,.data$fleet,.data$area) %>%
dplyr::summarise(amount=sum(.data$amount)),
by = c("year", "area", "fleet")) %>%
dplyr::group_by(.data$year,.data$step,.data$area,.data$fleet) %>%
dplyr::mutate(amount = ifelse(is.na(.data$amount),0,.data$amount),
harv.rate = .data$amount/.data$harv.bio)
if(!is.null(fleet.predict)){
d <-
predator.prey %>%
dplyr::filter(.data$predator %in% fleet.predict$fleet)
} else {
d <- predator.prey
}
harv.suit <-
d %>%
dplyr::group_by(.data$year,.data$step,.data$prey,.data$length) %>%
dplyr::filter(.data$biomass_consumed > 0) %>%
dplyr::summarise(suit = sum(.data$biomass_consumed*.data$suit)/sum(.data$biomass_consumed)) %>%
dplyr::rename(stock = .data$prey)
print('Merging input and output')
## merge data and estimates
if('surveyindices' %in% names(lik.dat$dat)){
sidat <-
out[names(lik.dat$dat$surveyindices)] %>%
purrr::set_names(.,names(.)) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::left_join(lik$surveyindices %>%
dplyr::select(.data$name,.data$stocknames,.data$sitype,.data$fittype),
by='name') %>%
dplyr::bind_rows(tibble::tibble(length=NA,
age=NA,
survey = NA,
fleet = NA)) %>%
dplyr::mutate(age = ifelse(.data$sitype == 'ages',.data$label,.data$age),
length = ifelse(.data$sitype %in% c('lengths','fleets'),.data$label,.data$length),
fleet = ifelse(.data$sitype == 'effort',.data$label,.data$fleet),
survey = ifelse(.data$sitype == 'acoustic',.data$label,.data$survey)) %>%
dplyr::left_join(lik.dat$dat$surveyindices %>%
purrr::set_names(.,names(.)) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::as_tibble() %>%
dplyr::rename(observed=.data$number) %>%
dplyr::bind_rows(tibble::tibble(name = NA, year = NA, step = NA,
area = NA,length = NA,age = NA,
fleet = NA,survey = NA,
upper = NA, lower = NA)) %>%
dplyr::filter(!is.na(.data$year)),
by = c("name", "year", "step", "area", "age", "length","fleet","survey")) %>%
dplyr::mutate(length = ifelse(.data$sitype %in% c('lengths','fleets'),
paste(.data$lower,.data$upper,sep=' - '),
.data$length)) %>%
dplyr::mutate(predicted = ifelse(grepl('loglinearfit',tolower(.data$fittype)),
exp(.data$intercept)*.data$number^.data$slope,
.data$intercept + .data$slope*.data$number)) %>%
dplyr::filter(!is.na(.data$name))
} else {
sidat <- NULL
}
if('catchdistribution' %in% names(lik.dat$dat)){
dat.names <- names(lik.dat$dat$catchdistribution)
aggs <-
dat.names %>%
purrr::set_names(.,.) %>%
purrr::map(~attr(lik.dat$dat$catchdistribution[[.]],'len.agg')) %>%
dplyr::bind_rows(.id='name') %>%
tibble::as_tibble()
catchdist.fleets <-
lik.dat$dat$catchdistribution %>%
purrr::set_names(.,names(.)) %>%
purrr::map(. %>% dplyr::mutate(age = as.character(.data$age))) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::right_join(out[dat.names] %>%
purrr::set_names(.,dat.names) %>%
purrr::map(. %>% dplyr::mutate(age = as.character(.data$age))) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::left_join(aggs,by=c('name','length')) ,
by=c('name','length', 'year',
'step', 'area','age','upper','lower')) %>%
dplyr::ungroup() %>%
dplyr::group_by(.data$name,.data$year, .data$step, .data$area) %>%
dplyr::mutate(total.catch = sum(.data$number.x,na.rm=TRUE),
total.pred = sum(.data$number.y,na.rm=TRUE),
observed = .data$number.x/sum(.data$number.x,na.rm=TRUE),
predicted = .data$number.y/sum(.data$number.y,na.rm=TRUE)) %>%
dplyr::ungroup() %>%
dplyr::group_by(.data$name,.data$length,.data$age) %>%
dplyr::mutate(upper = as.double(max(ifelse(is.na(.data$upper),0.0,
.data$upper))),
lower = as.double(max(ifelse(is.na(.data$lower),0.0,
.data$lower))),
avg.length = as.numeric((.data$lower+.data$upper)/2),
residuals = as.numeric(.data$observed - .data$predicted)) %>%
dplyr::inner_join(lik$catchdistribution %>%
dplyr::select(.data$name,.data$fleetnames,.data$stocknames),
by = 'name')
} else {
catchdist.fleets <- NULL
}
if(sum(grepl('.std',names(out),fixed = TRUE))>0){
if(is.null(f.age.range)){
f.age.range <-
stock.prey %>%
dplyr::group_by(.data$stock) %>%
dplyr::summarise(age.min = max(.data$age),age.max=max(.data$age))
}
f.by.year <-
stock.prey %>%
dplyr::left_join(f.age.range,by="stock") %>%
dplyr::group_by(.data$stock,.data$year,.data$area) %>%
dplyr::summarise(catch=sum(.data$biomass_consumed),
num.catch=sum(.data$number_consumed),
F=mean(.data$mortality[.data$age>=.data$age.min&.data$age<=.data$age.max]))
res.by.year <-
stock.full %>%
dplyr::filter(.data$step %in% steps) %>%
dplyr::left_join(harv.suit,
by = c("stock", "year", "step", "length")) %>%
dplyr::group_by(.data$stock,.data$year,.data$area,.data$step) %>%
dplyr::summarise(total.number = sum(.data$number),
total.biomass = sum(.data$number*.data$mean_weight),
harv.biomass = sum(.data$number*.data$suit*.data$mean_weight),
ssb = sum(.data$mean_weight*logit(mat.par[1],
mat.par[2],
.data$length)*
.data$number)) %>%
dplyr::left_join(f.by.year %>%
dplyr::mutate(area = .data$area),
by = c("stock","year","area")) %>%
dplyr::left_join(stock.recruitment %>%
dplyr::group_by(.data$stock,.data$year,.data$area) %>%
dplyr::summarise(recruitment = sum(.data$recruitment)),
by = c('stock','year','area')) %>%
dplyr::ungroup()
} else {
res.by.year <- NULL
}
# if('catchinkilos' %in% names(lik.dat$dat)){
# catchdist <-
# ldply(names(lik.dat$dat$catchinkilos),
# function(x){
#
# }
# }
if('stockdistribution' %in% names(lik.dat$dat)){
dat.names <- names(lik.dat$dat$stockdistribution)
aggs <-
dat.names %>%
purrr::set_names(.,.) %>%
purrr::map(~attr(lik.dat$dat$stockdistribution[[.]],'len.agg')) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::as_tibble()
stockdist <-
lik.dat$dat$stockdistribution %>%
purrr::set_names(.,names(.)) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::right_join(out[dat.names] %>%
purrr::set_names(.,dat.names) %>%
dplyr::bind_rows(.id='name') %>%
dplyr::left_join(aggs,by=c('name','length')),
by=c('name','length', 'year',
'step', 'area','age',
'stock','upper','lower'),
suffix = c('.y','.x')) %>%
dplyr::group_by(.data$name, .data$year, .data$step,
.data$area, .data$age, .data$length) %>%
dplyr::mutate(pred.ratio = .data$number.x/sum(.data$number.x,na.rm=TRUE),
obs.ratio = .data$number.y/sum(.data$number.y)) %>%
dplyr::ungroup() %>%
dplyr::mutate(length = (.data$lower+.data$upper)/2) %>%
dplyr::inner_join(lik$stockdistribution %>%
dplyr::select(.data$name,.data$fleetnames,.data$stocknames) %>%
dplyr::distinct(),
by='name')
} else {
stockdist <- NULL
}
if('stomachcontent' %in% names(lik.dat$dat)){
pred.agg <-
names(lik.dat$dat$stomachcontent) %>%
purrr::set_names(.,.) %>%
purrr::map(~attr(lik.dat$dat$stomachcontent[[.]],'pred.agg')) %>%
dplyr::bind_rows(.id='component')
prey.agg <-
names(lik.dat$dat$stomachcontent) %>%
purrr::set_names(.,.) %>%
purrr::map(~attr(lik.dat$dat$stomachcontent[[.]],'prey.agg')) %>%
dplyr::bind_rows(.id='component')
stomachcontent <-
lik.dat$dat$stomachcontent %>%
dplyr::bind_rows(.id = 'component') %>%
dplyr::select(-c(.data$prey.lower,.data$prey.upper,
.data$lower,.data$upper)) %>%
dplyr::right_join(
out[lik.dat$dat$stomachcontent %>%
names()] %>%
dplyr::bind_rows(.id = 'component') %>%
dplyr::left_join(prey.agg) %>%
dplyr::left_join(pred.agg),
by = c('component','predator','prey','year','step','area')) %>%
dplyr::group_by(.data$component,.data$year,.data$step,.data$predator) %>%
dplyr::mutate(observed=.data$ratio/sum(.data$ratio,na.rm=TRUE),
predicted=.data$number/sum(.data$number,na.rm=TRUE),
prey.length = (.data$prey.lower+.data$prey.upper)/2,
pred.length = (.data$lower+.data$upper)/2) %>%
dplyr::as_tibble()
} else {
stomachcontent <- NULL
}
if('catchstatistics' %in% names(lik.dat$dat)){
catchstatistics <-
lik.dat$dat$catchstatistics %>%
names() %>%
purrr::set_names(.,.) %>%
purrr::map(function(x){
out[[x]] %>%
dplyr::rename(fitted_mean = .data$mean,
fitted_number = .data$number) %>%
dplyr::left_join(lik.dat$dat$catchstatistics[[x]],
by = c("year", "step", "area", "age")) %>%
dplyr::rename(observed_mean = .data$mean,
observed_number = .data$number)
}) %>%
dplyr::bind_rows(.id = 'name') %>%
dplyr::as_tibble()
} else {
catchstatistics <- NULL
}
out <-
list(sidat = sidat, resTable = resTable, nesTable = nesTable,
suitability = predator.prey %>%
dplyr::select(.data$year,.data$step,stock=.data$prey,
fleet=.data$predator,.data$length,.data$suit),# gss.suit,
#stock.growth = stock.growth,
stock.recruitment = stock.recruitment,
res.by.year = res.by.year,
stomachcontent = stomachcontent,
likelihoodsummary = out$likelihoodsummary %>% dplyr::as_tibble(),
catchdist.fleets = catchdist.fleets,
stockdist = stockdist,
#out.fit=out,
SS = SS,
stock.full = stock.full,
stock.std = stock.std,
stock.prey = stock.prey,
fleet.info = fleet.info,
predator.prey = predator.prey,
params = params,
catchstatistics = catchstatistics)
class(out) <- c('gadget.fit',class(out))
save(out,file=sprintf('%s/WGTS.Rdata',wgts))
setwd(old.dir)
return(out)
}
#' @param ... any number of gadget.fit objects
#'
#' @rdname gadget.fit
#' @description \code{bind.gadget.fit} combines two (or more) gadget.fit objects into a
#' single list. Each object in the list will have an extra column indicating the model,
#' if the list is named those names will be used to identify the model.
#' @export
bind.gadget.fit <- function(...){
tmp <-
list(...) %>%
purrr::transpose() %>%
purrr::map(~purrr::map_if(.,is.null,data.frame)) %>%
purrr::map(~purrr::map_if(.,is.factor,as.character)) %>%
purrr::map(purrr::safely(~dplyr::bind_rows(.,.id='model'))) %>%
purrr::map('result')
class(tmp) <- c('gadget.fit',class(tmp))
return(tmp)
}
get_gadget_recruitment <- function(stocks,params){
stocks %>%
purrr::map(function(x){
if(x$doesrenew$doesrenew){
if('normalparamfile' %in% names(x$doesrenew)){
key <- 'normalparamfile'
col.names <- c('year','step','area','age','number','mean','stddev','alpha','beta')
} else if('normalcondfile' %in% names(x$doesrenew)){
key <- 'normalcondfile'
col.names <- c('year','step','area','age','number','mean','stddev','relcond')
} else {
key <- 'numberfile'
col.names <- c('year','step','area','age','length','number','weight')
}
utils::capture.output(x$doesrenew[[key]] %>% print()) %>%
utils::read.table(text = ., comment.char = ';',sep = '\t',fill = TRUE,stringsAsFactors = FALSE) %>%
stats::na.omit() %>%
purrr::set_names(.,col.names)
} else{
NULL
}
}) %>%
dplyr::bind_rows(.id = 'stock') %>%
dplyr::mutate_at(.vars=dplyr::vars(-.data$stock),
~purrr::map(.,function(x)
tryCatch(parse.gadget.formulae(x) %>%
eval(params %>%
dplyr::select(.data$value) %>%
t() %>%
purrr::set_names(.,dimnames(.)[[2]]) %>%
as.list()),
error = function(e) 0)))
}
#' Retro fit
#'
#' Compile the fit objects from retrospective runs
#' @param pre Folder containing the retro run
#' @return gadget fit list
#' @export
gadget.retro.fit <- function(pre='RETRO'){
tmp_func <- purrr::lift(bind.gadget.fit,.unnamed = TRUE)
list.files(pre,pattern = 'R[0-9]+') %>%
purrr::set_names(.,.) %>%
purrr::map(function(x){
print(x)
gadget.fit(main.file = sprintf('%s/%s/main',pre,x),
params.file = sprintf('%s/params.retro.%s',pre,gsub('R','',x)),
wgts = NULL)
}) %>%
tmp_func()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.