In einarhjorleifsson/fishvice: A collection of various personal functions used in my work on fish stock assessment

knitr::opts_chunk$set(
  message = FALSE, 
  warning = FALSE,
  collapse = TRUE,
  error = TRUE,
  comment = "#>"
)

library(FLCore)
library(FLXSA)
library(ggplotFL)
library(patchwork)
library(tidyverse)
library(fishvice)

Convert rbx to flr

rbx <- fishvice::mup_rbx("/u2/reikn/Tac/2022/01-cod/ass/mup/smx", scale = 1000)
rbx$rby <- rbx$rby %>% filter(year %in% 1985:2022)
rbx$rbya <- rbx$rbya %>% filter(year %in% 1985:2022)
# set maturity such that SSB is actually B4+
if(TRUE) {
  rbx$rbya <-
    rbx$rbya %>% 
    mutate(mat = replace_na(mat, 0),
           mat = ifelse(age >= 4, 1, 0),
           ssbW = cW,
           ssbW = replace_na(ssbW, 0))
}


fls <- rbx %>% rbx_to_flstock(pf = 0, pm = 0, fage = c(5, 10))
fls %>% plot()

Run a default xsa

... i.e. no fine tweeking of the default xsa setup

years <- 1985:2021
fls <- fls %>% trim(year = years)

ind.smb <- rbx_to_flindex(rbx, "U1") %>% trim(year = years)
ind.smh <- rbx_to_flindex(rbx, "U2", time = c(10/12, 10.5/12)) %>% trim(year = 1996:2021)
ind <- FLIndices(ind.smb)
cntr <- FLXSA.control(maxit = 70, rage = 5, qage = 14)
xsa   <- FLXSA(fls, ind, cntr)
#diagnostics(xsa)
fls.new  <- fls + xsa
plot(fls.new) +
  geom_vline(xintercept = 2021)

Compare muppet and xsa

... check terminal biomass for xsa, not the same as the plot above

fls.new %>% 
  fishvice::flstock_to_rbya() %>% 
  group_by(year) %>% 
  summarise(fbar = mean(f[age %in% 5:10]),
            bio = sum(n[age %in% 4:14] * cW[age %in% 4:14])) %>% 
  mutate(run = "xsa") %>% 
  bind_rows(rbx$rby %>% 
              filter(year %in% 1985:2022) %>% 
              select(year, fbar, bio) %>% 
              mutate(run = "mup")) %>% 
  gather(var, val, -c(year, run)) %>% 
  ggplot(aes(year, val, colour = run)) +
  theme_bw(base_size = 15) +
  geom_line(lwd = 1) +
  scale_colour_brewer(palette = "Set1") +
  facet_wrap(~ var, scales = "free_y")

einarhjorleifsson/fishvice documentation built on Jan. 4, 2024, 8:43 p.m.