---
title: "4. Projections and reference points"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{4. Projections and reference points}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE,
  fig.width = 7,
  fig.height = 5
)
```

Once a model has been estimated, Rceattle can project the population forward under alternative fishing mortality scenarios and recruitment assumptions. This vignette covers:

- Fixed-F projections (including mean historical F)
- Manual and stochastic recruitment in the projection period
- Random recruitment effects
- Visualising projection uncertainty across scenarios

## Setup

```{r}
library(Rceattle)

data(BS2017SS)
data(BS2017MS)

# Set the terminal projection year
BS2017SS$projyr <- 2060
BS2017MS$projyr <- 2060

# Plot input data
plot_data(BS2017SS)

# Fit baseline single- and multi-species models
ss_run <- fit_mod(
  data_list    = BS2017SS,
  inits        = NULL,
  file         = NULL,
  estimateMode = 0,
  random_rec   = FALSE,
  msmMode      = 0,
  fit_control = fit_control(
    phase        = TRUE,
    verbose      = 1))

ms_run <- fit_mod(
  data_list    = BS2017MS,
  inits        = ss_run$estimated_params,
  file         = NULL,
  estimateMode = 0,
  niter        = 5,
  random_rec   = FALSE,
  msmMode      = 1,
  suitMode     = 0,
  fit_control = fit_control(
    verbose      = 1))
```

## Fixed-F projections

`estimateMode = 2` runs only the projection period (holding all estimated parameters fixed) under a specified harvest control rule. `build_hcr(HCR = 2, Ftarget = ...)` sets a constant F for each species.

The `proj_F_prop` column in `fleet_control` controls how projected F is split across fleets when there are multiple fisheries targeting the same species (values must sum to 1 per species).

```{r}
BS2017MS$fleet_control$proj_F_prop <- 1  # One fishery fleet per species

ms_run_proj <- fit_mod(
  data_list    = BS2017MS,
  inits        = ms_run$estimated_params,
  file         = NULL,
  estimateMode = 2,
  HCR          = build_hcr(HCR = "ConstantF",
                            Ftarget = c(0.2342936, 0.513, 0.0774777)),
  niter        = 5,
  random_rec   = FALSE,
  msmMode      = 1,
  suitMode     = 0,
  fit_control = fit_control(
    verbose      = 1))

plot_catch(ms_run_proj, incl_proj = TRUE)
```

## Recruitment in the projection period

### Manual recruitment deviations

Projection recruitment deviations are stored in the parameter list and can be replaced directly before a forward-projection run. `estimateMode = 3` evaluates the model without re-estimating any parameters, so only the updated recruitment deviations affect the projection.

```{r}
nyrs      <- BS2017MS$endyr - BS2017MS$styr + 1
nyrs_proj <- BS2017MS$projyr - BS2017MS$styr + 1
yrs_proj  <- (nyrs + 1):nyrs_proj

# Replace future rec_devs with draws from N(0, 0.707)
ms_run$estimated_params$rec_dev[, yrs_proj] <- stats::rnorm(
  n    = length(ms_run$estimated_params$rec_dev[, yrs_proj]),
  mean = 0,
  sd   = 0.707
)

ms_run_proj2 <- fit_mod(
  data_list    = BS2017MS,
  inits        = ms_run$estimated_params,
  file         = NULL,
  estimateMode = 3,   # Evaluate only — do not re-estimate
  niter        = 5,
  random_rec   = FALSE,
  msmMode      = 1,
  suitMode     = 0,
  fit_control = fit_control(
    verbose      = 1))
```

### Stochastic recruitment via `sample_rec()`

`sample_rec()` bootstraps historical recruitment deviations into the projection period. The optional `rec_trend` argument adds a linear trend: `rec_trend = -0.5` imposes a 50% decline in mean recruitment by the end of the projection.

```{r}
ms_run_proj3 <- sample_rec(
  Rceattle     = ms_run,
  sample_rec   = TRUE,
  update_model = TRUE,
  rec_trend    = -0.5   # 50% decline in mean recruitment
)
```

## Comparing projection scenarios

```{r}
mod_list  <- list(ms_run, ms_run_proj, ms_run_proj2, ms_run_proj3)
mod_names <- c(
  "No fishing",
  "Mean historical F",
  "Mean F + stochastic recruitment",
  "Mean F + declining recruitment"
)

plot_biomass(Rceattle = mod_list, model_names = mod_names, incl_proj = TRUE)
plot_recruitment(Rceattle = mod_list, model_names = mod_names, incl_proj = TRUE)
plot_catch(Rceattle = mod_list, model_names = mod_names, incl_proj = TRUE)
```

## Random recruitment effects

Setting `random_rec = TRUE` treats log-recruitment deviations as random effects marginalised via the Laplace approximation rather than as penalised fixed effects. This propagates process uncertainty into projections automatically.

```{r}
ss_re <- fit_mod(
  data_list    = BS2017SS,
  inits        = NULL,
  file         = NULL,
  estimateMode = 0,
  random_rec   = TRUE,
  msmMode      = 0,
  fit_control = fit_control(
    phase        = FALSE,
    verbose      = 1))

# Compare fixed vs. random effects recruitment uncertainty
plot_recruitment(
  Rceattle    = list(ss_run, ss_re),
  model_names = c("Fixed effects", "Random effects"),
  add_ci      = TRUE
)
```

## Available harvest control rules

The `HCR` argument to `fit_mod()` accepts the output of `build_hcr()`. See `?build_hcr` for full parameter details.

| Value | String | Description | Multi-species? |
|---|---|---|---|
| 0 | `"NoFishing"` | No fishing — estimate hindcast only | Yes |
| 1 | `"CMSY"` | CMSY: maximise catch across all species (can constrain depletion ≥ `Plimit`) | Yes |
| 2 | `"ConstantF"` | Constant F set at `Ftarget` for each species | Yes |
| 3 | `"ConstantFSSB"` | F that achieves `Ftarget`% of SSB₀ at the end of the projection | Yes |
| 4 | `"ConstantFSPR"` | Constant F_SPR set at `Ftarget`; can be scaled by `Fmult` (NEFSC convention) | No |
| 5 | `"NPFMC"` | NPFMC Tier 3 SPR-based HCR (`Ftarget`, `Flimit`, `Ptarget`, `Alpha`) | No |
| 6 | `"PFMC"` | PFMC Category 1 40-10 ACL rule with uncertainty buffer (`Pstar`, `Sigma`) | Yes |
| 7 | `"SESSF"` | SESSF Tier 1 SPR-based HCR | No |

```{r}
# Example: NPFMC Tier 3 control rule (single-species only)
ss_run_tier3 <- fit_mod(
  data_list    = BS2017SS,
  inits        = ss_run$estimated_params,
  file         = NULL,
  estimateMode = 2,
  HCR          = build_hcr(HCR = "NPFMC",
                            Ftarget = 0.40,   # F40%
                            Flimit  = 0.35,   # F35%
                            Ptarget = 0.40,
                            Plimit  = 0.20,
                            Alpha   = 0.05),
  random_rec   = FALSE,
  msmMode      = 0,
  fit_control = fit_control(
    verbose      = 1))
```

For closed-loop testing of control rules across many simulated trajectories, see `vignette("hcrs-and-mses")`.
