---
title: "Operating model summary"
subtitle: "Comparison of OM historical period and Stock Synthesis output"
date: "`r Sys.Date()`"
---
<style type="text/css">
h1 { /* Header 1 */
  font-size: 24px;
}
</style>

```{r setup, include = FALSE, echo = FALSE}
  knitr::opts_chunk$set(collapse = TRUE, echo = FALSE, message = FALSE,
  fig.width = 6, fig.height = 4.5, out.width = "650px", comment = "#>")
```

# {.tabset}

## Spawning biomass

```{r, fig.cap = "Spawning biomass from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
if(season_as_years) {
  OM_years <- seas1_yind_full$assess_year[seas1_yind]
  year_lab <- "Stock Synthesis Year"
} else {
  OM_years <- mainyrs
  year_lab <- "Year"
}
ylim_SSB <- c(0, 1.1 * max(c(SSB, Hist@TSdata$SSB)))
matplot(OM_years, t(Hist@TSdata$SSB), xlab = year_lab, ylab = "Spawning Biomass", ylim = ylim_SSB, pch = 1, col = "black", typ = "o")
lines(mainyrs, SSB, col = "red", lwd = 3)
abline(h = 0, col = "grey")
```

```{r, fig.cap = "Spawning depletion from the operating model (black lines and dots) and Stock Synthesis (red lines).."}
ylim_dep <- c(0, 1.1 * max(SSB/SSB0, Hist@TSdata$SSB/Hist@Ref$SSB0))
matplot(OM_years, t(Hist@TSdata$SSB/Hist@Ref$SSB0), xlab = year_lab, ylab = "Spawning Depletion", ylim = ylim_dep, pch = 1, col = "black", typ = "o")
lines(mainyrs, SSB/SSB0, col = "red", lwd = 3)
abline(h = 0, col = "grey")
```

## Catch

```{r, fig.cap = "Catch from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
Catch_SS <- cbind(replist$timeseries[, 1:4], rowSums(replist$timeseries[, startsWith(names(replist$timeseries), "dead(B):")], na.rm = TRUE))
Catch_SS <- aggregate(Catch_SS[, 5], list(Year = Catch_SS$Yr), sum)
Catch_SS <- Catch_SS[vapply(Catch_SS$Year, "%in%", logical(1), mainyrs), 2]
if(season_as_years) Catch_SS <- aggregate(Catch_SS, list(Year = seas1_yind_full$true_year), sum)[, 2]

ylim_cat <- c(0, 1.1 * max(Catch_SS, Hist@TSdata$Catch))
matplot(OM_years, t(Hist@TSdata$Catch), xlab = year_lab, ylab = "Catch", ylim = ylim_cat, pch = 1, col = "black", typ = "o")
lines(OM_years, Catch_SS, col = "red", lwd = 3)
abline(h = 0, col = "grey")
```

```{r, fig.cap = "Relative catch from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
ylim_cat <- c(0, 1.1 * max(Catch_SS/Catch_SS[length(Catch_SS)], Hist@TSdata$Catch/Hist@TSdata$Catch[, OM@nyears]))
matplot(OM_years, t(Hist@TSdata$Catch/Hist@TSdata$Catch[, OM@nyears]), xlab = year_lab, ylab = "Catch Relative to today", ylim = ylim_cat, pch = 1, col = "black", typ = "o")
lines(OM_years, Catch_SS/Catch_SS[length(Catch_SS)], col = "red", lwd = 3)
abline(h = 0, col = "grey")
abline(h = 1, lty = 3)
```

## Abundance-at-age {.tabset}

### Recruitment

```{r, fig.cap = "Recruitment from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
N_at_age <- reshape2::melt(replist$natage[replist$natage$`Beg/Mid` == "B" & replist$natage$Seas == 1, -c(9:11)], list("Area", "Bio_Pattern", "Sex", "BirthSeas", "Platoon", "Morph", "Seas", "Yr"),
                           variable.name = "Age", value.name = "N")
N_at_age <- N_at_age[as.numeric(N_at_age$Age)-1 >= age_rec * ifelse(season_as_years, nseas, 1), ]  # Subset by age >= age_rec
N_at_age <- N_at_age[vapply(N_at_age$Yr, "%in%", logical(1), table = mainyrs), ]  # Subset by year
N_at_age <- N_at_age[vapply(N_at_age$Sex, "%in%", logical(1), table = gender), ] # Subset by gender

if(season_as_years) { # Sum across sub-ages, select years corresponding to true season 1
  N_at_age$true_Age <- seas1_aind_full$true_age[match(N_at_age$Age, seas1_aind_full$assess_age)]
  N_at_age <- summarise(group_by(N_at_age, Sex, Yr, true_Age), N = sum(N))
  
  N_at_age <- N_at_age[vapply(N_at_age$Yr, "%in%", logical(1), seas1_yind_full$assess_year[seas1_yind_full$nseas == 1]), ]
  N_at_age <- summarise(group_by(N_at_age, Yr, true_Age), N = sum(N))
  
  N_at_age <- reshape2::acast(N_at_age, list("Yr", "true_Age"), value.var = "N")
} else {
  N_at_age <- summarise(group_by(N_at_age, Yr, Age), N = sum(N)) # Sum over area, morphs, platoons, sex, etc.
  N_at_age <- reshape2::acast(N_at_age, list("Yr", "Age"), value.var = "N")
}

ylim_rec <- c(0, 1.1 * max(Hist@AtAge$Nage[, 1, ], N_at_age[, 1]))
matplot(OM_years, t(Hist@AtAge$Nage[, 1, ]), xlab = year_lab, ylab = "Recruitment (age 1)", ylim = ylim_rec, pch = 1, col = "black", typ = "o")
lines(OM_years, N_at_age[, 1], col = "red", lwd = 3)
abline(h = 0, col = "grey")
```

### Absolute abundance

```{r, fig.cap = "Annual abundance at age from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
plot_composition(OM_years, t(Hist@AtAge$Nage[1, , 1:nyears]), N_at_age, N = NULL, annual_ylab = "Numbers at age", annual_yscale = "raw")
```

### Relative abundance

```{r, fig.cap = "Annual proportions at age from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
plot_composition(OM_years, t(Hist@AtAge$Nage[1, , 1:nyears]), N_at_age, N = NULL, annual_ylab = "Proportions at age")
```

## Selectivity and F {.tabset}

### Apical F

```{r, fig.cap = "Apical F from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
Hist_F <- apply(Hist@AtAge$FM, c(1, 3), max)
ylim_F <- c(0, 1.1 * max(c(OM@cpars$Find, Hist_F)))
matplot(OM_years, t(Hist_F), typ = "o", xlab = "Year", ylab = "Apical F", ylim = ylim_F, pch = 1, col = "black")
lines(OM_years, OM@cpars$Find[1, ], col = "red", lwd = 3)
abline(h = 0, col = "grey")
```

### Selectivity 

```{r, fig.cap = "Annual selectivity at age from the operating model (black lines and dots) and Stock Synthesis (red lines)."}
plot_composition(OM_years, t(Hist@AtAge$Select[1, , 1:nyears]), t(OM@cpars$V[1, , 1:nyears]), N = NULL, annual_ylab = "Selectivity", annual_yscale = "raw")
```

## About

This report was generated on: `r Sys.time()`<br />
MSEtool version `r packageVersion("MSEtool")`<br />
DLMtool version `r packageVersion("DLMtool")`<br />
`r R.version.string`<br />

