R/rmd_SMSE.R

In salmonMSE: Management Strategy Evaluation for Salmon Species

rmd_natural_production <- function(s, sname) {

  c(
    paste("##", sname, "{.tabset}"),
    "\n",
    "### Natural production {.tabset}",
    "\n",
    "#### Parameters",
    "\n",
    "```{r, fig.cap=\"Histogram of kappa (productivity, recruits/spawner).\"}",
    "SOM <- SMSE@Misc$SOM",
    paste0("use_habitat <- SOM@Habitat[[", s, "]]@use_habitat"),
    paste0("SRrel <- SOM@Bio[[", s, "]]@SRrel"),
    "",
    paste0("if (!use_habitat) plot_SOM(SOM@Bio[[", s, "]], var = \"kappa\", xlab = \"Productivity\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Histogram of Smax (spawner abundance that maximizes smolt production).\"}",
    paste0("if (!use_habitat && SRrel == \"Ricker\") plot_SOM(SOM@Bio[[", s, "]], var = \"Smax\", xlab = expression(S[max]))"),
    "```",
    "\n",
    "```{r, fig.cap=\"Histogram of smolt capacity from egg production (Beverton-Holt asymptote).\"}",
    paste0("if (!use_habitat && SRrel == \"BH\") plot_SOM(SOM@Bio[[", s, "]], var = \"capacity\", xlab = \"Smolt capacity\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Smolt production (density-dependent stock-recruit relationship). The dotted vertical line denotes the median egg production at unfished one to one replacement.\"}",
    paste0("if (!use_habitat) plot_SRR(SOM, s = ", s, ", surv = FALSE, check = FALSE, quant = TRUE)"),
    "```",
    "\n",
    "```{r, fig.cap=\"Egg to smolt survival from the density-dependent stock-recruit relationship. The dotted vertical line denotes the median egg production at unfished replacement.\"}",
    paste0("if (!use_habitat) plot_SRR(SOM, s = ", s, ", surv = TRUE, check = FALSE, quant = TRUE)"),
    "```",
    "```{r, fig.cap=\"Juvenile instantaneous natural mortality rate of natural origin fish (median and 95 percent interval across simulations).\"}",
    paste0("if (SOM@Bio[[", s, "]]@n_g > 1) {"),
    paste0("  plot_Mjuv_LHG(SOM@Bio[[", s, "]]@Mjuv_NOS[, , SOM@proyears, ])"),
    "} else {",
    paste0("  plot_SOM(SOM@Bio[[", s, "]], \"Mjuv_NOS\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage-1, proyears = SOM@proyears, ylab = \"Juvenile natural mortality rate\")"),
    "}",
    "```",
    "",
    "```{r, fig.cap=\"Juvenile natural survival of natural origin fish (median and 95 percent interval across simulations).\"}",
    paste0("if (SOM@Bio[[", s, "]]@n_g > 1) {"),
    paste0("  plot_Mjuv_LHG(SOM@Bio[[", s, "]]@Mjuv_NOS[, , SOM@proyears, ], surv = TRUE, ylab = \"Juvenile survival\")"),
    "} else {",
    paste0("  plot_SOM(SOM@Bio[[", s, "]], \"Mjuv_NOS\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage-1, proyears = SOM@proyears, surv = TRUE, ylab = \"Juvenile natural survival\")"),
    "}",
    "```",
    "",
    "```{r, fig.cap=\"Proportion mature at age (median and 95 percent interval across simulations).\"}",
    paste0("plot_SOM(SOM@Bio[[", s, "]], \"p_mature\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage, proyears = SOM@proyears, ylab = \"Proportion mature\")"),
    "```",
    "",
    "```{r, fig.cap=\"Fecundity at age (median and 95 percent interval across simulations).\"}",
    paste0("plot_SOM(SOM@Bio[[", s, "]], \"fec\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage, proyears = SOM@proyears, ylab = \"Fecundity\")"),
    "```",
    "\n",
    "#### Time series",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin spawners time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"NOS\", s = ", s, ", quant = TRUE)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin spawners, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"NOS\", s = ", s, ")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery origin spawners time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"HOS\", s = ", s, ", quant = TRUE)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery origin spawners, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"HOS\", s = ", s, ")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of hatchery spawners (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"pHOS_census\", s = ", s, ", quant = TRUE, ylab = \"P(HOS)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of hatchery spawners, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"pHOS_census\", s = ", s, ", xlab = \"P(HOS)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of wild spawners (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"p_wild\", s = ", s, ", quant = TRUE, ylab = \"P(WILD)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of wild spawners, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"p_wild\", s = ", s, ", xlab = \"P(WILD)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, natural origin (annual median and 95 percent interval). Losses include en-route mortality and broodtake.\"}",
    paste0("plot_statevar_ts(SMSE, var = \"pNOSesc\", s = ", s, ", quant = TRUE, ylab = \"Spawner/escapement (natural)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, natural origin, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"pNOSesc\", s = ", s, ", xlab = \"Spawner/escapement (natural)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, hatchery origin (annual median and 95 percent interval). Losses include en-route mortality, broodtake, and removal from spawning grounds.\"}",
    paste0("plot_statevar_ts(SMSE, var = \"pHOSesc\", s = ", s, ", quant = TRUE, ylab = \"Spawner/escapement (hatchery)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, hatchery origin, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"pHOSesc\", s = ", s, ", xlab = \"Spawner/escapement (hatchery)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Smolt production (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"Smolt\", s = ", s, ", quant = TRUE, ylab = \"Smolt production\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Smolt production, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"Smolt\", s = ", s, ", xlab = \"Smolt production\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Egg to smolt survival (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"ESS\", s = ", s, ", quant = TRUE, ylab = \"Egg-smolt survival\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Egg to smolt survival, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"ESS\", s = ", s, ", xlab = \"Egg-smolt survival\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    ""
  )

}

rmd_hatchery_production <- function(s, sname) {
  c(
    "### Hatchery production {.tabset}",
    "",
    "#### Parameters",
    "",
    "```{r results=\"asis\"}",
    paste0("knitr::kable(make_hatchery_table(SMSE, s = ", s, "), row.names = FALSE, caption = \"Hatchery enhancement settings.\")"),
    "```",
    "",
    "```{r results=\"asis\"}",
    paste0("knitr::kable(make_fitness_table(SMSE, s = ", s, "), row.names = FALSE, caption = \"Fitness parameters.\")"),
    "```",
    "",
    "```{r, fig.cap=\"Juvenile instantaneous natural mortality rate of hatchery origin fish (median and 95 percent interval across simulations).\"}",
    paste0("if (SOM@Hatchery[[", s, "]]@n_r > 1) {"),
    paste0("  plot_Mjuv_RS(SOM@Hatchery[[", s, "]]@Mjuv_HOS[, , SOM@proyears, ])"),
    "} else {",
    paste0("  plot_SOM(SOM@Hatchery[[", s, "]], \"Mjuv_HOS\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage-1, proyears = SOM@proyears, ylab = \"Juvenile natural mortality rate\")"),
    "}",
    "```",
    "",
    "```{r, fig.cap=\"Juvenile natural survival of hatchery origin fish (median and 95 percent interval across simulations).\"}",
    paste0("if (SOM@Hatchery[[", s, "]]@n_r > 1) {"),
    paste0("  plot_Mjuv_RS(SOM@Hatchery[[", s, "]]@Mjuv_HOS[, , SOM@proyears, ], surv = TRUE)"),
    "} else {",
    paste0("  plot_SOM(SOM@Hatchery[[", s, "]], \"Mjuv_HOS\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage-1, proyears = SOM@proyears, surv = TRUE, ylab = \"Juvenile natural survival\")"),
    "}",
    "```",
    "",

    "```{r, fig.cap=\"Maturity of hatchery origin fish (median and 95 percent interval across simulations).\"}",
    paste0("if (SOM@Hatchery[[", s, "]]@n_r > 1) {"),
    paste0("  plot_Mjuv_RS(SOM@Hatchery[[", s, "]]@p_mature_HOS[, , SOM@proyears, ], ylab = \"Proportion mature\")"),
    "} else {",
    paste0("  plot_SOM(SOM@Hatchery[[", s, "]], \"p_mature_HOS\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage, proyears = SOM@proyears, ylab = \"Proportion mature\")"),
    "}",
    "```",
    "",
    "```{r, fig.cap=\"Broodtake fecundity.\"}",
    paste0("plot_SOM(SOM@Hatchery[[", s, "]], \"fec_brood\", type = \"age\", maxage = SOM@Bio[[", s, "]]@maxage, nsim = SOM@nsim, proyears = SOM@proyears, ylab = \"Brood fecundity\")"),
    "```",
    "\n",
    "#### Time series",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin broodtake (NOB) time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"NOB\", s = ", s, ", quant = TRUE, ylab = \"NOB\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin broodtake (NOB), histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"NOB\", s = ", s, ", xlab = \"NOB\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Local hatchery origin broodtake (HOB) time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"HOB\", s = ", s, ", quant = TRUE, ylab = \"Local HOB\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Local hatchery origin broodtake (HOB), histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"HOB\", s = ", s, ", xlab = \"Local HOB\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Imported hatchery origin broodtake (HOB) time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"HOB_import\", s = ", s, ", quant = TRUE, ylab = \"Imported HOB\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Imported hatchery origin broodtake (HOB), histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"HOB_import\", s = ", s, ", xlab = \"Imported HOB\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion natural origin broodtake NOB/(NOB + HOB) time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"pNOB\", s = ", s, ", quant = TRUE, ylab = \"P(NOB)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion natural origin broodtake, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"pNOB\", s = ", s, ", xlab = \"P(NOB)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion broodtake (ratio of total broodtake to total escapement) time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"pbrood\", s = ", s, ", quant = TRUE, ylab = \"P(brood)\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion broodtake, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"pbrood\", s = ", s, ", xlab = \"P(brood)\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery smolt releases time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"Smolt_Rel\", s = ", s, ", quant = TRUE, ylab = \"Smolt releases\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery smolt releases, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"Smolt_Rel\", s = ", s, ", xlab = \"Smolt releases\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Fitness (of spawners to natural environment) time series (annual median and 95 percent interval).\"}",
    "SMSEf <- SMSE",
    "SMSEf@fitness <- array(SMSE@fitness[, , 1, ], c(SMSE@nsim, SMSE@nstocks, SMSE@proyears))",
    paste0("plot_statevar_ts(SMSEf, var = \"fitness\", s = ", s, ", quant = TRUE, ylab = \"Fitness\")"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Fitness (of spawners to natural environment), histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSEf, var = \"fitness\", s = ", s, ", xlab = \"Fitness\")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"PNI time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"PNI\", s = ", s, ", quant = TRUE)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"PNI, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"PNI\", s = ", s, ")"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    ""
  )
}

rmd_fishery <- function(s, sname) {
  c(
    "### Fishery {.tabset}",
    "\n",
    "#### Parameters",
    "\n",
    "```{r results=\"asis\"}",
    paste0("knitr::kable(make_harvest_table(SMSE, s = ", s, "), row.names = FALSE, caption = \"Harvest parameters.\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Vulnerability at age for the preterminal fishery (median and 95 percent interval across simulations).\"}",
    paste0("plot_SOM(SOM@Harvest[[", s, "]], \"vulPT\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage, ylab = \"Preterminal vulnerability\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Vulnerability at age for the terminal fishery (median and 95 percent interval across simulations).\"}",
    paste0("plot_SOM(SOM@Harvest[[", s, "]], \"vulT\", type = \"age\", nsim = SOM@nsim, maxage = SOM@Bio[[", s, "]]@maxage, ylab = \"Terminal vulnerability\")"),
    "```",
    "\n",
    "#### Time series",
    "\n",
    "```{r, fig.cap=\"Kept catch in the preterminal (PT) and terminal (T) fishery (annual medians across simulations).\"}",
    paste0("plot_fishery(SMSE, s = ", s, ", type = \"catch\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Harvest rate (from kept catch and abundance) in the preterminal (PT) and terminal (T) fishery (annual medians across simulations).\"}",
    paste0("plot_fishery(SMSE, s = ", s, ", type = \"harvest\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Exploitation rate (from kept catch and dead discards and abundance) in the preterminal (PT) and terminal (T) fishery (annual medians across simulations).\"}",
    paste0("plot_fishery(SMSE, s = ", s, ", type = \"ex\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"NOS/SMSY time series (annual median and 95 percent interval).\"}",
    paste0("plot_statevar_ts(SMSE, var = \"NOS/SMSY\", s = ", s, ", quant = TRUE, ylab = expression(NOS/S[MSY]))"),
    "```",
    "\n",
    "```{r, fig.cap=\"NOS/SMSY, histogram in last projection year.\"}",
    paste0("plot_statevar_hist(SMSE, var = \"NOS/SMSY\", s = ", s, ", xlab = expression(NOS/S[MSY]))"),
    "```",
    "\n",
    "```{r, fig.cap=\"Kobe figure (annual medians) with exploitation rate for the preterminal fishery.\"}",
    paste0("plot_Kobe(SMSE, s = ", s, ", type = \"PT\")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Kobe figure (annual medians) with exploitation rate for the terminal fishery.\"}",
    paste0("plot_Kobe(SMSE, s = ", s, ", type = \"T\")"),
    "```",
    "\n",
    "### Comparison",
    "\n",
    "```{r, fig.cap=\"Composition of spawners (annual medians across simulations).\"}",
    paste0("plot_spawners(SMSE, s = ", s, ", prop = FALSE)"),
    "```",
    "\n",
    "```{r, fig.cap=\"Composition of spawners (proportions from annual medians across simulations).\"}",
    paste0("plot_spawners(SMSE, s = ", s, ", prop = TRUE)"),
    "```",
    "\n",
    "```{r, fig.cap=\"Ratio of spawners and broodstock to escapement (proportions from annual medians across simulations).\"}",
    paste0("plot_escapement(SMSE, s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Fitness and other metrics related to hatchery production (annual medians across simulations).\"}",
    paste0("plot_fitness(SMSE, s = ", s, ")"),
    "```",
    "\n"
  )
}

rmd_lhg <- function(s, sname) {
  c(
    "#### Life history groups",
    "\n",
    "```{r, fig.cap=\"Smolt production by life history group (median across simulations).\"}",
    paste0("plot_LHG(SMSE, var = \"Smolt\", type = \"abs\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Natural origin spawners by life history group (median across simulations).\"}",
    paste0("plot_LHG(SMSE, var = \"NOS\", type = \"abs\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of life history groups at the smolt life stage.\"}",
    paste0("plot_LHG(SMSE, var = \"Smolt\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of life history groups of natural spawners.\"}",
    paste0("plot_LHG(SMSE, var = \"NOS\", s = ", s, ")"),
    "```",
    "\n"
  )
}

rmd_rs <- function(s, sname) {
  c(
    "#### Release strategies",
    "\n",
    "```{r, fig.cap=\"Hatchery releases at smolt life stage by release strategy (median across simulations).\"}",
    paste0("plot_RS(SMSE, var = \"Smolt\", ylab = \"Hatchery smolts\", type = \"abs\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Hatchery origin spawners by release strategy  (median across simulations).\"}",
    paste0("plot_RS(SMSE, var = \"HOS\", type = \"abs\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of release strategy at the smolt life stage.\"}",
    paste0("plot_RS(SMSE, var = \"Smolt\", ylab = \"Proportion hatchery smolts\", s = ", s, ")"),
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of release strategy of hatchery spawners.\"}",
    paste0("plot_RS(SMSE, var = \"HOS\", s = ", s, ")"),
    "```",
    "\n"
  )
}

make_rmd_stock <- function(s, sname, n_g = 1, n_r = 1) {

  natural_production <- rmd_natural_production(s, sname)
  if (n_g > 1) {
    lhg <- rmd_lhg(s, sname)
  } else {
    lhg <- NULL
  }

  hatchery_production <- rmd_hatchery_production(s, sname)
  if (n_r > 1) {
    rs <- rmd_rs(s, sname)
  } else {
    rs <- NULL
  }

  fishery <- rmd_fishery(s, sname)

  c(natural_production, lhg, hatchery_production, rs, fishery)
}

make_rmd_stock_comparison <- function() {

  c(
    "## Among-population comparison",
    "\n",
    "```{r, fig.cap=\"Stray matrix of hatchery origin fish at the recruitment stage.\"}",
    "plot_stray(SMSE@Misc$SOM@stray)",
    "```",
    "\n",
    "```{r, fig.cap=\"Natural origin spawners time series (annual medians) across population.\"}",
    "plot_statevar_ts(SMSE, var = \"NOS\", s = 1:SMSE@nstocks)",
    "```",
    "\n",
    "```{r, fig.cap=\"Hatchery origin spawners time series (annual medians) across population.\"}",
    "plot_statevar_ts(SMSE, var = \"HOS\", s = 1:SMSE@nstocks)",
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of hatchery spawners (annual medians) across populations..\"}",
    "plot_statevar_ts(SMSE, var = \"pHOScensus\", s = 1:SMSE@nstocks, ylab = \"P(HOS)\")",
    "```",
    "\n",
    "```{r, fig.cap=\"Proportion of wild spawners (annual medians) across populations.\"}",
    "plot_statevar_ts(SMSE, var = \"p_wild\", s = 1:SMSE@nstocks, ylab = \"P(WILD)\")",
    "```",
    "\n",
    "```{r, fig.cap=\"Natural smolt production (annual medians) across populations.\"}",
    "plot_statevar_ts(SMSE, var = \"Smolt\", s = 1:SMSE@nstocks, ylab = \"Natural smolt production\")",
    "```",
    "\n",
    "```{r, fig.cap=\"Hatchery smolt releases (annual medians) across populations.\"}",
    "plot_statevar_ts(SMSE, var = \"Smolt_Rel\", s = 1:SMSE@nstocks, ylab = \"Hatchery smolt releases\")",
    "```",
    "\n",
    "```{r, fig.cap=\"Fitness (of spawners to natural environment) (annual medians) across populations.\"}",
    "plot_statevar_ts(SMSEf, var = \"fitness\", s = 1:SMSE@nstocks, ylab = \"Fitness\")",
    "```",
    "\n",
    "```{r, fig.cap=\"PNI (proportionate natural influence) (annual medians) across populations.\"}",
    "plot_statevar_ts(SMSE, var = \"PNI\", s = 1:SMSE@nstocks)",
    "```",
    "\n"
  )
}


make_rmd_compare <- function(s, sname) {

  nat <- c(paste("##", sname, "{.tabset}"),
    "### Natural production",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin spawners time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"NOS\", s = ", s, ", quant = TRUE, names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin spawners, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"NOS\", s = ", s, ", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery origin spawners time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"HOS\", s = ", s, ", quant = TRUE, names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery origin spawners, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"HOS\", s = ", s, ", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of hatchery spawners (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"pHOS_census\", s = ", s, ", quant = TRUE, ylab = \"P(HOS)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of hatchery spawners, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"pHOS_census\", s = ", s, ", xlab = \"P(HOS)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of wild spawners (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"p_wild\", s = ", s, ", quant = TRUE, ylab = \"P(WILD)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion of wild spawners, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"p_wild\", s = ", s, ", xlab = \"P(WILD)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, natural origin (annual median and 95 percent interval). Losses include en-route mortality and broodtake.\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"pNOSesc\", s = ", s, ", quant = TRUE, ylab = \"Spawner/escapement (natural)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, natural origin, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"pNOSesc\", s = ", s, ", xlab = \"Spawner/escapement (natural)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, hatchery origin (annual median and 95 percent interval). Losses include en-route mortality, broodtake, and removal from spawning grounds.\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"pHOSesc\", s = ", s, ", quant = TRUE, ylab = \"Spawner/escapement (hatchery)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Spawner to escapement ratio, hatchery origin, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"pHOSesc\", s = ", s, ", xlab = \"Spawner/escapement (hatchery)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Smolt production (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"Smolt\", s = ", s, ", quant = TRUE, ylab = \"Smolt production\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Smolt production, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"Smolt\", s = ", s, ", xlab = \"Smolt production\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    ""
  )

  hat <- c(
    "### Hatchery production",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin broodtake (NOB) time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"NOB\", s = ", s, ", quant = TRUE, ylab = \"NOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Natural origin broodtake (NOB), histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"NOB\", s = ", s, ", xlab = \"NOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Local hatchery origin broodtake (HOB) time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"HOB\", s = ", s, ", quant = TRUE, ylab = \"Local HOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Local hatchery origin broodtake (HOB), histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"HOB\", s = ", s, ", xlab = \"Local HOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Imported hatchery origin broodtake (HOB) time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"HOB_import\", s = ", s, ", quant = TRUE, ylab = \"Imported HOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Imported hatchery origin broodtake (HOB), histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"HOB_import\", s = ", s, ", xlab = \"Imported HOB\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion natural origin broodtake NOB/(NOB + HOB) time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"pNOB\", s = ", s, ", quant = TRUE, ylab = \"P(NOB)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion natural origin broodtake, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"pNOB\", s = ", s, ", xlab = \"P(NOB)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion broodtake (ratio of total broodtake to total escapement) time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"pbrood\", s = ", s, ", quant = TRUE, ylab = \"P(brood)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Proportion broodtake, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"pbrood\", s = ", s, ", xlab = \"P(brood)\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery smolt releases time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"Smolt_Rel\", s = ", s, ", quant = TRUE, ylab = \"Smolt releases\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Hatchery smolt releases, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"Smolt_Rel\", s = ", s, ", xlab = \"Smolt releases\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Fitness (of spawners to natural environment) time series (annual median and 95 percent interval).\"}",
    "SMSEf_list <- lapply(SMSE_list, function(SMSE) {",
    "  SMSE@fitness <- array(SMSE@fitness[, , 1, ], c(SMSE@nsim, SMSE@nstocks, SMSE@proyears))",
    "  return(SMSE)",
    "})",
    paste0("compare_statevar_ts(SMSEf_list, var = \"fitness\", s = ", s, ", quant = TRUE, ylab = \"Fitness\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"Fitness (of spawners to natural environment), histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSEf_list, var = \"fitness\", s = ", s, ", xlab = \"Fitness\", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    "",
    "<div style=\"display: flex;\">",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"PNI time series (annual median and 95 percent interval).\"}",
    paste0("compare_statevar_ts(SMSE_list, var = \"PNI\", s = ", s, ", quant = TRUE, names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "<div style=\"width: 50%;\">",
    "",
    "```{r, fig.cap=\"PNI, histogram in last projection year.\"}",
    paste0("compare_statevar_hist(SMSE_list, var = \"PNI\", s = ", s, ", names = names, col_vec = col_vec)"),
    "```",
    "",
    "</div>",
    "",
    "</div>",
    ""
  )

  c(nat, hat)

}