Introduction Data & base model Sensitivities Grid Projections

More

Getting started GT skipping Review B10+ and TRO OMMP16 Model notes ADMB vs SBT Data weighting Length frequency Packages sbt vignettes sbtdata vignettes

ESC31 Projections

Publication and diagnostic status: There is one canonical projection report, 5_proj.html. Its production publication uses all 2,000 balanced posterior draws. The same source can be rendered cache-only with 100 draws for an abbreviated smoke review. The retained format-5 200-draw caches are legacy evidence and require refresh before they can pass the current audit; there is no separate diagnostic HTML. All nine selected grid cells, the compatible balanced posterior, the replacement 108-fit direct-get_M MLE grid, and its 2,000-draw resample pass their current validation gates. This report uses the completed and validated 2,000-draw production caches for both allocation scenarios, NoUAM, and the sampled direct-M MLE-grid comparison.

This page runs projections through the run_projection_scenario() wrapper in the Projection Dynamics section below. That wrapper passes the balanced posterior draws, projected recruitment deviates, projected selectivity, fixed pre-CTP catch inputs, catch-allocation splits, and monitoring-sample settings to sbt::run_projections().

The selected posterior is conditioned by the approved preventive harvest wall used during model fitting. Projection dynamics do not add that fitting wall to future states or TAC paths: they advance the population with the package projection dynamics and retain the separate feasibility-continuation and positive-state checks. The wall contract is nevertheless carried in every projection source identity so a posterior fitted under a different wall cannot be reused silently.

In the sbt package, run_projections() is implemented in R/projections.R. That file also provides project_rec_devs(), project_selectivity(), project_ctp_schedule(), and the internal helpers that build catch arrays, advance projected population dynamics, simulate future CPUE, GT, POP, and HSP observations, and cache timing/results. Projection dynamics use the fitted RTMB assessment model from R/model.R, with shared population dynamics in R/dynamics.R, selectivity in R/selectivity.R, recruitment in R/recruitment.R, and observation likelihood/simulation helpers in R/likelihoods.R. CTP TAC updates are calculated through CTP() in R/CTP.R, which also contains the CKMR CTP objective ctp_ckmrmod().

options(readr.show_col_types = FALSE)

library(knitr)
library(tidyverse)
library(RTMB)
library(SparseNUTS)
library(scales)
library(sbt)

theme_set(theme_bw())

`%||%` <- function(x, y) {
  if (is.null(x) || length(x) == 0L) y else x
}

format_decimal <- function(x, digits = 3) {
  out <- rep(NA_character_, length(x))
  finite <- is.finite(x)
  x[finite & abs(x) < 0.5 * 10^(-digits)] <- 0
  out[finite] <- formatC(
    x[finite],
    format = "f",
    digits = digits,
    big.mark = ","
  )
  out[!finite & !is.na(x)] <- as.character(x[!finite & !is.na(x)])
  out
}

zero_y_breaks <- function(n = 5) {
  force(n)
  function(x) {
    rng <- range(c(0, x), finite = TRUE)
    if (!all(is.finite(rng))) return(0)
    breaks <- pretty(rng, n = n)
    sort(unique(c(0, breaks[breaks >= 0 & breaks <= rng[[2L]]])))
  }
}

scale_y_zero <- function(labels = waiver(), n = 5) {
  scale_y_continuous(
    limits = c(0, NA),
    breaks = zero_y_breaks(n),
    labels = labels,
    expand = expansion(mult = c(0, 0.05))
  )
}

projection_quantiles <- function(x) {
  tibble(
    lower = quantile(x, 0.025, na.rm = TRUE),
    median = median(x, na.rm = TRUE),
    upper = quantile(x, 0.975, na.rm = TRUE)
  )
}

projection_lapply <- function(x, fun, cores = 1L) {
  cores <- min(max(as.integer(cores), 1L), length(x))
  if (.Platform$OS.type == "unix" && cores > 1L) {
    out <- parallel::mclapply(x, fun, mc.cores = cores, mc.preschedule = FALSE)
    if (!any(vapply(out, inherits, logical(1), "try-error"))) {
      return(out)
    }
  }
  lapply(x, fun)
}

read_rds_cache <- function(path) {
  if (!file.exists(path)) return(NULL)
  tryCatch(
    suppressWarnings(readRDS(path)),
    error = function(e) NULL
  )
}

load_rdata_cache <- function(path) {
  if (!file.exists(path)) return(NULL)

  tryCatch({
    cache_env <- new.env(parent = emptyenv())
    loaded <- suppressWarnings(load(path, envir = cache_env))
    if (!length(loaded)) stop("the file contains no objects")
    cache_env
  }, error = function(e) NULL)
}

atomic_cache_write <- function(path, writer, rename_file = file.rename) {
  dir.create(dirname(path), recursive = TRUE, showWarnings = FALSE)
  cache_dir <- normalizePath(dirname(path), mustWork = TRUE)
  path <- file.path(cache_dir, basename(path))
  temporary <- tempfile(
    pattern = paste0(".", basename(path), "-"),
    tmpdir = cache_dir
  )
  on.exit(unlink(temporary), add = TRUE)

  writer(temporary)
  temporary_size <- file.size(temporary)
  if (!file.exists(temporary) || is.na(temporary_size) || temporary_size <= 0) {
    stop("Cache writer did not create a valid temporary file for ", path, call. = FALSE)
  }

  rename_ok <- function(from, to) {
    isTRUE(suppressWarnings(rename_file(from, to)))
  }
  if (rename_ok(temporary, path)) return(invisible(path))
  if (!file.exists(path)) {
    stop("Could not move the completed cache into place: ", path, call. = FALSE)
  }

  backup <- tempfile(
    pattern = paste0(".", basename(path), "-backup-"),
    tmpdir = cache_dir
  )
  if (!rename_ok(path, backup)) {
    stop("Could not stage the existing cache for replacement: ", path, call. = FALSE)
  }
  if (!rename_ok(temporary, path)) {
    restored <- rename_ok(backup, path)
    if (restored) {
      stop("Could not replace the cache; the previous file was restored: ", path,
           call. = FALSE)
    }
    stop(
      "Could not replace or restore the cache; the previous file remains at: ",
      backup,
      call. = FALSE
    )
  }
  if (unlink(backup) != 0L) {
    warning("The previous cache could not be removed from: ", backup, call. = FALSE)
  }
  invisible(path)
}

atomic_save_rds <- function(object, path) {
  atomic_cache_write(path, function(temporary) saveRDS(object, temporary))
}

cache_record_is_current <- function(cache, signature, fields) {
  is.list(cache) &&
    all(c("signature", fields) %in% names(cache)) &&
    identical(cache$signature, signature)
}

cache_record_is_compatible <- function(
    cache, signature, fields, legacy_fingerprint = NULL) {
  if (!is.list(cache) ||
      !all(c("signature", fields) %in% names(cache))) {
    return(FALSE)
  }
  if (identical(cache$signature, signature)) return(TRUE)
  if (is.null(legacy_fingerprint) ||
      !is.list(legacy_fingerprint) ||
      !all(c("signature_md5", "payload_md5") %in%
        names(legacy_fingerprint))) {
    return(FALSE)
  }

  payload <- if (length(fields) == 1L) {
    cache[[fields]]
  } else {
    cache[fields]
  }
  identical(
    esc31_object_md5(cache$signature),
    legacy_fingerprint$signature_md5
  ) && identical(
    esc31_object_md5(payload),
    legacy_fingerprint$payload_md5
  )
}

Configuration

This document reads the available ESC31 grid MCMC outputs and builds a balanced projection posterior with approximately equal random draws from each completed grid cell. A normal render reads the accepted 2,000-draw production caches. ESC31_PROJECTION_N_ITER=100 selects the abbreviated smoke caches and ESC31_PROJECTION_N_ITER=200 selects the retained extended diagnostic caches; all three settings target the single canonical 5_proj.html. The format-5 200-draw caches predate the raw CTP-input snapshots and expanded optimizer record required by the independent audit. They therefore remain legacy diagnostic evidence but cannot produce a current 200-draw report until that arm is refreshed. The 100-draw and 2,000-draw caches carry the complete audit contract. Rendering remains cache-only: it cannot create projected recruitment, selectivity, CPUE-q, historical-overlay, or dynamics caches unless the corresponding ESC31_RUN_PROJECTIONS, ESC31_RUN_NO_UAM_PROJECTIONS, or ESC31_RUN_MLE_GRID_PROJECTIONS flag is set explicitly after its source artifact is accepted.

if (basename(getwd()) == "ESC31") {
  esc_dir <- normalizePath(".")
} else if (dir.exists("ESC31")) {
  esc_dir <- normalizePath("ESC31")
} else {
  stop("Could not find the ESC31 directory.", call. = FALSE)
}

projection_first_yr <- 2022
projection_last_yr <- 2035
projection_years <- projection_first_yr:projection_last_yr
projection_n_iter <- as.integer(Sys.getenv("ESC31_PROJECTION_N_ITER", "2000"))
projection_mcmc_source_draws <- as.integer(Sys.getenv("ESC31_PROJECTION_MCMC_SOURCE_DRAWS", "2000"))
projection_seed <- 102
projection_draw_seed <- 42L
no_uam_projection_draw_seed <- 43L
mle_grid_projection_draw_seed <- 45L
if (!projection_n_iter %in% c(100L, 200L, 2000L)) {
  stop(
    "ESC31_PROJECTION_N_ITER must be exactly 100 or 200 for diagnostic ",
    "review, or 2000 for the accepted production run.",
    call. = FALSE
  )
}
projection_diagnostic_run <- projection_n_iter < 2000L
if (!identical(projection_mcmc_source_draws, 2000L)) {
  stop("ESC31_PROJECTION_MCMC_SOURCE_DRAWS must be exactly 2000.", call. = FALSE)
}

projection_env_flag <- function(name, default = FALSE) {
  value <- tolower(trimws(Sys.getenv(name, if (default) "true" else "false")))
  if (!value %in% c("true", "false", "1", "0", "yes", "no")) {
    stop(name, " must be true or false.", call. = FALSE)
  }
  value %in% c("true", "1", "yes")
}
run_projection_computation <- projection_env_flag("ESC31_RUN_PROJECTIONS")
run_no_uam_projection_computation <- projection_env_flag(
  "ESC31_RUN_NO_UAM_PROJECTIONS",
  default = run_projection_computation
)
run_mle_grid_projection_computation <- projection_env_flag(
  "ESC31_RUN_MLE_GRID_PROJECTIONS"
)
require_projection_computation <- function(stage) {
  if (!run_projection_computation) {
    stop(
      "The compatible ", stage, " cache is unavailable. Set ",
      "ESC31_RUN_PROJECTIONS=true only after accepting the approved grid.",
      call. = FALSE
    )
  }
}
require_no_uam_projection_computation <- function(stage) {
  if (!run_no_uam_projection_computation) {
    stop(
      "The compatible NoUAM ", stage, " cache is unavailable. Set ",
      "ESC31_RUN_NO_UAM_PROJECTIONS=true to build the comparison arm.",
      call. = FALSE
    )
  }
}
require_mle_grid_projection_computation <- function(stage) {
  if (!run_mle_grid_projection_computation) {
    stop(
      "The compatible sampled-MLE-grid ", stage,
      " cache is unavailable. Set ",
      "ESC31_RUN_MLE_GRID_PROJECTIONS=true to build this report-specific ",
      "comparison arm.",
      call. = FALSE
    )
  }
}

run_dir <- file.path(esc_dir, "runs")
dir.create(run_dir, recursive = TRUE, showWarnings = FALSE)
projection_run_dir <- file.path(run_dir, "projections")
dir.create(projection_run_dir, recursive = TRUE, showWarnings = FALSE)

model_file <- file.path(run_dir, "esc31_base.rds")
no_uam_model_file <- file.path(
  run_dir, "sens", "esc31_sens_no_uam.sbt.rds"
)
combined_mcmc_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_mcmc_", projection_mcmc_source_draws, ".rda")
)
projection_run_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_run_", projection_n_iter, ".rds")
)
projection_run_scenario_2_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_run_scenario_2_", projection_n_iter, ".rds")
)
projection_recdev_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_recdev_", projection_n_iter, ".rds")
)
projection_selectivity_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_selectivity_", projection_n_iter, ".rds")
)
projection_cpue_q_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_cpue_q_", projection_n_iter, ".rds")
)
projection_history_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_history_", projection_n_iter, ".rds")
)
no_uam_projection_run_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_run_no_uam_", projection_n_iter, ".rds")
)
no_uam_projection_recdev_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_recdev_no_uam_", projection_n_iter, ".rds")
)
no_uam_projection_selectivity_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_selectivity_no_uam_", projection_n_iter, ".rds")
)
mle_grid_projection_run_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_run_mle_grid_", projection_n_iter, ".rds")
)
mle_grid_projection_recdev_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_recdev_mle_grid_", projection_n_iter, ".rds")
)
mle_grid_projection_selectivity_file <- file.path(
  projection_run_dir,
  paste0(
    "esc31_projection_selectivity_mle_grid_",
    projection_n_iter,
    ".rds"
  )
)

# These fingerprints bind the pre-split upstream caches to their exact stored
# signatures and payloads. They permit a one-time compatibility read after the
# projection workflow identity was narrowed, without allowing an arbitrary
# stale or modified cache through. Future caches use component-specific source
# identities and do not need this bridge.
legacy_projection_input_fingerprints <- list(
  `100` = list(
    base_recdev = list(
      signature_md5 = "c32f00ca5181898cf79702e7ce1f919f",
      payload_md5 = "ddadfcd653c004afa65d4b9236236ce5"
    ),
    base_selectivity = list(
      signature_md5 = "a1139d297b20f22fc422217fa1e11a1c",
      payload_md5 = "b818419e8ea538d12cbb5b1dbe2c9ce3"
    ),
    base_cpue_q = list(
      signature_md5 = "41752b47d075313bfb92c2b3fed7f033",
      payload_md5 = "359470ee0fe4726ce8a060555017be84"
    ),
    no_uam_recdev = list(
      signature_md5 = "56fd00ce36b12da8d272a24087968ceb",
      payload_md5 = "05d462e5f2beeb471bdcd7776931742c"
    ),
    no_uam_selectivity = list(
      signature_md5 = "4627e4dee2b6562d8076888e8eafda33",
      payload_md5 = "2965b01970321b464e7ec8217dac267d"
    ),
    mle_grid_recdev = list(
      signature_md5 = "027ccd0ffd2f0fcbb94d1fa23e8f7df5",
      payload_md5 = "394aa57fa8d92e7bf814e0543cc871d6"
    ),
    mle_grid_selectivity = list(
      signature_md5 = "6d6dc19c667b0bc076e30a9bec0bb622",
      payload_md5 = "f716e2309ec0400c9ffe2f3171dd1c39"
    ),
    history = list(
      signature_md5 = "0bc9c4076130ee7d4522fbc3ee7b4ba0",
      payload_md5 = "4c359de4ec20f3fe1c6d1d8e6d9ff323"
    )
  ),
  `200` = list(
    base_recdev = list(
      signature_md5 = "48cfeacb093b0fb87f02a66af088768f",
      payload_md5 = "dfea54739d522de5808fe7226b6f33a1"
    ),
    base_selectivity = list(
      signature_md5 = "4a81f4e897cacaffa91ca020791433b2",
      payload_md5 = "dcc94db2decaf1675cbb27112786b5e6"
    ),
    base_cpue_q = list(
      signature_md5 = "aa656bc9d73b09da7d8b6b1de3d30ae2",
      payload_md5 = "4da36d7fe3c507feae6d91fcb981c9c0"
    ),
    no_uam_recdev = list(
      signature_md5 = "9dfbd378205765810b1a1da21c200ea9",
      payload_md5 = "234ea0eb6057709fec95089d11024c2c"
    ),
    no_uam_selectivity = list(
      signature_md5 = "fec3cddc29fb2015d60caed51c980b5d",
      payload_md5 = "c1a9da24a1c378119fb8803a78c54fdc"
    ),
    mle_grid_recdev = list(
      signature_md5 = "5fb581989448fed5b6132ec71cedf17b",
      payload_md5 = "f63269414cbeebfc348993dcd04d8141"
    ),
    mle_grid_selectivity = list(
      signature_md5 = "90a2acf342f348fb07bbcf117f11989d",
      payload_md5 = "a72c242a963ef61d22f1e26c999808d0"
    ),
    history = list(
      signature_md5 = "bccf19718207e3f65a4ee3a7132e93d2",
      payload_md5 = "0a76a865f16d4ef4962dbb55ee698cf0"
    )
  ),
  `2000` = list(
    base_recdev = list(
      signature_md5 = "0d380e1765aa022bb2a938755018871e",
      payload_md5 = "c6c0037d826a5f63f0300e7f399e78d4"
    ),
    base_selectivity = list(
      signature_md5 = "c788fd1409ddbe51191a73c90a0f40ad",
      payload_md5 = "b37e176b134616e29fd3c224c1074bfc"
    ),
    base_cpue_q = list(
      signature_md5 = "872f36ee79cf5993cac016bc59bd0034",
      payload_md5 = "c25e96eb170f1a6f06a36ae5362fba45"
    ),
    no_uam_recdev = list(
      signature_md5 = "ab1c689de6b46b5eec7275dd74fec608",
      payload_md5 = "71afbf5d3e9b94f38da8e20ecc3aca02"
    ),
    no_uam_selectivity = list(
      signature_md5 = "2b3509faadad24ffc80d071c245c168e",
      payload_md5 = "53f97faffc705e6a6ab02bc8ac3d1ae9"
    ),
    mle_grid_recdev = list(
      signature_md5 = "3f6cec2be331c0c7db87cf9a225e619c",
      payload_md5 = "7f705dbc3251b1a2ceacddedafb226d4"
    ),
    mle_grid_selectivity = list(
      signature_md5 = "af7a5945f21f790689813a96ecb98bd6",
      payload_md5 = "c23d5f35e04d4946c013569308fc430e"
    ),
    history = list(
      signature_md5 = "7faf568f33e0554593815c3e56443531",
      payload_md5 = "6e1d67a56144b574d1977e87175c72e6"
    )
  )
)
legacy_projection_input_fingerprint <- function(cache_name) {
  fingerprints <- legacy_projection_input_fingerprints[[
    as.character(projection_n_iter)
  ]]
  if (is.null(fingerprints) || !cache_name %in% names(fingerprints)) {
    return(NULL)
  }
  fingerprints[[cache_name]]
}

# Format 8 already contains the complete CTP-input and optimiser audit. The
# format-9 change only stages not-yet-decided future catch as zero and records
# realized catch separately. These exact completed caches used the fail-closed
# rule, incurred no catch shortfall, and remained well inside the harvest wall.
# They are therefore accepted as immutable reporting artifacts, but are not
# relabelled or rewritten as format 9.
legacy_projection_reporting_fingerprints <- list(
  `100` = list(
    base = list(
      file_md5 = "4049c6dffc8033a63e17af8a6afdf8d4",
      run_signature = "75db66f213ebe90385c4e6bc28b9940b",
      maximum_raw_harvest = 0.406548081876978,
      maximum_harvest_penalty = 0
    ),
    scenario_2 = list(
      file_md5 = "fc4b578bf341ef88cb43bfd3614c208a",
      run_signature = "02fe2aaf66c053c54157737891fda301",
      maximum_raw_harvest = 0.334868051865903,
      maximum_harvest_penalty = 0
    ),
    no_uam = list(
      file_md5 = "9ef3148db7ee9753fa8182957d979174",
      run_signature = "591d221b0f92f4b33377b51cab01ed5d",
      maximum_raw_harvest = 0.312555452138961,
      maximum_harvest_penalty = 0
    ),
    mle_grid = list(
      file_md5 = "6ddee07b92df8c069093a38c44e0b5aa",
      run_signature = "e57be63bb480ba320cefb86c4d81a922",
      maximum_raw_harvest = 0.344945143496497,
      maximum_harvest_penalty = 0
    )
  ),
  `2000` = list(
    base = list(
      file_md5 = "2dac1a0bcce8e2d113aa629ba20ae613",
      run_signature = "5ba112cdb8348355cbb2dc4468732ede",
      maximum_raw_harvest = 0.563855117458,
      maximum_harvest_penalty = 0
    ),
    scenario_2 = list(
      file_md5 = "3961386869a7674b02138832ea861ac9",
      run_signature = "81881842ba20e89778b631c1915f7bef",
      maximum_raw_harvest = 0.563855117458,
      maximum_harvest_penalty = 0
    ),
    no_uam = list(
      file_md5 = "49c3da0c16437cb7749ac43ff5db6a4e",
      run_signature = "78890dbefe88aa7044717ea85e6e1320",
      maximum_raw_harvest = 0.560658532424,
      maximum_harvest_penalty = 0
    ),
    mle_grid = list(
      file_md5 = "d343014258d314856f2135810232a1a7",
      run_signature = "03ee9a12ed9aef32e0e3aa42151b798d",
      maximum_raw_harvest = 0.353244342607,
      maximum_harvest_penalty = 0
    )
  )
)
legacy_projection_reporting_fingerprint <- function(arm) {
  fingerprints <- legacy_projection_reporting_fingerprints[[
    as.character(projection_n_iter)
  ]]
  if (is.null(fingerprints) || !arm %in% names(fingerprints)) {
    return(NULL)
  }
  fingerprints[[arm]]
}

# These are compatibility signatures for 100- and 200-draw caches created
# under the immediately preceding signature contract. The projection package
# still checks the scientific inputs and all fields required by the applicable
# cache schema before migrating them; package and dependency version strings
# alone no longer make a cache stale. The format-5 200-draw outputs lack the
# later raw CTP-input snapshots and expanded optimizer audit, so the current
# report retains them only as legacy evidence. The direct-M arm was completed
# during the same transition and therefore has one 200-draw legacy signature.
legacy_projection_run_signatures <- list(
  `100` = c(
    base = "38b6b3ed9a9854e06e9eedfc048ac268",
    scenario_2 = "ac4260d007e59eb8e29704f3f2b43f54",
    no_uam = "b8bc2a975ce145c3b4a88a50e2a27624"
  ),
  `200` = c(
    base = "631717cd9b0177190a78e8bd35d9d69a",
    scenario_2 = "d8d26499abb71b1c1fdfca4e2d6348cd",
    no_uam = "f1cd526e9e78ced5752aab98e23a63a9",
    mle_grid = "1f7a1943dbacb316920b9358c30c082b"
  )
)
legacy_projection_cache_signature <- function(arm) {
  signatures <-
    legacy_projection_run_signatures[[as.character(projection_n_iter)]]
  if (is.null(signatures) || !arm %in% names(signatures)) {
    return(character())
  }
  unname(signatures[[arm]])
}

if (!exists("refit_model", inherits = FALSE)) refit_model <- FALSE
rebuild_projection_dynamics <- FALSE
rebuild_projection_inputs <- rebuild_projection_dynamics
rebuild_projection_history <- FALSE
projection_detected_cores <- parallel::detectCores(logical = FALSE)
if (!is.finite(projection_detected_cores) || projection_detected_cores < 1L) {
  projection_detected_cores <- 1L
}
projection_default_cores <- min(4L, projection_detected_cores)
projection_cores <- as.integer(Sys.getenv("ESC31_PROJECTION_CORES", projection_default_cores))
projection_cores <- min(max(projection_cores, 1L), projection_detected_cores)
projection_state_diagnostic_cores <- suppressWarnings(as.integer(Sys.getenv(
  "ESC31_STATE_DIAGNOSTIC_CORES",
  as.character(min(16L, projection_detected_cores))
)))
if (length(projection_state_diagnostic_cores) != 1L ||
    is.na(projection_state_diagnostic_cores) ||
    projection_state_diagnostic_cores < 1L) {
  stop("ESC31_STATE_DIAGNOSTIC_CORES must be one positive integer.",
       call. = FALSE)
}

active_fixed_projection_tac_scenario <- "scenario_1"
projection_arm_labels <- c(
  base = "Balanced MCMC grid — nominal fleet allocations",
  indonesia_allocation = paste(
    "Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia"
  ),
  no_uam = "NoUAM — nominal fleet allocations (NCNM removed)",
  mle_grid = "Sampled direct-M MLE grid — nominal fleet allocations"
)
projection_scenario_labels <- c(
  scenario_1 = unname(projection_arm_labels[["base"]]),
  scenario_2 = unname(projection_arm_labels[["indonesia_allocation"]])
)
projection_scenario_colors <- setNames(
  c("#D55E00", "#0072B2"),
  unname(projection_scenario_labels)
)
fixed_projection_tac_scenarios <- tribble(
  ~scenario, ~fixed_year_index, ~year, ~total, ~LL1, ~LL2, ~Indonesia, ~Australia,
  "scenario_1", 1L, 2026L, 22670.91, 14325.59, 1534.35, 1248.48, 5562.49,
  "scenario_1", 2L, 2027L, 23647.00, 15062.80, 1652.52, 1370.44, 5561.24,
  "scenario_1", 3L, 2028L, 23647.00, 15062.80, 1652.52, 1370.44, 5561.24,
  "scenario_1", 4L, 2029L, 23647.00, 15062.80, 1652.52, 1370.44, 5561.24,
  "scenario_2", 1L, 2026L, 22670.91, 14325.59, 1534.35, 1248.48, 5562.49,
  "scenario_2", 2L, 2027L, 23647.00, 13151.36, 1442.82, 4196.53, 4856.29,
  "scenario_2", 3L, 2028L, 23647.00, 13151.36, 1442.82, 4196.53, 4856.29,
  "scenario_2", 4L, 2029L, 23647.00, 13151.36, 1442.82, 4196.53, 4856.29
)
fixed_tac_total_check <- fixed_projection_tac_scenarios |>
  select(.data$scenario, .data$year, .data$total) |>
  pivot_wider(names_from = .data$scenario, values_from = .data$total)
fixed_tac_allocation_sum <- rowSums(as.matrix(
  fixed_projection_tac_scenarios[
    c("LL1", "LL2", "Indonesia", "Australia")
  ]
))
if (any(abs(
  fixed_tac_allocation_sum - fixed_projection_tac_scenarios$total
) > 1e-8)) {
  stop(
    "Each fixed-year fleet allocation must sum exactly to its total TAC.",
    call. = FALSE
  )
}
if (any(abs(
  fixed_tac_total_check$scenario_1 -
    fixed_tac_total_check$scenario_2
) > 1e-8)) {
  stop(
    "The two allocation scenarios must have identical fixed total TAC in every year.",
    call. = FALSE
  )
}

# The legacy projection control converts nominal TAC allocations to removals.
# NoUAM removes the LL1 NCNM component but retains the separate Australian
# surface-fishery correction used in both conditioning specifications.
base_projection_removal_multiplier_f <- c(
  LL1 = 1.11,
  LL2 = 1,
  LL3 = 1,
  LL4 = 1,
  Indonesia = 1,
  Australia = 1.20
)
no_uam_projection_removal_multiplier_f <- c(
  LL1 = 1,
  LL2 = 1,
  LL3 = 1,
  LL4 = 1,
  Indonesia = 1,
  Australia = 1.20
)
projection_removal_contract <- list(
  interpretation = paste(
    "Nominal future TAC allocation multiplied by fishery-specific factors;",
    "historical assessment catches are retained without re-multiplication"
  ),
  base = base_projection_removal_multiplier_f,
  no_uam = no_uam_projection_removal_multiplier_f,
  no_uam_change = "Remove the LL1 NCNM addition in conditioning and projections",
  retained_correction = "Australia surface-fishery multiplier 1.20"
)

make_fixed_projection_tac_inputs <- function(scenario_name) {
  fixed_projection_tac_df <- fixed_projection_tac_scenarios |>
    filter(.data$scenario == .env$scenario_name) |>
    arrange(.data$fixed_year_index)
  if (!nrow(fixed_projection_tac_df)) {
    stop("No fixed projection TAC inputs found for ", scenario_name, call. = FALSE)
  }

  fixed_projection_tac <- fixed_projection_tac_df$total
  projection_tac_split_source <- fixed_projection_tac_df |>
    filter(.data$fixed_year_index > 1L)
  projection_tac_split_base_f <- c(
    colSums(as.matrix(projection_tac_split_source[c("LL1", "LL2")])) / sum(projection_tac_split_source$total),
    LL3 = 0,
    LL4 = 0,
    colSums(as.matrix(projection_tac_split_source[c("Indonesia", "Australia")])) / sum(projection_tac_split_source$total)
  )
  projection_tac_split_yf <- matrix(
    rep(projection_tac_split_base_f, each = length(projection_years)),
    nrow = length(projection_years),
    dimnames = list(Year = projection_years, Fishery = names(projection_tac_split_base_f))
  )
  fixed_projection_tac_split_yf <- fixed_projection_tac_df |>
    transmute(
      year = .data$year,
      LL1 = .data$LL1 / .data$total,
      LL2 = .data$LL2 / .data$total,
      LL3 = 0,
      LL4 = 0,
      Indonesia = .data$Indonesia / .data$total,
      Australia = .data$Australia / .data$total
    )
  projection_tac_split_yf[match(fixed_projection_tac_split_yf$year, projection_years), ] <-
    as.matrix(fixed_projection_tac_split_yf[names(projection_tac_split_base_f)])

  list(
    scenario = scenario_name,
    fixed_projection_tac_df = fixed_projection_tac_df,
    fixed_catch_n_years = nrow(fixed_projection_tac_df),
    fixed_projection_tac = fixed_projection_tac,
    projection_tac_split_base_f = projection_tac_split_base_f,
    projection_tac_split_yf = projection_tac_split_yf,
    fixed_projection_tac_split_yf = fixed_projection_tac_split_yf
  )
}

fixed_projection_tac_inputs <- make_fixed_projection_tac_inputs(active_fixed_projection_tac_scenario)
fixed_projection_tac_scenario_2_inputs <- make_fixed_projection_tac_inputs("scenario_2")
fixed_projection_tac_df <- fixed_projection_tac_inputs$fixed_projection_tac_df
fixed_catch_n_years <- fixed_projection_tac_inputs$fixed_catch_n_years
fixed_projection_tac <- fixed_projection_tac_inputs$fixed_projection_tac
projection_tac_split_base_f <- fixed_projection_tac_inputs$projection_tac_split_base_f
projection_tac_split_yf <- fixed_projection_tac_inputs$projection_tac_split_yf
ctp_tac_schedule <- 3L
ctp_tac_calculation_lag <- 2L
ctp_catch_data_lag <- 3L
ctp_cpue_data_lag <- 3L
ctp_gt_data_lag <- 4L
ctp_ckmr_data_lag <- 8L

gt_skip_years <- integer(0)
gt_projection_nrel <- 5000
gt_projection_nsam <- 10000
hsp_projection_n_per_cohort <- 1500
hsp_projection_nC <- hsp_projection_n_per_cohort^2
pop_projection_n_juvenile <- 1500
pop_projection_n_adult <- 1500
pop_projection_min_cohort <- 2002L
pop_projection_adult_ages <- 5:30

risk_relative_biomass_threshold <- 0.2
first_ctp_implementation_year <- max(fixed_projection_tac_df$year) + 1L
short_term_risk_years <- first_ctp_implementation_year:min(
  first_ctp_implementation_year + ctp_tac_schedule - 1L,
  projection_last_yr + 1L
)
long_term_risk_years <- seq.int(
  max(short_term_risk_years) + 1L,
  projection_last_yr + 1L
)

Projection Posterior

The projection simulations use 2,000 draws from the balanced 2,000-draw nine-cell MCMC grid posterior. The draw plan is reported in the grid page.

source(file.path(esc_dir, "esc31_workflow.R"))
source(file.path(esc_dir, "esc31_inputs.R"), local = TRUE)
source(file.path(esc_dir, "esc31_grid_contract.R"), local = TRUE)
source(file.path(esc_dir, "esc31_mle_grid_projection.R"), local = TRUE)
if (run_projection_computation ||
    run_no_uam_projection_computation ||
    run_mle_grid_projection_computation) {
  esc31_require_ll34_conditioning_approval(
    base_model_contract,
    stage = "projection computation"
  )
  esc31_require_grid_specification_approval(
    esc31_grid_specification,
    stage = "projection computation"
  )
}
expected_base_identity <- esc31_base_run_identity(
  esc_dir,
  base_model_contract = base_model_contract
)
base_fit_candidate <- if (file.exists(model_file)) {
  tryCatch(
    sbt_fit_read(
      model_file,
      strict = TRUE,
      rebuild = FALSE,
      verify_validation = FALSE
    ),
    error = function(e) NULL
  )
} else {
  NULL
}

if (refit_model || !esc31_base_fit_is_current(base_fit_candidate, expected_base_identity)) {
  stop(
    "The current base MLE is unavailable. Render 2_base.qmd explicitly before ",
    "running or reviewing the projection workflow.",
    call. = FALSE
  )
}

base_fit <- sbt_fit_read(
  model_file,
  strict = TRUE,
  rebuild = FALSE,
  verify_validation = FALSE
)
base_fit_md5 <- unname(tools::md5sum(model_file))
projection_harvest_wall_identity <- esc31_grid_harvest_wall_identity(
  base_model_contract
)
base_wall_data <- lapply(
  names(projection_harvest_wall_identity$data),
  function(name) base_fit$data[[name]]
)
names(base_wall_data) <- names(projection_harvest_wall_identity$data)
if (!identical(base_wall_data, projection_harvest_wall_identity$data) ||
    !identical(
      base_fit$provenance$metadata$run_identity$config$
        base_model_contract$harvest_wall,
      projection_harvest_wall_identity$decision
    )) {
  stop(
    "The current base fit does not carry the approved harvest-wall contract.",
    call. = FALSE
  )
}
grid_prerequisite_check <-
  esc31_validate_grid_prerequisite_acceptance(
    esc_dir = esc_dir,
    base_run_signature =
      base_fit$provenance$metadata$run_identity$signature,
    harvest_wall_identity = projection_harvest_wall_identity
  )
grid_prerequisite_validation <- grid_prerequisite_check$provenance

projection_map <- base_fit$map
projection_map$par_log_h <- NULL
projection_map$par_log_psi <- NULL
obj_projection <- MakeADFun(
  func = cmb(sbt_model, base_fit$data),
  parameters = base_fit$parameters,
  map = projection_map
)

mle_grid_projection_source <- esc31_load_mle_grid_projection_source(
  esc_dir = esc_dir,
  base_fit = base_fit,
  n_iter = projection_n_iter,
  seed = mle_grid_projection_draw_seed,
  state_cores = projection_state_diagnostic_cores
)
mle_grid_projection_fit <- mle_grid_projection_source$fit
mle_grid_projection_data <- mle_grid_projection_source$data
mle_grid_projection_object <- mle_grid_projection_source$object
mle_grid_projection_iters <- mle_grid_projection_source$iters
mle_grid_projection_post <- mle_grid_projection_source$post
mle_grid_projection_draw_metadata <-
  mle_grid_projection_source$metadata
mle_grid_projection_draw_plan <- mle_grid_projection_source$draw_plan
mle_grid_projection_historical_state <-
  mle_grid_projection_source$historical_state

if (!file.exists(no_uam_model_file)) {
  stop(
    "The accepted NoUAM fit is unavailable at ",
    no_uam_model_file,
    ". Render 3_sens.qmd before setting up its projection comparison.",
    call. = FALSE
  )
}
no_uam_fit <- sbt_fit_read(
  no_uam_model_file,
  strict = TRUE,
  rebuild = FALSE,
  verify_validation = FALSE
)
no_uam_metadata <- no_uam_fit$provenance$metadata
no_uam_validation <- sbt_fit_validation(
  no_uam_fit,
  scope = "mcmc",
  verify = FALSE,
  require_pass = TRUE
)
no_uam_diagnostics <- no_uam_fit$fit$diagnostics$mcmc
no_uam_contract_passes <-
  identical(no_uam_metadata$sensitivity_id, "noUAM") &&
  identical(no_uam_metadata$catch_UA_zeroed, TRUE) &&
  is.finite(no_uam_metadata$removed_total) &&
  no_uam_metadata$removed_total > 0 &&
  identical(
    no_uam_metadata$base_fit_scientific_signature,
    base_fit$provenance$metadata$
      mcmc_source_fit_implementation_signature
  ) &&
  identical(
    no_uam_metadata$base_run_signature,
    base_fit$provenance$metadata$run_identity$signature
  ) &&
  isTRUE(no_uam_validation$passes) &&
  is.data.frame(no_uam_diagnostics) &&
  nrow(no_uam_diagnostics) == 1L &&
  isTRUE(no_uam_diagnostics$passes) &&
  identical(as.integer(no_uam_diagnostics$chains), 4L) &&
  identical(as.integer(no_uam_diagnostics$samples_per_chain), 750L) &&
  identical(as.integer(no_uam_diagnostics$divergences), 0L) &&
  identical(as.integer(no_uam_diagnostics$max_treedepth_hits), 0L) &&
  isTRUE(no_uam_diagnostics$state_passes) &&
  identical(as.integer(no_uam_diagnostics$state_invalid_draws), 0L) &&
  isTRUE(all(is.na(no_uam_fit$map$par_log_h))) &&
  isTRUE(all(is.na(no_uam_fit$map$par_log_psi)))
if (!no_uam_contract_passes) {
  stop(
    "The saved NoUAM fit does not satisfy the accepted sensitivity and ",
    "posterior contract.",
    call. = FALSE
  )
}
no_uam_obj_projection <- MakeADFun(
  func = cmb(sbt_model, no_uam_fit$data),
  parameters = no_uam_fit$parameters,
  map = no_uam_fit$map
)
no_uam_mcmc <- as_tmbfit(no_uam_fit)

if (!file.exists(combined_mcmc_file)) {
  stop(
    "Could not find ", basename(combined_mcmc_file),
    ". Render 4_grid.qmd first so the balanced grid posterior is available.",
    call. = FALSE
  )
}

combined_saved <- load_rdata_cache(combined_mcmc_file)
required_combined_objects <- c(
  "projection_mcmc", "draw_metadata", "draw_plan", "combined_run_metadata",
  "combined_payload_checksum", "saved_source_signature"
)
if (is.null(combined_saved) || !all(vapply(
      required_combined_objects,
      exists,
      logical(1),
      envir = combined_saved,
      inherits = FALSE
    ))) {
  stop(
    "The combined posterior is missing required posterior or provenance objects. ",
    "Rerun 4_grid.qmd before projections.",
    call. = FALSE
  )
}
combined_run_metadata <- combined_saved$combined_run_metadata
saved_base_reuse <- combined_run_metadata$base_cell_reuse
current_base_posterior_payload_checksum <-
  sbt:::.grid_mcmc_posterior_payload_checksum(as_tmbfit(base_fit))
normalize_base_reuse_decision <- function(decision) {
  if (is.list(decision)) decision$source_posterior_file <- NULL
  decision
}
base_reuse_matches_current <-
  is.list(saved_base_reuse) &&
  identical(saved_base_reuse$schema_version, 1L) &&
  identical(
    normalize_base_reuse_decision(saved_base_reuse$decision),
    normalize_base_reuse_decision(esc31_grid_specification$base_cell_reuse)
  ) &&
  identical(
    as.integer(saved_base_reuse$cell),
    as.integer(esc31_grid_specification$coordinates$reused_base_cell)
  ) &&
  esc31_saved_fit_signature_matches(
    saved_base_reuse$source_fit_scientific_signature,
    base_fit
  ) &&
  identical(
    saved_base_reuse$source_mle_run_signature,
    base_fit$provenance$metadata$run_identity$signature
  ) &&
  identical(
    saved_base_reuse$source_mcmc_run_signature,
    base_fit$provenance$metadata$mcmc_run_identity$signature
  ) &&
  identical(
    saved_base_reuse$source_posterior_payload_checksum,
    current_base_posterior_payload_checksum
  ) &&
  identical(saved_base_reuse$source_sampler, base_fit$mcmc$settings$sampler) &&
  identical(
    saved_base_reuse$source_acceptance_thresholds,
    base_fit$mcmc$settings$acceptance_thresholds
  )

# The portable base fit can gain derived report summaries without changing its
# accepted posterior payload. Match the exact payload, sampler, thresholds, and
# run identities above before treating its changed whole-file checksum as
# audit-only, following the same rule used by 4_grid.qmd.
saved_prerequisite_for_comparison <-
  combined_run_metadata$prerequisite_acceptance
if (isTRUE(base_reuse_matches_current) &&
    is.list(saved_prerequisite_for_comparison) &&
    is.list(grid_prerequisite_validation) &&
    identical(names(saved_prerequisite_for_comparison$artifact_checksums), "base") &&
    identical(names(grid_prerequisite_validation$artifact_checksums), "base")) {
  saved_prerequisite_for_comparison$artifact_checksums <-
    grid_prerequisite_validation$artifact_checksums
}
expected_grid_samples_per_chain <- as.integer(unname(
  esc31_grid_mixed_source_selection$retained_draws_per_chain[
    as.character(seq_len(9L))
  ]
))
expected_grid_state_draws <- as.integer(
  esc31_grid_mixed_source_selection$chains *
    expected_grid_samples_per_chain
)
expected_combined_grid_sampler_contract <- list(
  chains = esc31_grid_mixed_source_selection$chains,
  num_warmup = esc31_grid_mixed_source_selection$warmup,
  num_samples_by_cell = expected_grid_samples_per_chain,
  metric = esc31_grid_mixed_source_selection$metric,
  adapt_delta = esc31_grid_mixed_source_selection$adapt_delta,
  max_treedepth = esc31_grid_mixed_source_selection$max_treedepth,
  init = esc31_grid_mixed_source_selection$init
)
combined_state_fields <- c(
  "Cell", "state_passes", "state_draws_expected", "state_draws_evaluated",
  "state_invalid_draws", "state_non_finite_draws", "state_invalid_cells",
  "state_max_raw_harvest", "state_min_number",
  "state_max_harvest_penalty", "state_checksum",
  "posterior_payload_checksum"
)
if (!is.list(combined_run_metadata) ||
    !is.list(combined_saved$projection_mcmc) ||
    !identical(combined_run_metadata$schema_version, 9L) ||
    !identical(combined_run_metadata$stage, "combined_grid_posterior") ||
    !identical(
      combined_run_metadata$grid_specification,
      esc31_grid_specification
    ) ||
    !esc31_grid_prerequisite_acceptance_equivalent(
      saved_prerequisite_for_comparison,
      grid_prerequisite_validation
    ) ||
    !is.list(combined_run_metadata$prerequisite_review_confirmation) ||
    !isTRUE(
      combined_run_metadata$prerequisite_review_confirmation$confirmed
    ) ||
    !identical(
      combined_run_metadata$prerequisite_review_confirmation$scope,
      "accepted_base_mcmc_only"
    ) ||
    !identical(
      combined_run_metadata$prerequisite_review_confirmation$mechanism,
      "esc31_base_only_grid_launch_v1"
    ) ||
    !identical(
      combined_run_metadata$implementation,
      "grid_mcmc_to_tmbfit_v11_mixed_source_reuse_accepted_base_cell5_harvest_wall_bound"
    ) ||
    !identical(
      combined_run_metadata$grid_mixed_source_selection,
      esc31_grid_mixed_source_selection
    ) ||
    !is.data.frame(combined_run_metadata$grid_source_artifacts) ||
    nrow(combined_run_metadata$grid_source_artifacts) != 9L ||
    !identical(
      as.integer(combined_run_metadata$grid_source_artifacts$Cell),
      seq_len(9L)
    ) ||
    !identical(
      as.integer(
        combined_run_metadata$grid_source_artifacts$
          retained_draws_per_chain
      ),
      expected_grid_samples_per_chain
    ) ||
    !identical(
      as.character(
        combined_run_metadata$grid_source_artifacts$
          source_grid_run_signature
      ),
      as.character(
        combined_run_metadata$grid_sampler_run_signature
      )
    ) ||
    !identical(
      as.character(
        combined_run_metadata$grid_source_artifacts$
          posterior_payload_checksum
      ),
      as.character(
        combined_run_metadata$biological_state_summary$
          posterior_payload_checksum
      )
    ) ||
    !identical(
      as.character(
        combined_run_metadata$grid_source_artifacts$state_checksum
      ),
      as.character(
        combined_run_metadata$biological_state_summary$state_checksum
      )
    ) ||
    !identical(
      combined_run_metadata$base_cell_reuse$decision,
      esc31_grid_specification$base_cell_reuse
    ) ||
    !identical(
      combined_run_metadata$grid_sampler_contract,
      expected_combined_grid_sampler_contract
    ) ||
    !is.character(combined_run_metadata$grid_sampler_run_signature) ||
    length(combined_run_metadata$grid_sampler_run_signature) != 9L ||
    anyNA(combined_run_metadata$grid_sampler_run_signature) ||
    any(!grepl(
      "^[0-9a-f]{32}$",
      combined_run_metadata$grid_sampler_run_signature
    )) ||
    !identical(
      combined_run_metadata$biological_state_contract,
      biological_state_contract()
    ) ||
    !identical(
      combined_run_metadata$harvest_wall_contract,
      projection_harvest_wall_identity$executable
    ) ||
    !identical(
      combined_run_metadata$harvest_wall_decision,
      projection_harvest_wall_identity$decision
    ) ||
    !identical(
      combined_run_metadata$harvest_wall_data,
      projection_harvest_wall_identity$data
    ) ||
    !identical(
      combined_run_metadata$harvest_wall_application,
      "fitted_grid_objective_conditioning"
    ) ||
    !is.data.frame(combined_run_metadata$biological_state_summary) ||
    nrow(combined_run_metadata$biological_state_summary) != 9L ||
    !all(combined_state_fields %in%
      names(combined_run_metadata$biological_state_summary)) ||
    !is.character(combined_run_metadata$biological_state_checksum) ||
    length(combined_run_metadata$biological_state_checksum) != 1L ||
    !grepl("^[0-9a-f]{32}$",
      combined_run_metadata$biological_state_checksum) ||
    !identical(
      combined_run_metadata$biological_state_checksum,
      esc31_object_md5(list(
        contract = combined_run_metadata$biological_state_contract,
        diagnostics = combined_run_metadata$biological_state_summary
      ))
    ) ||
    !identical(
      combined_run_metadata$biological_state_draws_per_cell,
      expected_grid_state_draws
    ) ||
    !identical(
      combined_run_metadata$mle_gradient_limit,
      esc31_grid_specification$optimizer$max_gradient
    ) ||
    !is.data.frame(combined_run_metadata$mle_gradient_summary) ||
    nrow(combined_run_metadata$mle_gradient_summary) != 9L ||
    !identical(
      names(combined_run_metadata$mle_gradient_summary),
      c("Cell", "max_gradient", "file_md5")
    ) ||
    !identical(
      as.integer(combined_run_metadata$mle_gradient_summary$Cell),
      seq_len(9L)
    ) ||
    any(!grepl(
      "^[0-9a-f]{32}$",
      combined_run_metadata$mle_gradient_summary$file_md5
    )) ||
    any(!is.finite(
      combined_run_metadata$mle_gradient_summary$max_gradient
    )) ||
    any(combined_run_metadata$mle_gradient_summary$max_gradient >
      combined_run_metadata$mle_gradient_limit) ||
    !is.character(combined_run_metadata$mle_gradient_checksum) ||
    length(combined_run_metadata$mle_gradient_checksum) != 1L ||
    !identical(
      combined_run_metadata$mle_gradient_checksum,
      esc31_object_md5(list(
        limit = combined_run_metadata$mle_gradient_limit,
        cells = combined_run_metadata$mle_gradient_summary
      ))
    ) ||
    any(c("diagnostic_exception", "provenance_exception") %in%
      names(combined_run_metadata)) ||
    !identical(as.integer(combined_run_metadata$draws), 2000L) ||
    !is.character(combined_run_metadata$grid_run_signature) ||
    length(combined_run_metadata$grid_run_signature) != 1L ||
    is.na(combined_run_metadata$grid_run_signature) ||
    !nzchar(combined_run_metadata$grid_run_signature)) {
  stop(
    "The combined posterior was not produced by the no-exceptions grid workflow. ",
    "Run a complete approved grid that passes every diagnostic and provenance check.",
    call. = FALSE
  )
}
projection_grid_prerequisite_provenance <- list(
  comparison =
    "validated_current_manifest_vs_manifest_used_to_build_combined_posterior",
  substantive_acceptance_identical =
    esc31_grid_prerequisite_acceptance_equivalent(
      saved_prerequisite_for_comparison,
      grid_prerequisite_validation
    ),
  base_posterior_payload_identical = base_reuse_matches_current,
  combined_base_file_md5 = combined_run_metadata$base_fit_md5,
  current_base_file_md5 = base_fit_md5,
  combined_posterior_substantive =
    combined_run_metadata$prerequisite_acceptance[
      setdiff(
        names(combined_run_metadata$prerequisite_acceptance),
        c(
          "manifest_signature",
          "validator_identity_signature",
          "base_posterior_source_checksum"
        )
      )
    ],
  current_validated_substantive =
    grid_prerequisite_validation[
      setdiff(
        names(grid_prerequisite_validation),
        c(
          "manifest_signature",
          "validator_identity_signature",
          "base_posterior_source_checksum"
        )
      )
    ]
)
combined_state_summary <- combined_run_metadata$biological_state_summary
combined_state_contract <- combined_run_metadata$biological_state_contract
if (anyNA(combined_state_summary[combined_state_fields]) ||
    !identical(sort(as.integer(combined_state_summary$Cell)), seq_len(9L)) ||
    any(!combined_state_summary$state_passes) ||
    any(combined_state_summary$state_draws_expected !=
      combined_run_metadata$biological_state_draws_per_cell) ||
    any(combined_state_summary$state_draws_evaluated !=
      combined_state_summary$state_draws_expected) ||
    any(combined_state_summary$state_invalid_draws != 0L) ||
    any(combined_state_summary$state_non_finite_draws != 0L) ||
    any(combined_state_summary$state_invalid_cells != 0L) ||
    any(!is.finite(combined_state_summary$state_max_raw_harvest)) ||
    any(combined_state_summary$state_max_raw_harvest >
      combined_state_contract$hrate_limit +
        combined_state_contract$hrate_tolerance) ||
    any(!is.finite(combined_state_summary$state_min_number)) ||
    any(combined_state_summary$state_min_number <=
      combined_state_contract$number_tolerance) ||
    any(!is.finite(combined_state_summary$state_max_harvest_penalty)) ||
    any(abs(combined_state_summary$state_max_harvest_penalty) >
      combined_state_contract$penalty_tolerance) ||
    any(!grepl("^[0-9a-f]{32}$", combined_state_summary$state_checksum)) ||
    any(!grepl(
      "^[0-9a-f]{32}$",
      combined_state_summary$posterior_payload_checksum
    ))) {
  stop(
    "The combined posterior does not have complete passing biological-state ",
    "diagnostics for all nine source cells.",
    call. = FALSE
  )
}
if (!isTRUE(base_reuse_matches_current) ||
    !identical(
      combined_run_metadata$base_run_signature,
      base_fit$provenance$metadata$run_identity$signature
    )) {
  stop(
    "The combined posterior does not match the current base posterior payload ",
    "and inputs. ",
    "Rerun 4_grid.qmd before projections.",
    call. = FALSE
  )
}
projection_mcmc <- combined_saved$projection_mcmc
draw_metadata <- combined_saved$draw_metadata
draw_plan <- combined_saved$draw_plan
saved_source_signature <- combined_saved$saved_source_signature
saved_payload_checksum <- combined_saved$combined_payload_checksum
if (!is.array(projection_mcmc$samples) ||
    length(dim(projection_mcmc$samples)) != 3L ||
    dim(projection_mcmc$samples)[2L] != 1L) {
  stop("The combined posterior must contain one canonical chain.", call. = FALSE)
}
combined_projection_post <- extract_samples(projection_mcmc)
recalculated_payload_checksum <- esc31_object_md5(list(
  samples = projection_mcmc$samples,
  warmup = projection_mcmc$warmup,
  sample_names = projection_mcmc$sample_names,
  draw_metadata = as.data.frame(draw_metadata),
  draw_plan = as.data.frame(draw_plan)
))
recalculated_source_signature <- esc31_object_md5(list(
  metadata_signature = esc31_object_md5(combined_run_metadata),
  payload_checksum = recalculated_payload_checksum
))
if (!identical(saved_payload_checksum, recalculated_payload_checksum) ||
    !identical(saved_source_signature, recalculated_source_signature)) {
  stop(
    "The combined posterior payload or provenance signature is invalid.",
    call. = FALSE
  )
}
if (!is.data.frame(draw_metadata) || !is.data.frame(draw_plan) ||
    nrow(draw_metadata) != 2000L || nrow(combined_projection_post) != 2000L ||
    nrow(draw_plan) != 9L ||
    !all(c("draw", "cell", "h", "psi") %in% names(draw_metadata)) ||
    !all(c("cell", "h", "psi", "draws") %in% names(draw_plan)) ||
    anyNA(draw_metadata[c("draw", "cell", "h", "psi")]) ||
    anyNA(draw_plan[c("cell", "h", "psi", "draws")]) ||
    anyDuplicated(draw_plan$cell) ||
    !identical(sort(as.integer(draw_plan$cell)), seq_len(9L))) {
  stop("The combined posterior does not contain the required 2,000-draw plan.", call. = FALSE)
}
expected_projection_grid <- expand.grid(
  h = c(0.6, 0.7, 0.8),
  psi = c(1.5, 1.75, 2)
)
ordered_draw_plan <- draw_plan[match(seq_len(9L), draw_plan$cell), , drop = FALSE]
if (any(!draw_metadata$cell %in% seq_len(9L)) ||
    !isTRUE(all.equal(
      as.numeric(ordered_draw_plan$h), expected_projection_grid$h,
      tolerance = 1e-12
    )) ||
    !isTRUE(all.equal(
      as.numeric(ordered_draw_plan$psi), expected_projection_grid$psi,
      tolerance = 1e-12
    )) ||
    any(abs(draw_metadata$h - expected_projection_grid$h[draw_metadata$cell]) > 1e-12) ||
    any(abs(draw_metadata$psi - expected_projection_grid$psi[draw_metadata$cell]) > 1e-12)) {
  stop("Combined-posterior cells do not match the accepted h/psi grid.", call. = FALSE)
}
ordered_draw_metadata <- draw_metadata[
  order(as.integer(draw_metadata$draw)), , drop = FALSE
]
combined_projection_names <- colnames(combined_projection_post)
if (!all(c("par_log_h", "par_log_psi") %in%
      combined_projection_names) ||
    !identical(
      as.integer(ordered_draw_metadata$draw),
      seq_len(nrow(combined_projection_post))
    ) ||
    !isTRUE(all.equal(
      exp(as.numeric(combined_projection_post[, "par_log_h"])),
      as.numeric(ordered_draw_metadata$h),
      tolerance = 1e-12
    )) ||
    !isTRUE(all.equal(
      exp(as.numeric(combined_projection_post[, "par_log_psi"])),
      as.numeric(ordered_draw_metadata$psi),
      tolerance = 1e-12
    ))) {
  stop(
    "Combined-posterior rows are not aligned with their h/psi draw metadata.",
    call. = FALSE
  )
}
metadata_cell_counts <- as.integer(table(factor(
  draw_metadata$cell,
  levels = seq_len(9L)
)))
planned_cell_counts <- draw_plan$draws[match(seq_len(9L), draw_plan$cell)]
if (anyNA(planned_cell_counts) ||
    !identical(metadata_cell_counts, as.integer(planned_cell_counts))) {
  stop("Combined-posterior cell counts do not match the saved draw plan.", call. = FALSE)
}
if (projection_n_iter > nrow(combined_projection_post)) {
  stop("Requested ", projection_n_iter, " projection draws but found only ", nrow(combined_projection_post), call. = FALSE)
}
if (sum(draw_plan$draws) != projection_mcmc_source_draws) {
  stop("Balanced grid draw plan does not sum to ", projection_mcmc_source_draws, call. = FALSE)
}

projection_iters <- sbt::select_balanced_grid_iters(
  draw_metadata,
  n_draws = projection_n_iter,
  seed = projection_draw_seed,
  expected_cells = seq_len(9L)
)
projection_draw_metadata <- as_tibble(draw_metadata) |>
  slice(match(projection_iters, .data$draw)) |>
  mutate(projection_draw = row_number(), .before = 1)
projection_draw_plan <- projection_draw_metadata |>
  group_by(.data$cell, .data$h, .data$psi) |>
  summarise(
    draws = n(),
    source_rows = paste(sort(.data$draw), collapse = ", "),
    .groups = "drop"
  ) |>
  arrange(.data$cell)

if (length(projection_iters) != projection_n_iter || anyDuplicated(projection_iters)) {
  stop("Projection row selection did not return the requested unique draws.", call. = FALSE)
}
if (!identical(as.integer(projection_draw_plan$cell), seq_len(9L)) ||
    diff(range(projection_draw_plan$draws)) > 1L) {
  stop("Projection draws are not balanced across all nine grid cells.", call. = FALSE)
}
if (projection_n_iter == projection_mcmc_source_draws &&
    (!identical(
       sort(projection_iters),
       seq_len(projection_mcmc_source_draws)
     ) ||
     !identical(
       as.integer(projection_draw_plan$draws),
       as.integer(draw_plan$draws)
     ))) {
  stop(
    "The 2,000-draw production projection plan must use every combined ",
    "posterior row exactly once and preserve its nine-cell allocation.",
    call. = FALSE
  )
}
projection_post <- combined_projection_post[projection_iters, , drop = FALSE]
projection_historical_state <- diagnose_sbt_states(
  obj_projection,
  data = base_fit$data,
  posterior = projection_mcmc,
  iters = projection_iters,
  chains = 1L,
  cores = projection_state_diagnostic_cores
)
if (!esc31_biological_state_diagnostics_pass(
      projection_historical_state, projection_n_iter
    )) {
  stop(
    "At least one selected projection draw has an invalid historical ",
    "biological state.",
    call. = FALSE
  )
}
projection_workflow_identity <- esc31_workflow_identity(
  esc_dir,
  workflow_files = "esc31_grid_contract.R",
  sbt_entry_points = esc31_sbt_entry_points("projections"),
  config = list(
    stage = "projections",
    grid_specification = esc31_grid_specification,
    base_fit_md5 = base_fit_md5,
    base_run_signature =
      base_fit$provenance$metadata$run_identity$signature,
    combined_file_md5 = unname(tools::md5sum(combined_mcmc_file)),
    combined_source_signature = saved_source_signature,
    combined_payload_checksum = saved_payload_checksum,
    projection_iters = projection_iters,
    projection_draw_seed = projection_draw_seed,
    projection_draw_plan = projection_draw_plan,
    grid_prerequisite_provenance =
      projection_grid_prerequisite_provenance,
    historical_state_checksum =
      projection_historical_state$summary$checksum,
    biological_state_contract = biological_state_contract(),
    harvest_wall_contract = projection_harvest_wall_identity$executable,
    harvest_wall_decision = projection_harvest_wall_identity$decision,
    harvest_wall_application =
      "source_posterior_conditioning_only_not_projection_dynamics",
    draw_selection = "balanced_nine_cell_seeded_v1",
    seed = projection_seed,
    projection_workflow_version =
      "esc31_projections_v8_terminal_selectivity_nominal_tac_removals"
  )
)
projection_source_signature <- esc31_object_md5(list(
  workflow_signature = projection_workflow_identity$scientific_signature,
  base_fit_md5 = base_fit_md5,
  base_run_signature = base_fit$provenance$metadata$run_identity$signature,
  combined_file_md5 = unname(tools::md5sum(combined_mcmc_file)),
  combined_source_signature = saved_source_signature,
  combined_payload_checksum = saved_payload_checksum,
  combined_run_metadata = combined_run_metadata,
  harvest_wall_contract = projection_harvest_wall_identity$executable,
  harvest_wall_decision = projection_harvest_wall_identity$decision,
  harvest_wall_application =
    "source_posterior_conditioning_only_not_projection_dynamics",
  draw_plan = draw_plan,
  draw_metadata = draw_metadata,
  projection_draw_plan = projection_draw_plan,
  projection_draw_metadata = projection_draw_metadata,
  projection_iters = projection_iters,
  grid_prerequisite_provenance =
    projection_grid_prerequisite_provenance,
  selected_posterior_md5 = esc31_object_md5(projection_post),
  historical_state_diagnostics = projection_historical_state$summary
))

mle_grid_projection_workflow_identity <- esc31_workflow_identity(
  esc_dir,
  workflow_files = c(
    "esc31_grid_contract.R",
    "esc31_mle_grid_projection.R"
  ),
  sbt_entry_points = esc31_sbt_entry_points("projections"),
  config = list(
    stage = "sampled_mle_grid_projections",
    base_projection_workflow_signature =
      projection_workflow_identity$scientific_signature,
    base_fit_md5 = base_fit_md5,
    grid_file_md5 = mle_grid_projection_source$file_md5,
    source_signature = mle_grid_projection_source$source_signature,
    projection_iters = mle_grid_projection_iters,
    draw_seed = mle_grid_projection_draw_seed,
    draw_plan = mle_grid_projection_draw_plan,
    M_switch = mle_grid_projection_data$M_switch,
    mortality_function = "get_M",
    seed = projection_seed,
    implementation =
      "esc31_report_specific_sampled_mle_grid_projection_v1"
  )
)
mle_grid_projection_source_signature <- esc31_object_md5(list(
  workflow_signature =
    mle_grid_projection_workflow_identity$scientific_signature,
  source_signature = mle_grid_projection_source$source_signature,
  grid_file_md5 = mle_grid_projection_source$file_md5,
  selected_posterior_md5 =
    esc31_object_md5(mle_grid_projection_post),
  draw_metadata = mle_grid_projection_draw_metadata,
  draw_plan = mle_grid_projection_draw_plan,
  historical_state_diagnostics =
    mle_grid_projection_historical_state$summary
))

no_uam_samples <- no_uam_mcmc$samples
no_uam_warmup <- as.integer(no_uam_mcmc$warmup)
if (!is.array(no_uam_samples) || length(dim(no_uam_samples)) != 3L ||
    dim(no_uam_samples)[2L] != 4L ||
    no_uam_warmup < 0L || no_uam_warmup >= dim(no_uam_samples)[1L] ||
    projection_n_iter %% dim(no_uam_samples)[2L] != 0L) {
  stop(
    "The NoUAM posterior must contain four equal retained chains and the ",
    "projection draw count must be divisible by four.",
    call. = FALSE
  )
}
no_uam_retained_per_chain <- dim(no_uam_samples)[1L] - no_uam_warmup
no_uam_draws_per_chain <- projection_n_iter %/% dim(no_uam_samples)[2L]
set.seed(no_uam_projection_draw_seed)
no_uam_iteration_subset <- sort(sample(
  seq_len(no_uam_retained_per_chain),
  no_uam_draws_per_chain,
  replace = FALSE
))
no_uam_chains <- seq_len(dim(no_uam_samples)[2L])
no_uam_projection_iters <- unlist(lapply(no_uam_chains, function(chain) {
  (chain - 1L) * no_uam_retained_per_chain + no_uam_iteration_subset
}), use.names = FALSE)
no_uam_projection_draw_plan <- crossing(
  chain = no_uam_chains,
  retained_iteration = no_uam_iteration_subset
) |>
  arrange(.data$chain, .data$retained_iteration) |>
  mutate(
    source_row =
      (.data$chain - 1L) * no_uam_retained_per_chain +
      .data$retained_iteration,
    projection_draw = row_number(),
    .before = 1
  )
if (!identical(no_uam_projection_draw_plan$source_row, no_uam_projection_iters) ||
    length(no_uam_projection_iters) != projection_n_iter ||
    anyDuplicated(no_uam_projection_iters)) {
  stop("The NoUAM chain-balanced projection draw plan is invalid.", call. = FALSE)
}
no_uam_projection_post_all <- SparseNUTS::extract_samples(no_uam_mcmc)
no_uam_projection_post <- no_uam_projection_post_all[
  no_uam_projection_iters, , drop = FALSE
]
no_uam_projection_historical_state <- diagnose_sbt_states(
  no_uam_obj_projection,
  data = no_uam_fit$data,
  posterior = no_uam_mcmc,
  iters = no_uam_iteration_subset,
  chains = no_uam_chains,
  cores = projection_state_diagnostic_cores
)
if (!esc31_biological_state_diagnostics_pass(
      no_uam_projection_historical_state, projection_n_iter
    )) {
  stop(
    "At least one selected NoUAM projection draw has an invalid historical ",
    "biological state.",
    call. = FALSE
  )
}
no_uam_projection_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  workflow_signature = projection_workflow_identity$scientific_signature,
  fit_md5 = unname(tools::md5sum(no_uam_model_file)),
  fit_metadata = no_uam_metadata,
  mcmc_validation_id = no_uam_validation$id,
  mcmc_diagnostics = no_uam_diagnostics,
  draw_seed = no_uam_projection_draw_seed,
  draw_plan = no_uam_projection_draw_plan,
  selected_posterior_md5 = esc31_object_md5(no_uam_projection_post),
  historical_state_diagnostics =
    no_uam_projection_historical_state$summary,
  implementation = "no_uam_four_chain_balanced_projection_source_v1"
))

# Upstream projection inputs have component-specific identities. In
# particular, recruitment and selectivity are not invalidated by changes to
# downstream CTP catch staging in run_projections().
projection_recdev_code_identity <- esc31_sbt_function_identity(
  "project_rec_devs"
)
projection_selectivity_code_identity <- esc31_sbt_function_identity(c(
  "sbt_model", "sbt_obj", "project_selectivity"
))
posterior_component_md5 <- function(posterior, pattern) {
  columns <- grepl(pattern, colnames(posterior))
  if (!any(columns)) {
    stop("The selected posterior has no columns matching: ", pattern,
         call. = FALSE)
  }
  esc31_object_md5(posterior[, columns, drop = FALSE])
}

projection_recdev_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  model_fit_md5 = base_fit_md5,
  posterior_file_md5 = unname(tools::md5sum(combined_mcmc_file)),
  posterior_source_signature = saved_source_signature,
  posterior_payload_checksum = saved_payload_checksum,
  selected_rdev_md5 = posterior_component_md5(
    projection_post, "par_rdev_y"
  ),
  data_years = c(base_fit$data$first_yr, base_fit$data$last_yr),
  code_identity = projection_recdev_code_identity$signature
))
projection_selectivity_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  model_fit_md5 = base_fit_md5,
  posterior_file_md5 = unname(tools::md5sum(combined_mcmc_file)),
  posterior_source_signature = saved_source_signature,
  posterior_payload_checksum = saved_payload_checksum,
  selected_posterior_md5 = esc31_object_md5(projection_post),
  model_data_md5 = esc31_object_md5(base_fit$data),
  code_identity = projection_selectivity_code_identity$signature
))
projection_cpue_q_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  model_fit_md5 = base_fit_md5,
  posterior_file_md5 = unname(tools::md5sum(combined_mcmc_file)),
  selected_cpue_parameter_md5 = posterior_component_md5(
    projection_post, "par_log_cpue_q|par_cpue_creep"
  ),
  cpue_years = base_fit$data$cpue_years,
  first_yr = base_fit$data$first_yr,
  implementation = "esc31_cpue_q_v3"
))

no_uam_projection_recdev_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  model_fit_md5 = unname(tools::md5sum(no_uam_model_file)),
  mcmc_validation_id = no_uam_validation$id,
  draw_plan = no_uam_projection_draw_plan,
  selected_rdev_md5 = posterior_component_md5(
    no_uam_projection_post, "par_rdev_y"
  ),
  data_years = c(no_uam_fit$data$first_yr, no_uam_fit$data$last_yr),
  code_identity = projection_recdev_code_identity$signature
))
no_uam_projection_selectivity_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  model_fit_md5 = unname(tools::md5sum(no_uam_model_file)),
  mcmc_validation_id = no_uam_validation$id,
  draw_plan = no_uam_projection_draw_plan,
  selected_posterior_md5 = esc31_object_md5(no_uam_projection_post),
  model_data_md5 = esc31_object_md5(no_uam_fit$data),
  code_identity = projection_selectivity_code_identity$signature
))

mle_grid_projection_recdev_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  source_signature = mle_grid_projection_source$source_signature,
  grid_file_md5 = mle_grid_projection_source$file_md5,
  draw_plan = mle_grid_projection_draw_plan,
  selected_rdev_md5 = posterior_component_md5(
    mle_grid_projection_post, "par_rdev_y"
  ),
  data_years = c(
    mle_grid_projection_data$first_yr,
    mle_grid_projection_data$last_yr
  ),
  code_identity = projection_recdev_code_identity$signature
))
mle_grid_projection_selectivity_source_signature <- esc31_object_md5(list(
  schema_version = 1L,
  source_signature = mle_grid_projection_source$source_signature,
  grid_file_md5 = mle_grid_projection_source$file_md5,
  draw_plan = mle_grid_projection_draw_plan,
  selected_posterior_md5 = esc31_object_md5(mle_grid_projection_post),
  model_data_md5 = esc31_object_md5(mle_grid_projection_data),
  code_identity = projection_selectivity_code_identity$signature
))

projection_draw_plan |>
  transmute(
    Cell = .data$cell,
    h = formatC(.data$h, format = "f", digits = 2),
    psi = formatC(.data$psi, format = "f", digits = 2),
    Draws = .data$draws,
    `Combined-posterior source rows` = .data$source_rows
  ) |>
  kable()

Table 1: Balanced projection draw set by grid cell. Source rows are the retained indices in the combined posterior.

Cell	h	psi	Draws	Combined-posterior source rows
1	0.60	1.50	222	1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222
2	0.70	1.50	222	223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444
3	0.80	1.50	222	445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666
4	0.60	1.75	222	667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888
5	0.70	1.75	222	889, 890, 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024, 1025, 1026, 1027, 1028, 1029, 1030, 1031, 1032, 1033, 1034, 1035, 1036, 1037, 1038, 1039, 1040, 1041, 1042, 1043, 1044, 1045, 1046, 1047, 1048, 1049, 1050, 1051, 1052, 1053, 1054, 1055, 1056, 1057, 1058, 1059, 1060, 1061, 1062, 1063, 1064, 1065, 1066, 1067, 1068, 1069, 1070, 1071, 1072, 1073, 1074, 1075, 1076, 1077, 1078, 1079, 1080, 1081, 1082, 1083, 1084, 1085, 1086, 1087, 1088, 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108, 1109, 1110
6	0.80	1.75	223	1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118, 1119, 1120, 1121, 1122, 1123, 1124, 1125, 1126, 1127, 1128, 1129, 1130, 1131, 1132, 1133, 1134, 1135, 1136, 1137, 1138, 1139, 1140, 1141, 1142, 1143, 1144, 1145, 1146, 1147, 1148, 1149, 1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170, 1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182, 1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194, 1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, 1203, 1204, 1205, 1206, 1207, 1208, 1209, 1210, 1211, 1212, 1213, 1214, 1215, 1216, 1217, 1218, 1219, 1220, 1221, 1222, 1223, 1224, 1225, 1226, 1227, 1228, 1229, 1230, 1231, 1232, 1233, 1234, 1235, 1236, 1237, 1238, 1239, 1240, 1241, 1242, 1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255, 1256, 1257, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266, 1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1276, 1277, 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, 1288, 1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297, 1298, 1299, 1300, 1301, 1302, 1303, 1304, 1305, 1306, 1307, 1308, 1309, 1310, 1311, 1312, 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320, 1321, 1322, 1323, 1324, 1325, 1326, 1327, 1328, 1329, 1330, 1331, 1332, 1333
7	0.60	2.00	223	1334, 1335, 1336, 1337, 1338, 1339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348, 1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358, 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368, 1369, 1370, 1371, 1372, 1373, 1374, 1375, 1376, 1377, 1378, 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1386, 1387, 1388, 1389, 1390, 1391, 1392, 1393, 1394, 1395, 1396, 1397, 1398, 1399, 1400, 1401, 1402, 1403, 1404, 1405, 1406, 1407, 1408, 1409, 1410, 1411, 1412, 1413, 1414, 1415, 1416, 1417, 1418, 1419, 1420, 1421, 1422, 1423, 1424, 1425, 1426, 1427, 1428, 1429, 1430, 1431, 1432, 1433, 1434, 1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453, 1454, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462, 1463, 1464, 1465, 1466, 1467, 1468, 1469, 1470, 1471, 1472, 1473, 1474, 1475, 1476, 1477, 1478, 1479, 1480, 1481, 1482, 1483, 1484, 1485, 1486, 1487, 1488, 1489, 1490, 1491, 1492, 1493, 1494, 1495, 1496, 1497, 1498, 1499, 1500, 1501, 1502, 1503, 1504, 1505, 1506, 1507, 1508, 1509, 1510, 1511, 1512, 1513, 1514, 1515, 1516, 1517, 1518, 1519, 1520, 1521, 1522, 1523, 1524, 1525, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1533, 1534, 1535, 1536, 1537, 1538, 1539, 1540, 1541, 1542, 1543, 1544, 1545, 1546, 1547, 1548, 1549, 1550, 1551, 1552, 1553, 1554, 1555, 1556
8	0.70	2.00	222	1557, 1558, 1559, 1560, 1561, 1562, 1563, 1564, 1565, 1566, 1567, 1568, 1569, 1570, 1571, 1572, 1573, 1574, 1575, 1576, 1577, 1578, 1579, 1580, 1581, 1582, 1583, 1584, 1585, 1586, 1587, 1588, 1589, 1590, 1591, 1592, 1593, 1594, 1595, 1596, 1597, 1598, 1599, 1600, 1601, 1602, 1603, 1604, 1605, 1606, 1607, 1608, 1609, 1610, 1611, 1612, 1613, 1614, 1615, 1616, 1617, 1618, 1619, 1620, 1621, 1622, 1623, 1624, 1625, 1626, 1627, 1628, 1629, 1630, 1631, 1632, 1633, 1634, 1635, 1636, 1637, 1638, 1639, 1640, 1641, 1642, 1643, 1644, 1645, 1646, 1647, 1648, 1649, 1650, 1651, 1652, 1653, 1654, 1655, 1656, 1657, 1658, 1659, 1660, 1661, 1662, 1663, 1664, 1665, 1666, 1667, 1668, 1669, 1670, 1671, 1672, 1673, 1674, 1675, 1676, 1677, 1678, 1679, 1680, 1681, 1682, 1683, 1684, 1685, 1686, 1687, 1688, 1689, 1690, 1691, 1692, 1693, 1694, 1695, 1696, 1697, 1698, 1699, 1700, 1701, 1702, 1703, 1704, 1705, 1706, 1707, 1708, 1709, 1710, 1711, 1712, 1713, 1714, 1715, 1716, 1717, 1718, 1719, 1720, 1721, 1722, 1723, 1724, 1725, 1726, 1727, 1728, 1729, 1730, 1731, 1732, 1733, 1734, 1735, 1736, 1737, 1738, 1739, 1740, 1741, 1742, 1743, 1744, 1745, 1746, 1747, 1748, 1749, 1750, 1751, 1752, 1753, 1754, 1755, 1756, 1757, 1758, 1759, 1760, 1761, 1762, 1763, 1764, 1765, 1766, 1767, 1768, 1769, 1770, 1771, 1772, 1773, 1774, 1775, 1776, 1777, 1778
9	0.80	2.00	222	1779, 1780, 1781, 1782, 1783, 1784, 1785, 1786, 1787, 1788, 1789, 1790, 1791, 1792, 1793, 1794, 1795, 1796, 1797, 1798, 1799, 1800, 1801, 1802, 1803, 1804, 1805, 1806, 1807, 1808, 1809, 1810, 1811, 1812, 1813, 1814, 1815, 1816, 1817, 1818, 1819, 1820, 1821, 1822, 1823, 1824, 1825, 1826, 1827, 1828, 1829, 1830, 1831, 1832, 1833, 1834, 1835, 1836, 1837, 1838, 1839, 1840, 1841, 1842, 1843, 1844, 1845, 1846, 1847, 1848, 1849, 1850, 1851, 1852, 1853, 1854, 1855, 1856, 1857, 1858, 1859, 1860, 1861, 1862, 1863, 1864, 1865, 1866, 1867, 1868, 1869, 1870, 1871, 1872, 1873, 1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883, 1884, 1885, 1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894, 1895, 1896, 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1914, 1915, 1916, 1917, 1918, 1919, 1920, 1921, 1922, 1923, 1924, 1925, 1926, 1927, 1928, 1929, 1930, 1931, 1932, 1933, 1934, 1935, 1936, 1937, 1938, 1939, 1940, 1941, 1942, 1943, 1944, 1945, 1946, 1947, 1948, 1949, 1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958, 1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000

projection_historical_state$summary |>
  transmute(
    `Draws expected` = draws_expected,
    `Draws checked` = draws_evaluated,
    `Invalid draws` = invalid_draws,
    `Non-finite draws` = non_finite_draws,
    `Invalid state cells` = invalid_state_cells,
    `Maximum raw harvest` = max_raw_harvest_rate,
    `Minimum abundance` = min_number,
    `Maximum continuation penalty` = max_harvest_penalty
  ) |>
  mutate(across(
    c(
      `Draws expected`,
      `Draws checked`,
      `Invalid draws`,
      `Non-finite draws`,
      `Invalid state cells`
    ),
    ~ formatC(as.integer(.x), format = "d", big.mark = ",")
  )) |>
  mutate(across(where(is.numeric), ~ format_decimal(.x, digits = 4))) |>
  kable()

Table 2: Independent biological-state audit of the selected historical posterior draws.

Draws expected	Draws checked	Invalid draws	Non-finite draws	Invalid state cells	Maximum raw harvest	Minimum abundance	Maximum continuation penalty
2,000	2,000	0	0	0	0.8593	575.2417	0.0000

NoUAM comparison posterior

The NoUAM comparison uses the accepted sensitivity posterior in which the NCNM catch additions were removed before fitting. The same removal is now carried into its projections: unlike the base arm, its future LL1 nominal TAC allocation is not multiplied by 1.11. Both arms retain the separate 1.20 Australian surface-fishery correction. They otherwise use the same projection years, fixed nominal TAC totals and allocations, observation programme, CTP schedule, and projection seed. The NoUAM posterior, historical state, recruitment deviations, and selectivity remain separate from the nine-cell base posterior.

no_uam_projection_draw_plan |>
  group_by(.data$chain) |>
  summarise(
    Draws = n(),
    `Retained iterations` =
      paste(.data$retained_iteration, collapse = ", "),
    .groups = "drop"
  ) |>
  rename(Chain = chain) |>
  kable()

Table 3: Chain-balanced NoUAM projection subset.

Chain	Draws	Retained iterations
1	500	1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 20, 21, 22, 23, 25, 27, 28, 29, 30, 31, 35, 37, 38, 39, 40, 42, 43, 44, 45, 46, 50, 51, 52, 53, 55, 56, 58, 62, 63, 66, 69, 70, 71, 72, 73, 74, 75, 76, 77, 79, 80, 81, 83, 84, 85, 87, 88, 89, 90, 91, 92, 93, 94, 95, 98, 99, 102, 105, 106, 108, 110, 111, 112, 113, 117, 118, 119, 120, 121, 122, 124, 125, 126, 127, 128, 129, 130, 133, 134, 136, 138, 139, 141, 142, 143, 146, 148, 149, 150, 151, 152, 154, 156, 159, 160, 161, 163, 164, 166, 167, 168, 169, 170, 171, 173, 174, 175, 176, 179, 180, 182, 183, 184, 187, 190, 191, 192, 194, 195, 196, 197, 198, 203, 204, 205, 206, 207, 208, 209, 211, 213, 214, 215, 217, 219, 221, 222, 223, 225, 227, 228, 229, 230, 231, 234, 235, 236, 239, 241, 244, 245, 246, 248, 250, 252, 253, 255, 256, 257, 258, 260, 261, 263, 265, 268, 271, 272, 273, 274, 275, 276, 277, 282, 288, 290, 291, 292, 294, 295, 296, 300, 301, 303, 304, 305, 307, 308, 310, 311, 312, 314, 315, 316, 317, 319, 320, 322, 323, 324, 325, 327, 328, 331, 332, 334, 335, 336, 337, 340, 341, 344, 345, 346, 347, 348, 349, 353, 354, 356, 362, 363, 366, 367, 369, 370, 371, 372, 374, 375, 376, 379, 382, 383, 384, 385, 386, 389, 390, 391, 392, 394, 396, 397, 398, 399, 401, 402, 403, 404, 405, 406, 407, 410, 412, 414, 418, 419, 420, 422, 423, 424, 426, 427, 428, 430, 431, 432, 433, 435, 437, 440, 441, 443, 446, 447, 448, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 463, 468, 469, 470, 471, 473, 474, 475, 476, 477, 478, 479, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 496, 499, 500, 501, 502, 503, 504, 505, 508, 509, 512, 513, 515, 516, 517, 519, 520, 521, 522, 523, 524, 525, 526, 530, 531, 532, 533, 536, 537, 538, 539, 542, 543, 544, 545, 547, 549, 550, 551, 552, 555, 557, 560, 561, 564, 565, 567, 568, 569, 570, 571, 573, 575, 576, 577, 578, 581, 582, 583, 586, 587, 588, 590, 592, 593, 594, 595, 596, 597, 601, 602, 603, 604, 606, 607, 609, 610, 612, 614, 615, 616, 617, 618, 619, 622, 623, 625, 627, 628, 629, 630, 631, 632, 633, 634, 637, 642, 643, 644, 645, 650, 652, 653, 654, 656, 657, 658, 659, 660, 663, 664, 665, 666, 667, 668, 670, 671, 673, 676, 677, 678, 680, 682, 684, 685, 686, 688, 689, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 717, 719, 720, 721, 723, 724, 725, 726, 727, 728, 729, 730, 736, 738, 739, 741, 742, 743, 744, 745, 748, 749, 750
2	500	1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 20, 21, 22, 23, 25, 27, 28, 29, 30, 31, 35, 37, 38, 39, 40, 42, 43, 44, 45, 46, 50, 51, 52, 53, 55, 56, 58, 62, 63, 66, 69, 70, 71, 72, 73, 74, 75, 76, 77, 79, 80, 81, 83, 84, 85, 87, 88, 89, 90, 91, 92, 93, 94, 95, 98, 99, 102, 105, 106, 108, 110, 111, 112, 113, 117, 118, 119, 120, 121, 122, 124, 125, 126, 127, 128, 129, 130, 133, 134, 136, 138, 139, 141, 142, 143, 146, 148, 149, 150, 151, 152, 154, 156, 159, 160, 161, 163, 164, 166, 167, 168, 169, 170, 171, 173, 174, 175, 176, 179, 180, 182, 183, 184, 187, 190, 191, 192, 194, 195, 196, 197, 198, 203, 204, 205, 206, 207, 208, 209, 211, 213, 214, 215, 217, 219, 221, 222, 223, 225, 227, 228, 229, 230, 231, 234, 235, 236, 239, 241, 244, 245, 246, 248, 250, 252, 253, 255, 256, 257, 258, 260, 261, 263, 265, 268, 271, 272, 273, 274, 275, 276, 277, 282, 288, 290, 291, 292, 294, 295, 296, 300, 301, 303, 304, 305, 307, 308, 310, 311, 312, 314, 315, 316, 317, 319, 320, 322, 323, 324, 325, 327, 328, 331, 332, 334, 335, 336, 337, 340, 341, 344, 345, 346, 347, 348, 349, 353, 354, 356, 362, 363, 366, 367, 369, 370, 371, 372, 374, 375, 376, 379, 382, 383, 384, 385, 386, 389, 390, 391, 392, 394, 396, 397, 398, 399, 401, 402, 403, 404, 405, 406, 407, 410, 412, 414, 418, 419, 420, 422, 423, 424, 426, 427, 428, 430, 431, 432, 433, 435, 437, 440, 441, 443, 446, 447, 448, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 463, 468, 469, 470, 471, 473, 474, 475, 476, 477, 478, 479, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 496, 499, 500, 501, 502, 503, 504, 505, 508, 509, 512, 513, 515, 516, 517, 519, 520, 521, 522, 523, 524, 525, 526, 530, 531, 532, 533, 536, 537, 538, 539, 542, 543, 544, 545, 547, 549, 550, 551, 552, 555, 557, 560, 561, 564, 565, 567, 568, 569, 570, 571, 573, 575, 576, 577, 578, 581, 582, 583, 586, 587, 588, 590, 592, 593, 594, 595, 596, 597, 601, 602, 603, 604, 606, 607, 609, 610, 612, 614, 615, 616, 617, 618, 619, 622, 623, 625, 627, 628, 629, 630, 631, 632, 633, 634, 637, 642, 643, 644, 645, 650, 652, 653, 654, 656, 657, 658, 659, 660, 663, 664, 665, 666, 667, 668, 670, 671, 673, 676, 677, 678, 680, 682, 684, 685, 686, 688, 689, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 717, 719, 720, 721, 723, 724, 725, 726, 727, 728, 729, 730, 736, 738, 739, 741, 742, 743, 744, 745, 748, 749, 750
3	500	1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 20, 21, 22, 23, 25, 27, 28, 29, 30, 31, 35, 37, 38, 39, 40, 42, 43, 44, 45, 46, 50, 51, 52, 53, 55, 56, 58, 62, 63, 66, 69, 70, 71, 72, 73, 74, 75, 76, 77, 79, 80, 81, 83, 84, 85, 87, 88, 89, 90, 91, 92, 93, 94, 95, 98, 99, 102, 105, 106, 108, 110, 111, 112, 113, 117, 118, 119, 120, 121, 122, 124, 125, 126, 127, 128, 129, 130, 133, 134, 136, 138, 139, 141, 142, 143, 146, 148, 149, 150, 151, 152, 154, 156, 159, 160, 161, 163, 164, 166, 167, 168, 169, 170, 171, 173, 174, 175, 176, 179, 180, 182, 183, 184, 187, 190, 191, 192, 194, 195, 196, 197, 198, 203, 204, 205, 206, 207, 208, 209, 211, 213, 214, 215, 217, 219, 221, 222, 223, 225, 227, 228, 229, 230, 231, 234, 235, 236, 239, 241, 244, 245, 246, 248, 250, 252, 253, 255, 256, 257, 258, 260, 261, 263, 265, 268, 271, 272, 273, 274, 275, 276, 277, 282, 288, 290, 291, 292, 294, 295, 296, 300, 301, 303, 304, 305, 307, 308, 310, 311, 312, 314, 315, 316, 317, 319, 320, 322, 323, 324, 325, 327, 328, 331, 332, 334, 335, 336, 337, 340, 341, 344, 345, 346, 347, 348, 349, 353, 354, 356, 362, 363, 366, 367, 369, 370, 371, 372, 374, 375, 376, 379, 382, 383, 384, 385, 386, 389, 390, 391, 392, 394, 396, 397, 398, 399, 401, 402, 403, 404, 405, 406, 407, 410, 412, 414, 418, 419, 420, 422, 423, 424, 426, 427, 428, 430, 431, 432, 433, 435, 437, 440, 441, 443, 446, 447, 448, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 463, 468, 469, 470, 471, 473, 474, 475, 476, 477, 478, 479, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 496, 499, 500, 501, 502, 503, 504, 505, 508, 509, 512, 513, 515, 516, 517, 519, 520, 521, 522, 523, 524, 525, 526, 530, 531, 532, 533, 536, 537, 538, 539, 542, 543, 544, 545, 547, 549, 550, 551, 552, 555, 557, 560, 561, 564, 565, 567, 568, 569, 570, 571, 573, 575, 576, 577, 578, 581, 582, 583, 586, 587, 588, 590, 592, 593, 594, 595, 596, 597, 601, 602, 603, 604, 606, 607, 609, 610, 612, 614, 615, 616, 617, 618, 619, 622, 623, 625, 627, 628, 629, 630, 631, 632, 633, 634, 637, 642, 643, 644, 645, 650, 652, 653, 654, 656, 657, 658, 659, 660, 663, 664, 665, 666, 667, 668, 670, 671, 673, 676, 677, 678, 680, 682, 684, 685, 686, 688, 689, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 717, 719, 720, 721, 723, 724, 725, 726, 727, 728, 729, 730, 736, 738, 739, 741, 742, 743, 744, 745, 748, 749, 750
4	500	1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 20, 21, 22, 23, 25, 27, 28, 29, 30, 31, 35, 37, 38, 39, 40, 42, 43, 44, 45, 46, 50, 51, 52, 53, 55, 56, 58, 62, 63, 66, 69, 70, 71, 72, 73, 74, 75, 76, 77, 79, 80, 81, 83, 84, 85, 87, 88, 89, 90, 91, 92, 93, 94, 95, 98, 99, 102, 105, 106, 108, 110, 111, 112, 113, 117, 118, 119, 120, 121, 122, 124, 125, 126, 127, 128, 129, 130, 133, 134, 136, 138, 139, 141, 142, 143, 146, 148, 149, 150, 151, 152, 154, 156, 159, 160, 161, 163, 164, 166, 167, 168, 169, 170, 171, 173, 174, 175, 176, 179, 180, 182, 183, 184, 187, 190, 191, 192, 194, 195, 196, 197, 198, 203, 204, 205, 206, 207, 208, 209, 211, 213, 214, 215, 217, 219, 221, 222, 223, 225, 227, 228, 229, 230, 231, 234, 235, 236, 239, 241, 244, 245, 246, 248, 250, 252, 253, 255, 256, 257, 258, 260, 261, 263, 265, 268, 271, 272, 273, 274, 275, 276, 277, 282, 288, 290, 291, 292, 294, 295, 296, 300, 301, 303, 304, 305, 307, 308, 310, 311, 312, 314, 315, 316, 317, 319, 320, 322, 323, 324, 325, 327, 328, 331, 332, 334, 335, 336, 337, 340, 341, 344, 345, 346, 347, 348, 349, 353, 354, 356, 362, 363, 366, 367, 369, 370, 371, 372, 374, 375, 376, 379, 382, 383, 384, 385, 386, 389, 390, 391, 392, 394, 396, 397, 398, 399, 401, 402, 403, 404, 405, 406, 407, 410, 412, 414, 418, 419, 420, 422, 423, 424, 426, 427, 428, 430, 431, 432, 433, 435, 437, 440, 441, 443, 446, 447, 448, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 463, 468, 469, 470, 471, 473, 474, 475, 476, 477, 478, 479, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 496, 499, 500, 501, 502, 503, 504, 505, 508, 509, 512, 513, 515, 516, 517, 519, 520, 521, 522, 523, 524, 525, 526, 530, 531, 532, 533, 536, 537, 538, 539, 542, 543, 544, 545, 547, 549, 550, 551, 552, 555, 557, 560, 561, 564, 565, 567, 568, 569, 570, 571, 573, 575, 576, 577, 578, 581, 582, 583, 586, 587, 588, 590, 592, 593, 594, 595, 596, 597, 601, 602, 603, 604, 606, 607, 609, 610, 612, 614, 615, 616, 617, 618, 619, 622, 623, 625, 627, 628, 629, 630, 631, 632, 633, 634, 637, 642, 643, 644, 645, 650, 652, 653, 654, 656, 657, 658, 659, 660, 663, 664, 665, 666, 667, 668, 670, 671, 673, 676, 677, 678, 680, 682, 684, 685, 686, 688, 689, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 717, 719, 720, 721, 723, 724, 725, 726, 727, 728, 729, 730, 736, 738, 739, 741, 742, 743, 744, 745, 748, 749, 750

no_uam_projection_historical_state$summary |>
  transmute(
    `Draws expected` = draws_expected,
    `Draws checked` = draws_evaluated,
    `Invalid draws` = invalid_draws,
    `Non-finite draws` = non_finite_draws,
    `Invalid state cells` = invalid_state_cells,
    `Maximum raw harvest` = max_raw_harvest_rate,
    `Minimum abundance` = min_number,
    `Maximum continuation penalty` = max_harvest_penalty
  ) |>
  mutate(across(
    c(
      `Draws expected`,
      `Draws checked`,
      `Invalid draws`,
      `Non-finite draws`,
      `Invalid state cells`
    ),
    ~ formatC(as.integer(.x), format = "d", big.mark = ",")
  )) |>
  mutate(across(where(is.numeric), ~ format_decimal(.x, digits = 4))) |>
  kable()

Table 4: Independent biological-state audit of the selected NoUAM posterior draws.

Draws expected	Draws checked	Invalid draws	Non-finite draws	Invalid state cells	Maximum raw harvest	Minimum abundance	Maximum continuation penalty
2,000	2,000	0	0	0	0.8730	378.6193	0.0000

Projection Inputs

Projected recruitment deviates are generated using automatically selected ARIMA models fitted to the posterior recruitment-deviate time series. Fitted selectivity is retained through 2025; from 2026 onward it is fixed at the mean selectivity-at-age over 2016-2025. Catch uses fitted observations through the final data year, a four-year pre-CTP nominal TAC block, and then CTP TACs split using the input six-fishery proportions, with no projected TAC assigned to LL3 or LL4. Nominal future TAC allocations and biological removals are kept separate, with the configured UAM multipliers applied only to future removals.

format_year_range <- function(years) {
  if (!length(years)) "none" else paste0(min(years), "-", max(years))
}
format_count <- function(x) {
  formatC(x, format = "f", digits = 0, big.mark = ",")
}
format_integer <- function(x) {
  formatC(x, format = "f", digits = 0)
}
format_elapsed_minutes <- function(seconds, digits = 1) {
  paste0(formatC(seconds / 60, format = "f", digits = digits, big.mark = ","), " minutes")
}
format_count_range <- function(x, missing = "&ndash;") {
  x <- x[!is.na(x)]
  if (!length(x)) {
    missing
  } else if (length(unique(x)) == 1L) {
    format_count(x[1])
  } else {
    paste0(format_count(min(x)), "&ndash;", format_count(max(x)))
  }
}
projected_point_color <- "#D55E00"
format_projected_cell <- function(x, source) {
  out <- as.character(x)
  ifelse(
    source == "Projected",
    paste0("<span style=\"color:", projected_point_color, "; font-weight:600;\">", out, "</span>"),
    out
  )
}
format_projected_count <- function(x, source) {
  format_projected_cell(format_count(x), source)
}
format_projected_integer <- function(x, source) {
  format_projected_cell(format_integer(x), source)
}

Recruitment deviates

projection_rec_samp_years <- base_fit$data$first_yr:(projection_first_yr - 1L)
recruitment_pin_year <- projection_first_yr
projection_recruitment_option <- "auto"
projection_recruitment_bootstrap <- TRUE
projection_recdev_signature <- paste(
  projection_recdev_source_signature,
  projection_seed,
  paste(projection_iters, collapse = ","),
  projection_first_yr,
  projection_last_yr,
  paste(projection_rec_samp_years, collapse = ","),
  recruitment_pin_year,
  projection_recruitment_option,
  projection_recruitment_bootstrap,
  "project_rec_devs_arima_v3",
  sep = "|"
)

projection_recdev_cache <- if (file.exists(projection_recdev_file) && !rebuild_projection_inputs) {
  read_rds_cache(projection_recdev_file)
} else {
  NULL
}

if (
  cache_record_is_compatible(
    projection_recdev_cache,
    projection_recdev_signature,
    "proj_rdev_ARIMA",
    legacy_projection_input_fingerprint("base_recdev")
  ) && is.list(projection_recdev_cache$proj_rdev_ARIMA) &&
    all(c("rdev_y", "proj_rdev_y", "arima_pars") %in%
      names(projection_recdev_cache$proj_rdev_ARIMA)) &&
    is.matrix(projection_recdev_cache$proj_rdev_ARIMA$rdev_y) &&
    is.matrix(projection_recdev_cache$proj_rdev_ARIMA$proj_rdev_y) &&
    is.matrix(projection_recdev_cache$proj_rdev_ARIMA$arima_pars)
) {
  proj_rdev_ARIMA <- projection_recdev_cache$proj_rdev_ARIMA
} else {
  require_projection_computation("projected recruitment-deviate")
  set.seed(projection_seed)
  proj_rdev_ARIMA <- project_rec_devs(
    data = base_fit$data,
    obj = obj_projection,
    mcmc = projection_mcmc,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = projection_rec_samp_years,
    iters = projection_iters,
    pin_year = recruitment_pin_year,
    option = projection_recruitment_option,
    bootstrap = projection_recruitment_bootstrap
  )
  atomic_save_rds(
    list(signature = projection_recdev_signature, proj_rdev_ARIMA = proj_rdev_ARIMA),
    projection_recdev_file
  )
}

projection_rdev_dynamics_y <- proj_rdev_ARIMA$proj_rdev_y[, as.character(projection_years), drop = FALSE]

The recruitment-deviation process standard deviation used when fitting the assessment is fixed at \(\sigma_R = 0.600\). This parameter sets the scale of the recruitment-deviation prior and the lognormal bias correction; it does not require the realised annual deviations in any posterior draw to have a sample standard deviation of exactly 0.6.

The selected projection method is ARIMA. For each posterior draw, sbt::project_rec_devs() uses BIC to choose the autoregressive, differencing, and moving-average orders, with each component limited to a maximum order of 3, then simulates innovations by resampling fitted residuals. The first projected value is pinned to the matching posterior draw for 2022.

fitted_sigma_r <- exp(as.numeric(base_fit$parameters$par_log_sigma_r))
if (length(fitted_sigma_r) != 1L || !is.finite(fitted_sigma_r)) {
  stop("The fitted recruitment-deviation SD must be one finite value.", call. = FALSE)
}

rdev_years <- as.integer(colnames(proj_rdev_ARIMA$rdev_y))
arima_training_rdev <- proj_rdev_ARIMA$rdev_y[
  , as.character(projection_rec_samp_years), drop = FALSE
]
recruitment_sd_draws <- tibble(
  iteration = seq_len(nrow(proj_rdev_ARIMA$rdev_y)),
  historical_sd = apply(proj_rdev_ARIMA$rdev_y, 1L, stats::sd),
  training_sd = apply(arima_training_rdev, 1L, stats::sd),
  projected_series_sd = apply(proj_rdev_ARIMA$proj_rdev_y, 1L, stats::sd)
)

summarise_recruitment_sd <- function(x) {
  tibble(
    lower = quantile(x, 0.025, na.rm = TRUE),
    median = median(x, na.rm = TRUE),
    upper = quantile(x, 0.975, na.rm = TRUE)
  )
}

recruitment_sd_table <- bind_rows(
  tibble(
    Quantity = "Fitting value: process SD (sigma_R)",
    lower = fitted_sigma_r,
    median = fitted_sigma_r,
    upper = fitted_sigma_r
  ),
  summarise_recruitment_sd(recruitment_sd_draws$historical_sd) |>
    mutate(Quantity = paste0(
      "Realised fitted deviations (", min(rdev_years), "–", max(rdev_years), ")"
    )),
  summarise_recruitment_sd(recruitment_sd_draws$training_sd) |>
    mutate(Quantity = paste0(
      "ARIMA training series (", min(projection_rec_samp_years), "–",
      max(projection_rec_samp_years), ")"
    )),
  summarise_recruitment_sd(recruitment_sd_draws$projected_series_sd) |>
    mutate(Quantity = paste0(
      "Realised projected series (", min(projection_years), "–",
      max(projection_years), ")"
    ))
) |>
  select(Quantity, lower, median, upper)

recruitment_sd_table |>
  mutate(across(c(lower, median, upper), ~ format_decimal(.x, 3))) |>
  rename(
    `2.5%` = lower,
    Median = median,
    `97.5%` = upper
  ) |>
  knitr::kable(align = c("l", "r", "r", "r"))

Table 5: Recruitment-deviation standard deviations. Posterior quantities are summarised across the 2,000 balanced-grid draws using the median and equal-tailed 95% interval. The fitting value is fixed, so its three entries are identical.

Quantity	2.5%	Median	97.5%
Fitting value: process SD (sigma_R)	0.600	0.600	0.600
Realised fitted deviations (1931–2025)	0.399	0.461	0.528
ARIMA training series (1931–2021)	0.395	0.454	0.519
Realised projected series (2022–2035)	0.227	0.403	0.664

The realised SD of a 14-year projected ARIMA series reflects serial dependence, the short projection horizon, residual resampling, and the value pinned in 2022. It is therefore not an estimate of the underlying innovation SD. The accepted projection cache does not retain the fitted ARIMA residuals or their SD, so that quantity cannot be recovered exactly without recreating the original time-series fits under their original software environment.

recruitment_sd_draws |>
  select(
    `Historical deviations (1931–2025)` = historical_sd,
    `ARIMA training series (1931–2021)` = training_sd,
    `Projected series (2022–2035)` = projected_series_sd
  ) |>
  pivot_longer(everything(), names_to = "Quantity", values_to = "SD") |>
  ggplot(aes(x = SD, fill = Quantity, color = Quantity)) +
  geom_density(alpha = 0.16, linewidth = 0.7) +
  geom_vline(
    xintercept = fitted_sigma_r,
    linetype = "dashed",
    color = "black",
    linewidth = 0.6
  ) +
  scale_x_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.04))) +
  labs(x = "Standard deviation", y = "Density", fill = NULL, color = NULL) +
  theme(legend.position = "bottom")

Figure 1: Distributions across balanced-grid posterior draws of the realised historical recruitment-deviation SD, the SD within the ARIMA training window, and the realised SD of the projected 2022–2035 series. The dashed line is the fixed sigma_R = 0.6 used during model fitting.

The accepted cache reliably retains the selected AR and MA orders, but not the non-seasonal differencing order. The order summary therefore reports only the two components that can be audited directly from the accepted artefact.

as_tibble(as.data.frame.table(
  proj_rdev_ARIMA$arima_pars[, c("AR", "MA"), drop = FALSE],
  responseName = "value"
)) |>
  mutate(value = factor(value)) |>
  count(statistic, value, name = "draws") |>
  group_by(statistic) |>
  mutate(share = draws / sum(draws)) |>
  ungroup() |>
  ggplot(aes(x = value, y = share)) +
  geom_col(fill = "#ECA82C", alpha = 0.85, width = 0.72) +
  geom_text(aes(label = percent(share, accuracy = 1)), vjust = -0.35, size = 3) +
  facet_wrap(vars(statistic), nrow = 1) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 1), expand = expansion(mult = c(0, 0.08))) +
  labs(x = "Selected order", y = "Projection draws")

Figure 2: Autoregressive (AR) and moving-average (MA) orders retained in the accepted recruitment-deviation projection cache. Bars show the share of projection draws selecting each order.

hist_rdev_df <- as_tibble(
  as.data.frame.table(proj_rdev_ARIMA$rdev_y, responseName = "value")
) |>
  transmute(
    iteration = as.integer(iteration),
    year = as.integer(as.character(year)),
    value = value,
    series = "Historical"
  )

proj_rdev_df <- as_tibble(
  as.data.frame.table(proj_rdev_ARIMA$proj_rdev_y, responseName = "value")
) |>
  transmute(
    iteration = as.integer(iteration),
    year = as.integer(as.character(year)),
    value = value,
    series = "Projected ARIMA"
  )

set.seed(42)
rdev_worms <- sample(unique(hist_rdev_df$iteration), min(8L, n_distinct(hist_rdev_df$iteration)))

bind_rows(hist_rdev_df, proj_rdev_df) |>
  ggplot(aes(x = year, y = value)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_vline(xintercept = recruitment_pin_year, linetype = "dotted", color = "grey35") +
  geom_line(
    data = ~ filter(.x, iteration %in% rdev_worms),
    aes(group = interaction(series, iteration), color = series),
    alpha = 0.18,
    linewidth = 0.35
  ) +
  stat_summary(
    aes(group = series, color = series),
    geom = "line",
    fun = median,
    linewidth = 0.75
  ) +
  stat_summary(
    aes(group = series, fill = series),
    geom = "ribbon",
    fun.data = ~ data.frame(
      ymin = quantile(.x, 0.025, na.rm = TRUE),
      ymax = quantile(.x, 0.975, na.rm = TRUE)
    ),
    alpha = 0.22,
    color = NA
  ) +
  stat_summary(
    data = filter(proj_rdev_df, year == recruitment_pin_year),
    geom = "point",
    fun = median,
    color = "black",
    size = 2
  ) +
  scale_color_manual(values = c("Historical" = "#00BFC4", "Projected ARIMA" = "#F8766D")) +
  scale_fill_manual(values = c("Historical" = "#00BFC4", "Projected ARIMA" = "#F8766D")) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Recruitment deviate", color = NULL, fill = NULL)

Figure 3: Historical and projected recruitment deviates. Historical values are from the balanced grid posterior; projected values are generated from ARIMA fits beginning in 2022 and pinned to each posterior draw’s 2022 recruitment deviate.

Selectivity

projection_selectivity_retained_through <- 2025L
projection_sel_samp_years <- 2016L:2025L
if (!identical(as.integer(base_fit$data$last_yr), 2025L) ||
    !all(projection_sel_samp_years %in%
      base_fit$data$first_yr:base_fit$data$last_yr)) {
  stop(
    "The projection contract requires a terminal 2025 assessment and selectivity sample years 2016-2025.",
    call. = FALSE
  )
}
projection_selectivity_signature <- paste(
  projection_selectivity_source_signature,
  paste(projection_iters, collapse = ","),
  projection_first_yr,
  projection_last_yr,
  paste(projection_sel_samp_years, collapse = ","),
  projection_selectivity_retained_through,
  "project_selectivity_retain_fitted_then_mean_last_v2",
  sep = "|"
)

projection_selectivity_cache <- if (file.exists(projection_selectivity_file) && !rebuild_projection_inputs) {
  read_rds_cache(projection_selectivity_file)
} else {
  NULL
}

if (
  cache_record_is_compatible(
    projection_selectivity_cache,
    projection_selectivity_signature,
    "proj_sel_fya",
    legacy_projection_input_fingerprint("base_selectivity")
  ) && is.array(projection_selectivity_cache$proj_sel_fya)
) {
  proj_sel_fya <- projection_selectivity_cache$proj_sel_fya
} else {
  require_projection_computation("projected selectivity")
  proj_sel_fya <- project_selectivity(
    data = base_fit$data,
    obj = obj_projection,
    mcmc = projection_mcmc,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = projection_sel_samp_years,
    iters = projection_iters,
    option = "mean_last",
    n_years = length(projection_sel_samp_years),
    retain_fitted_through = projection_selectivity_retained_through
  )
  atomic_save_rds(
    list(signature = projection_selectivity_signature, proj_sel_fya = proj_sel_fya),
    projection_selectivity_file
  )
}

projection_sel_dynamics_fya <- proj_sel_fya[, seq_len(base_fit$data$n_fishery), , , drop = FALSE]

no_uam_projection_rec_samp_years <-
  no_uam_fit$data$first_yr:(projection_first_yr - 1L)
no_uam_projection_sel_samp_years <- projection_sel_samp_years
if (!identical(as.integer(no_uam_fit$data$last_yr), 2025L) ||
    !all(no_uam_projection_sel_samp_years %in%
      no_uam_fit$data$first_yr:no_uam_fit$data$last_yr)) {
  stop(
    "The NoUAM projection contract requires fitted selectivity through 2025 and sample years 2016-2025.",
    call. = FALSE
  )
}
no_uam_projection_recdev_signature <- esc31_object_md5(list(
  source_signature = no_uam_projection_recdev_source_signature,
  seed = projection_seed,
  iters = no_uam_projection_iters,
  years = projection_years,
  sample_years = no_uam_projection_rec_samp_years,
  pin_year = recruitment_pin_year,
  option = projection_recruitment_option,
  bootstrap = projection_recruitment_bootstrap,
  implementation = "project_rec_devs_arima_v3"
))
no_uam_projection_recdev_cache <- if (
  file.exists(no_uam_projection_recdev_file) && !rebuild_projection_inputs
) {
  read_rds_cache(no_uam_projection_recdev_file)
} else {
  NULL
}
if (
  cache_record_is_compatible(
    no_uam_projection_recdev_cache,
    no_uam_projection_recdev_signature,
    "proj_rdev_ARIMA",
    legacy_projection_input_fingerprint("no_uam_recdev")
  ) &&
    is.list(no_uam_projection_recdev_cache$proj_rdev_ARIMA) &&
    all(c("rdev_y", "proj_rdev_y", "arima_pars") %in%
      names(no_uam_projection_recdev_cache$proj_rdev_ARIMA))
) {
  no_uam_proj_rdev_ARIMA <-
    no_uam_projection_recdev_cache$proj_rdev_ARIMA
} else {
  require_no_uam_projection_computation(
    "projected recruitment-deviate"
  )
  set.seed(projection_seed)
  no_uam_proj_rdev_ARIMA <- project_rec_devs(
    data = no_uam_fit$data,
    obj = no_uam_obj_projection,
    mcmc = no_uam_mcmc,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = no_uam_projection_rec_samp_years,
    iters = no_uam_projection_iters,
    pin_year = recruitment_pin_year,
    option = projection_recruitment_option,
    bootstrap = projection_recruitment_bootstrap
  )
  atomic_save_rds(
    list(
      signature = no_uam_projection_recdev_signature,
      proj_rdev_ARIMA = no_uam_proj_rdev_ARIMA
    ),
    no_uam_projection_recdev_file
  )
}
no_uam_projection_rdev_dynamics_y <-
  no_uam_proj_rdev_ARIMA$proj_rdev_y[
    , as.character(projection_years), drop = FALSE
  ]

no_uam_projection_selectivity_signature <- esc31_object_md5(list(
  source_signature = no_uam_projection_selectivity_source_signature,
  iters = no_uam_projection_iters,
  years = projection_years,
  sample_years = no_uam_projection_sel_samp_years,
  retain_fitted_through = projection_selectivity_retained_through,
  implementation = "project_selectivity_retain_fitted_then_mean_last_v2"
))
no_uam_projection_selectivity_cache <- if (
  file.exists(no_uam_projection_selectivity_file) &&
    !rebuild_projection_inputs
) {
  read_rds_cache(no_uam_projection_selectivity_file)
} else {
  NULL
}
if (
  cache_record_is_compatible(
    no_uam_projection_selectivity_cache,
    no_uam_projection_selectivity_signature,
    "proj_sel_fya",
    legacy_projection_input_fingerprint("no_uam_selectivity")
  ) &&
    is.array(no_uam_projection_selectivity_cache$proj_sel_fya)
) {
  no_uam_proj_sel_fya <-
    no_uam_projection_selectivity_cache$proj_sel_fya
} else {
  require_no_uam_projection_computation("projected selectivity")
  no_uam_proj_sel_fya <- project_selectivity(
    data = no_uam_fit$data,
    obj = no_uam_obj_projection,
    mcmc = no_uam_mcmc,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = no_uam_projection_sel_samp_years,
    iters = no_uam_projection_iters,
    option = "mean_last",
    n_years = length(no_uam_projection_sel_samp_years),
    retain_fitted_through = projection_selectivity_retained_through
  )
  atomic_save_rds(
    list(
      signature = no_uam_projection_selectivity_signature,
      proj_sel_fya = no_uam_proj_sel_fya
    ),
    no_uam_projection_selectivity_file
  )
}

mle_grid_projection_rec_samp_years <-
  mle_grid_projection_data$first_yr:(projection_first_yr - 1L)
mle_grid_projection_sel_samp_years <- projection_sel_samp_years
if (!identical(as.integer(mle_grid_projection_data$last_yr), 2025L) ||
    !all(mle_grid_projection_sel_samp_years %in%
      mle_grid_projection_data$first_yr:
        mle_grid_projection_data$last_yr)) {
  stop(
    "The sampled MLE-grid projection requires fitted selectivity through ",
    "2025 and sample years 2016-2025.",
    call. = FALSE
  )
}

mle_grid_projection_recdev_signature <- esc31_object_md5(list(
  source_signature = mle_grid_projection_recdev_source_signature,
  seed = projection_seed,
  iters = mle_grid_projection_iters,
  years = projection_years,
  sample_years = mle_grid_projection_rec_samp_years,
  pin_year = recruitment_pin_year,
  option = projection_recruitment_option,
  bootstrap = projection_recruitment_bootstrap,
  implementation = "project_rec_devs_arima_v3"
))
mle_grid_projection_recdev_cache <- if (
  file.exists(mle_grid_projection_recdev_file) &&
    !rebuild_projection_inputs
) {
  read_rds_cache(mle_grid_projection_recdev_file)
} else {
  NULL
}
if (
  cache_record_is_compatible(
    mle_grid_projection_recdev_cache,
    mle_grid_projection_recdev_signature,
    "proj_rdev_ARIMA",
    legacy_projection_input_fingerprint("mle_grid_recdev")
  ) &&
    is.list(mle_grid_projection_recdev_cache$proj_rdev_ARIMA) &&
    all(c("rdev_y", "proj_rdev_y", "arima_pars") %in%
      names(mle_grid_projection_recdev_cache$proj_rdev_ARIMA))
) {
  mle_grid_proj_rdev_ARIMA <-
    mle_grid_projection_recdev_cache$proj_rdev_ARIMA
} else {
  require_mle_grid_projection_computation(
    "projected recruitment-deviate"
  )
  set.seed(projection_seed)
  mle_grid_proj_rdev_ARIMA <- project_rec_devs(
    data = mle_grid_projection_data,
    obj = mle_grid_projection_object,
    mcmc = mle_grid_projection_fit,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = mle_grid_projection_rec_samp_years,
    iters = mle_grid_projection_iters,
    pin_year = recruitment_pin_year,
    option = projection_recruitment_option,
    bootstrap = projection_recruitment_bootstrap
  )
  atomic_save_rds(
    list(
      signature = mle_grid_projection_recdev_signature,
      proj_rdev_ARIMA = mle_grid_proj_rdev_ARIMA
    ),
    mle_grid_projection_recdev_file
  )
}
mle_grid_projection_rdev_dynamics_y <-
  mle_grid_proj_rdev_ARIMA$proj_rdev_y[
    , as.character(projection_years), drop = FALSE
  ]

mle_grid_projection_selectivity_signature <- esc31_object_md5(list(
  source_signature = mle_grid_projection_selectivity_source_signature,
  iters = mle_grid_projection_iters,
  years = projection_years,
  sample_years = mle_grid_projection_sel_samp_years,
  retain_fitted_through = projection_selectivity_retained_through,
  implementation = "project_selectivity_retain_fitted_then_mean_last_v2"
))
mle_grid_projection_selectivity_cache <- if (
  file.exists(mle_grid_projection_selectivity_file) &&
    !rebuild_projection_inputs
) {
  read_rds_cache(mle_grid_projection_selectivity_file)
} else {
  NULL
}
if (
  cache_record_is_compatible(
    mle_grid_projection_selectivity_cache,
    mle_grid_projection_selectivity_signature,
    "proj_sel_fya",
    legacy_projection_input_fingerprint("mle_grid_selectivity")
  ) &&
    is.array(mle_grid_projection_selectivity_cache$proj_sel_fya)
) {
  mle_grid_proj_sel_fya <-
    mle_grid_projection_selectivity_cache$proj_sel_fya
} else {
  require_mle_grid_projection_computation("projected selectivity")
  mle_grid_proj_sel_fya <- project_selectivity(
    data = mle_grid_projection_data,
    obj = mle_grid_projection_object,
    mcmc = mle_grid_projection_fit,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    samp_years = mle_grid_projection_sel_samp_years,
    iters = mle_grid_projection_iters,
    option = "mean_last",
    n_years = length(mle_grid_projection_sel_samp_years),
    retain_fitted_through = projection_selectivity_retained_through
  )
  atomic_save_rds(
    list(
      signature = mle_grid_projection_selectivity_signature,
      proj_sel_fya = mle_grid_proj_sel_fya
    ),
    mle_grid_projection_selectivity_file
  )
}

fishery_names <- c("LL1", "LL2", "LL3", "LL4", "Indonesia", "Australia", "CPUE")

proj_sel_df <- as_tibble(
  as.data.frame.table(proj_sel_fya, responseName = "value")
) |>
  transmute(
    iteration = as.integer(iteration),
    fishery = factor(fishery_names[as.integer(fishery)], levels = fishery_names),
    year = as.integer(as.character(year)),
    age = as.numeric(as.character(age)),
    value = value,
    period = if_else(
      year <= projection_selectivity_retained_through,
      "Fitted through 2025",
      "2016-2025 mean from 2026"
    )
  )

proj_sel_df |>
  group_by(fishery, year, age, period) |>
  summarise(value = median(value, na.rm = TRUE), .groups = "drop") |>
  ggplot(aes(
    x = age,
    y = value,
    group = year,
    color = factor(year),
    linetype = period
  )) +
  geom_line(linewidth = 0.5, alpha = 0.85) +
  facet_wrap(vars(fishery), scales = "free_y") +
  scale_y_zero() +
  scale_linetype_manual(values = c(
    "Fitted through 2025" = "solid",
    "2016-2025 mean from 2026" = "dashed"
  )) +
  labs(
    x = "Age",
    y = "Selectivity",
    color = "Projection year",
    linetype = NULL
  )

Figure 4: Median selectivity at age by fishery and projection year. Posterior fitted selectivity is retained for 2022-2025; from 2026 onward each draw uses its 2016-2025 mean selectivity.

CPUE q

projection_cpue_q_signature <- paste(
  projection_cpue_q_source_signature,
  paste(projection_iters, collapse = ","),
  projection_first_yr,
  projection_last_yr,
  "cpue_q_v3",
  sep = "|"
)

if (file.exists(projection_cpue_q_file) && !rebuild_projection_inputs) {
  projection_cpue_q_cache <- read_rds_cache(projection_cpue_q_file)
} else {
  projection_cpue_q_cache <- NULL
}

if (
  cache_record_is_compatible(
    projection_cpue_q_cache,
    projection_cpue_q_signature,
    "draws",
    legacy_projection_input_fingerprint("base_cpue_q")
  ) && is.data.frame(projection_cpue_q_cache$draws) &&
    all(c("draw", "year", "series", "q") %in%
      names(projection_cpue_q_cache$draws))
) {
  cpue_q_draws <- projection_cpue_q_cache$draws
} else {
  require_projection_computation("projected CPUE-q")
  cpue_adjust_years <- base_fit$data$cpue_years + base_fit$data$first_yr - 1L
  cpue_q_years <- min(cpue_adjust_years):projection_last_yr
  fixed_pars <- obj_projection$env$parList(obj_projection$env$last.par.best)
  cpue_log_q_draws <- if ("par_log_cpue_q" %in% names(projection_post)) {
    projection_post$par_log_cpue_q
  } else {
    rep(tail(as.numeric(fixed_pars$par_log_cpue_q), 1L), nrow(projection_post))
  }
  cpue_creep_draws <- if ("par_cpue_creep" %in% names(projection_post)) {
    projection_post$par_cpue_creep
  } else {
    rep(as.numeric(fixed_pars$par_cpue_creep), nrow(projection_post))
  }

  cpue_q_draws <- map_dfr(seq_len(nrow(projection_post)), function(i) {
    cpue_creep_i <- as.numeric(cpue_creep_draws[i])

    historical_adjust <- 1 + cpue_creep_i * (seq_along(cpue_adjust_years) - 1L)
    adjustment <- rep(NA_real_, length(cpue_q_years))
    names(adjustment) <- cpue_q_years
    adjustment[as.character(cpue_adjust_years)] <- historical_adjust

    projection_years_i <- cpue_q_years[cpue_q_years > max(cpue_adjust_years)]
    if (length(projection_years_i)) {
      adjustment[as.character(projection_years_i)] <-
        tail(historical_adjust, 1L) + cpue_creep_i * seq_along(projection_years_i)
    }

    tibble(
      draw = i,
      year = cpue_q_years,
      series = if_else(year <= max(cpue_adjust_years), "Historical", "Projected"),
      q = exp(as.numeric(cpue_log_q_draws[i])) * pmax(adjustment, 1e-8)
    )
  })

  atomic_save_rds(
    list(signature = projection_cpue_q_signature, draws = cpue_q_draws),
    projection_cpue_q_file
  )
}

cpue_q_draws |>
  group_by(series, year) |>
  summarise(projection_quantiles(q), .groups = "drop") |>
  ggplot(aes(x = year, y = median, color = series, fill = series)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_color_manual(values = c("Historical" = "#2F6F73", "Projected" = "#D55E00")) +
  scale_fill_manual(values = c("Historical" = "#72B7B2", "Projected" = "#D55E00")) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Effective CPUE q", color = NULL, fill = NULL)

Figure 5: Historical and projected effective CPUE catchability multiplier (q). Lines show posterior median and shaded ribbons show 95% intervals across the balanced grid posterior; projected values extend the final historical adjustment using the CPUE creep parameter.

CTP schedule

projection_input_signature <- esc31_object_md5(list(
  schema_version = 4L,
  source_signature = projection_source_signature,
  workflow_signature = projection_workflow_identity$scientific_signature,
  recdev_signature = projection_recdev_signature,
  selectivity_signature = projection_selectivity_signature,
  cpue_q_signature = projection_cpue_q_signature,
  projection_years = projection_years,
  projection_iters = projection_iters,
  seed = projection_seed,
  ctp = list(
    tac_schedule = ctp_tac_schedule,
    tac_calculation_lag = ctp_tac_calculation_lag,
    catch_data_lag = ctp_catch_data_lag,
    cpue_data_lag = ctp_cpue_data_lag,
    gt_data_lag = ctp_gt_data_lag,
    ckmr_data_lag = ctp_ckmr_data_lag
  ),
  monitoring = list(
    gt_skip_years = gt_skip_years,
    gt_projection_nrel = gt_projection_nrel,
    gt_projection_nsam = gt_projection_nsam,
    pop_projection_n_juvenile = pop_projection_n_juvenile,
    pop_projection_n_adult = pop_projection_n_adult,
    pop_projection_min_cohort = pop_projection_min_cohort,
    pop_projection_adult_ages = pop_projection_adult_ages,
    hsp_projection_nC = hsp_projection_nC
  ),
  fixed_tac_scenarios = fixed_projection_tac_scenarios,
  scenario_1_inputs = fixed_projection_tac_inputs,
  scenario_2_inputs = fixed_projection_tac_scenario_2_inputs,
  removal_contract = projection_removal_contract,
  projected_recruitment_md5 = esc31_object_md5(projection_rdev_dynamics_y),
  projected_selectivity_md5 = esc31_object_md5(proj_sel_fya),
  harvest_wall_contract = projection_harvest_wall_identity$executable,
  harvest_wall_decision = projection_harvest_wall_identity$decision,
  harvest_wall_application =
    "source_posterior_conditioning_only_not_projection_dynamics",
  implementation =
    "projection_inputs_v8_terminal_selectivity_nominal_tac_and_removals"
))

no_uam_projection_input_signature <- esc31_object_md5(list(
  schema_version = 2L,
  source_signature = no_uam_projection_source_signature,
  recdev_signature = no_uam_projection_recdev_signature,
  selectivity_signature = no_uam_projection_selectivity_signature,
  projection_years = projection_years,
  projection_iters = no_uam_projection_iters,
  projection_seed = projection_seed,
  ctp = list(
    tac_schedule = ctp_tac_schedule,
    tac_calculation_lag = ctp_tac_calculation_lag,
    catch_data_lag = ctp_catch_data_lag,
    cpue_data_lag = ctp_cpue_data_lag,
    gt_data_lag = ctp_gt_data_lag,
    ckmr_data_lag = ctp_ckmr_data_lag
  ),
  monitoring = list(
    gt_skip_years = gt_skip_years,
    gt_projection_nrel = gt_projection_nrel,
    gt_projection_nsam = gt_projection_nsam,
    pop_projection_n_juvenile = pop_projection_n_juvenile,
    pop_projection_n_adult = pop_projection_n_adult,
    pop_projection_min_cohort = pop_projection_min_cohort,
    pop_projection_adult_ages = pop_projection_adult_ages,
    hsp_projection_nC = hsp_projection_nC
  ),
  fixed_tac_scenario = fixed_projection_tac_inputs,
  removal_contract = projection_removal_contract,
  projected_recruitment_md5 =
    esc31_object_md5(no_uam_projection_rdev_dynamics_y),
  projected_selectivity_md5 =
    esc31_object_md5(no_uam_proj_sel_fya),
  comparison =
    "base_scenario_1_with_ncnm_vs_no_uam_without_ncnm_conditioning_and_projection",
  implementation =
    "no_uam_projection_inputs_v2_terminal_selectivity_and_no_future_ncnm"
))

mle_grid_projection_input_signature <- esc31_object_md5(list(
  schema_version = 1L,
  source_signature = mle_grid_projection_source_signature,
  recdev_signature = mle_grid_projection_recdev_signature,
  selectivity_signature = mle_grid_projection_selectivity_signature,
  projection_years = projection_years,
  projection_iters = mle_grid_projection_iters,
  projection_seed = projection_seed,
  ctp = list(
    tac_schedule = ctp_tac_schedule,
    tac_calculation_lag = ctp_tac_calculation_lag,
    catch_data_lag = ctp_catch_data_lag,
    cpue_data_lag = ctp_cpue_data_lag,
    gt_data_lag = ctp_gt_data_lag,
    ckmr_data_lag = ctp_ckmr_data_lag
  ),
  monitoring = list(
    gt_skip_years = gt_skip_years,
    gt_projection_nrel = gt_projection_nrel,
    gt_projection_nsam = gt_projection_nsam,
    pop_projection_n_juvenile = pop_projection_n_juvenile,
    pop_projection_n_adult = pop_projection_n_adult,
    pop_projection_min_cohort = pop_projection_min_cohort,
    pop_projection_adult_ages = pop_projection_adult_ages,
    hsp_projection_nC = hsp_projection_nC
  ),
  fixed_tac_scenario = fixed_projection_tac_inputs,
  removal_contract = projection_removal_contract,
  projected_recruitment_md5 =
    esc31_object_md5(mle_grid_projection_rdev_dynamics_y),
  projected_selectivity_md5 =
    esc31_object_md5(mle_grid_proj_sel_fya),
  comparison =
    "balanced_mcmc_grid_base_vs_full_objective_weighted_direct_m_mle_grid",
  implementation =
    "esc31_report_specific_mle_grid_projection_inputs_v1"
))

projection_actual_catch_years <- projection_years[projection_years <= base_fit$data$last_yr]
projection_future_years <- projection_years[projection_years > base_fit$data$last_yr]
fixed_catch_years <- head(projection_future_years, fixed_catch_n_years)
if (!identical(as.integer(fixed_projection_tac_df$year), as.integer(fixed_catch_years))) {
  stop(
    "Pre-CTP TAC input years must match the first four projection years after the final model year.",
    call. = FALSE
  )
}
mp_start_yr <- if (length(projection_future_years) > length(fixed_catch_years)) {
  projection_future_years[length(fixed_catch_years) + 1L]
} else {
  NA_integer_
}

ctp_schedule <- project_ctp_schedule(
  first_yr = projection_first_yr,
  last_yr = projection_last_yr,
  data_last_yr = base_fit$data$last_yr,
  fixed_catch_n_years = fixed_catch_n_years,
  tac_schedule = ctp_tac_schedule,
  tac_calculation_lag = ctp_tac_calculation_lag,
  catch_data_lag = ctp_catch_data_lag,
  cpue_data_lag = ctp_cpue_data_lag,
  gt_data_lag = ctp_gt_data_lag,
  ckmr_data_lag = ctp_ckmr_data_lag
)

Recruitment deviations begin in 2022 and use the selected ARIMA method; the first value is pinned to the posterior draw for 2022. Posterior fitted selectivity is retained through 2025. Selectivity from 2026 onward is the draw-specific mean over 2016–2025, not the 2012–2021 mean. This deterministic treatment is confirmed for final ESC31 production; no stochastic or shifted-selectivity scenario is part of the production run.

Catch inputs retain the conditioned assessment catches for 2022-2025. The fixed nominal-TAC period is 2026- 2029, with annual totals 22,671, 23,647, 23,647, 23,647, followed by CTP TACs from 2030. Both allocation scenarios use exactly those same fixed total TACs; only their fleet allocation differs in 2027-2029. CTP allocation proportions are calculated from the selected scenario’s non-carryover fixed-TAC years, excluding 2026, and projected LL3 and LL4 nominal allocations are zero. Biological removals are then calculated from those nominal allocations using the factors below. The CTP schedule follows the OMMP16 table of projected TAC changes and lags: TACs are calculated two years before implementation, catch and CPUE are available through TAC year minus 3, gene tagging through TAC year minus 4, and POPs/HSPs through TAC year minus 8 (Commission for the Conservation of Southern Bluefin Tuna 2026).

tibble(
  Fishery = names(base_projection_removal_multiplier_f),
  `Base removal multiplier` =
    as.numeric(base_projection_removal_multiplier_f),
  `NoUAM removal multiplier` =
    as.numeric(no_uam_projection_removal_multiplier_f),
  Interpretation = case_when(
    Fishery == "LL1" ~ "Base includes NCNM; NoUAM excludes it",
    Fishery == "Australia" ~
      "Separate surface-fishery correction retained in both",
    TRUE ~ "No additional projected removal"
  )
) |>
  mutate(across(contains("multiplier"), ~ format_decimal(.x, 2))) |>
  kable()

Table 6: Future nominal-TAC-to-removal multipliers. NoUAM removes the LL1 NCNM addition in both conditioning and projection; the separate Australian surface-fishery correction is retained.

Fishery	Base removal multiplier	NoUAM removal multiplier	Interpretation
LL1	1.11	1.00	Base includes NCNM; NoUAM excludes it
LL2	1.00	1.00	No additional projected removal
LL3	1.00	1.00	No additional projected removal
LL4	1.00	1.00	No additional projected removal
Indonesia	1.00	1.00	No additional projected removal
Australia	1.20	1.20	Separate surface-fishery correction retained in both

Scientific assumptions audit

The assumptions below were checked against OMMP16 and the retained legacy projection controls. Items labelled “workflow decision” are not choices made by OMMP16; they are stated explicitly so they cannot be mistaken for report instructions.

tribble(
  ~Component, ~OMMP16_or_legacy_requirement, ~Implemented_decision,
  "Posterior draws",
  "Not specified by OMMP16",
  "Production uses all 2,000 balanced nine-cell MCMC-grid draws",
  "Recruitment start",
  "Start in Y-4 for Y = 2026",
  "Begin in 2022 and pin 2022 to each fitted posterior draw",
  "Recruitment model",
  "OMMP16 does not choose AR(1), ARIMA order selection, or innovation treatment",
  "User decision: automatic BIC-selected ARIMA with residual bootstrap",
  "Selectivity",
  "Use the most recent 10 years",
  "Retain fitted selectivity through 2025; use the 2016-2025 mean from 2026",
  "Fixed TAC",
  "Four fixed future TAC years and the two Alloc treatments",
  "2026-2029 total TAC is identical in both scenarios; only fleet allocation differs",
  "Post-CTP fleet allocation",
  "Split TAC among fleets using allocation proportions; exact reference years are not selected",
  "Workflow decision: hold each scenario's 2027-2029 mean proportions constant; LL3 and LL4 receive zero projected TAC",
  "Within-year season allocation",
  "Not specified by OMMP16",
  "Workflow decision: use the terminal 2025 conditioned season proportions; these are identical in the base and NoUAM fits",
  "Base NCNM/UAM",
  "Legacy projection control applies fleet removal multipliers",
  "Base future removals use LL1 1.11 and Australia 1.20",
  "NoUAM",
  "Remove NCNM catches from conditioning and projections",
  "Conditioning catch additions are zero; future LL1 multiplier is 1.00",
  "Australian surface correction",
  "Retained legacy NoUAM control is 1, 1, 1, 1.2 over the four legacy fleets",
  "Australia 1.20 is retained in base and NoUAM as distinct from LL1 NCNM",
  "CTP timing",
  "TAC calculation lag 2; catch/CPUE lag 3; GT lag 4; CKMR lag 8",
  "Implemented exactly in the CTP schedule",
  "Monitoring programme",
  "Use the previous projection data-generation programme",
  "GT, HSP, and POP sample sizes, adult ages, and minimum cohort retain the legacy settings",
  "Adaptive GT rule",
  "Remove the former additional 5,000 GT samples when rebuilding is inadequate",
  "No adaptive GT sample increase is implemented",
  "Skipped GT years",
  "Projection code should allow gene-tagging years to be skipped; no years are prescribed",
  "Capability is present; the baseline comparison skips no GT years",
  "Selectivity-shift stress tests",
  "Workplan requests capability for recruitment-failure or smaller-size selectivity scenarios, but gives no numeric scenario",
  "Not selected for ESC31 production; retain fitted selectivity through 2025 and each draw's deterministic 2016-2025 mean from 2026 onward",
  "Projection horizon",
  "No terminal year selected by OMMP16",
  "Confirmed ESC31 decision: project dynamics through 2035 and report state through 2036",
  "Risk display",
  "No 0.2 threshold or short/long reporting windows selected by OMMP16",
  "Confirmed ESC31 decision: relative TRO threshold 0.2; 2030-2032 and 2033-2036 windows"
) |>
  rename(
    Component = Component,
    `Source requirement or status` = OMMP16_or_legacy_requirement,
    `Implemented decision` = Implemented_decision
  ) |>
  kable()

Table 7: Projection scientific assumptions, sources, and implemented decisions.

Component	Source requirement or status	Implemented decision
Posterior draws	Not specified by OMMP16	Production uses all 2,000 balanced nine-cell MCMC-grid draws
Recruitment start	Start in Y-4 for Y = 2026	Begin in 2022 and pin 2022 to each fitted posterior draw
Recruitment model	OMMP16 does not choose AR(1), ARIMA order selection, or innovation treatment	User decision: automatic BIC-selected ARIMA with residual bootstrap
Selectivity	Use the most recent 10 years	Retain fitted selectivity through 2025; use the 2016-2025 mean from 2026
Fixed TAC	Four fixed future TAC years and the two Alloc treatments	2026-2029 total TAC is identical in both scenarios; only fleet allocation differs
Post-CTP fleet allocation	Split TAC among fleets using allocation proportions; exact reference years are not selected	Workflow decision: hold each scenario’s 2027-2029 mean proportions constant; LL3 and LL4 receive zero projected TAC
Within-year season allocation	Not specified by OMMP16	Workflow decision: use the terminal 2025 conditioned season proportions; these are identical in the base and NoUAM fits
Base NCNM/UAM	Legacy projection control applies fleet removal multipliers	Base future removals use LL1 1.11 and Australia 1.20
NoUAM	Remove NCNM catches from conditioning and projections	Conditioning catch additions are zero; future LL1 multiplier is 1.00
Australian surface correction	Retained legacy NoUAM control is 1, 1, 1, 1.2 over the four legacy fleets	Australia 1.20 is retained in base and NoUAM as distinct from LL1 NCNM
CTP timing	TAC calculation lag 2; catch/CPUE lag 3; GT lag 4; CKMR lag 8	Implemented exactly in the CTP schedule
Monitoring programme	Use the previous projection data-generation programme	GT, HSP, and POP sample sizes, adult ages, and minimum cohort retain the legacy settings
Adaptive GT rule	Remove the former additional 5,000 GT samples when rebuilding is inadequate	No adaptive GT sample increase is implemented
Skipped GT years	Projection code should allow gene-tagging years to be skipped; no years are prescribed	Capability is present; the baseline comparison skips no GT years
Selectivity-shift stress tests	Workplan requests capability for recruitment-failure or smaller-size selectivity scenarios, but gives no numeric scenario	Not selected for ESC31 production; retain fitted selectivity through 2025 and each draw’s deterministic 2016-2025 mean from 2026 onward
Projection horizon	No terminal year selected by OMMP16	Confirmed ESC31 decision: project dynamics through 2035 and report state through 2036
Risk display	No 0.2 threshold or short/long reporting windows selected by OMMP16	Confirmed ESC31 decision: relative TRO threshold 0.2; 2030-2032 and 2033-2036 windows

fixed_projection_tac_df |>
  transmute(
    Year = .data$year,
    LL1 = format_count(.data$LL1),
    LL2 = format_count(.data$LL2),
    Indonesia = format_count(.data$Indonesia),
    Australia = format_count(.data$Australia),
    Total = format_count(.data$total)
  ) |>
  kable()

Table 8: Active fixed pre-CTP nominal TAC allocations by year.

Year	LL1	LL2	Indonesia	Australia	Total
2026	14,326	1,534	1,248	5,562	22,671
2027	15,063	1,653	1,370	5,561	23,647
2028	15,063	1,653	1,370	5,561	23,647
2029	15,063	1,653	1,370	5,561	23,647

ctp_schedule |>
  transmute(
    `TAC Year` = .data$year,
    `TAC change?` = .data$tac_change,
    `TAC was calculated` = if_else(is.na(.data$tac_calculation_year), "–", as.character(.data$tac_calculation_year)),
    `Catch data` = if_else(is.na(.data$catch_data_year), "–", as.character(.data$catch_data_year)),
    `CPUE data` = if_else(is.na(.data$cpue_data_year), "–", as.character(.data$cpue_data_year)),
    `Gene tagging (age 2)` = if_else(is.na(.data$gt_data_year), "–", as.character(.data$gt_data_year)),
    `POPs and HSPs` = if_else(is.na(.data$ckmr_data_year), "–", as.character(.data$ckmr_data_year))
  ) |>
  kable(escape = FALSE)

Table 9: OMMP16 schedule of projected TAC changes and lags in data availability.

TAC Year	TAC change?	TAC was calculated	Catch data	CPUE data	Gene tagging (age 2)	POPs and HSPs
2026	hardwired	2022	–	–	–	–
2027	hardwired	2025	2024	2024	2023	2019
2028	hardwired	–	2025	2025	2024	2020
2029	hardwired	–	2026	2026	2025	2021
2030	Yes	2028	2027	2027	2026	2022
2031	No	–	2028	2028	2027	2023
2032	No	–	2029	2029	2028	2024
2033	Yes	2031	2030	2030	2029	2025
2034	No	–	2031	2031	2030	2026
2035	No	–	2032	2032	2031	2027

Projection Dynamics

The figures below separate conditioned history from simulated future output. Black denotes observed or conditioned historical context, purple identifies the fixed nominal-TAC period, and orange denotes future CTP-controlled inputs or outputs. Shaded ribbons are equal-tailed 95% intervals across projection draws. Where a blue dynamics-output line is shown, it is an accounting check that must lie on the configured biological-removal input rather than a separate scenario.

Projection run

load_legacy_projection_reporting_cache <- function(
    inputs, cache_file, arm, fit, post, iters, rdev_y, sel_fya,
    removal_multiplier_f) {
  fingerprint <- legacy_projection_reporting_fingerprint(arm)
  if (is.null(fingerprint) || !file.exists(cache_file) ||
      !identical(
        unname(tools::md5sum(cache_file)),
        fingerprint$file_md5
      )) {
    return(NULL)
  }

  cached <- read_rds_cache(cache_file)
  required_fields <- c(
    "cache_format_version", "run_signature", "projection_config",
    "projection_post", "projection_iters", "projection_years",
    "projection_dynamics_first_yr", "projection_dynamics_years",
    "fixed_catch_years", "mp_start_yr", "ctp_schedule",
    "projection_rdev_dynamics_y", "proj_sel_fya",
    "projection_sel_dynamics_fya", "projection_removal_multiplier_f",
    "projection_n_dyn", "proj_dyn_base", "proj_dyn",
    "ctp_nominal_catch_iysf", "ctp_catch_iysf", "ctp_total_tac_iy",
    "ctp_total_removal_iy", "ctp_status", "ctp_input_snapshots", "timing"
  )
  if (!is.list(cached) ||
      length(setdiff(required_fields, names(cached))) ||
      !identical(cached$cache_format_version, 8L) ||
      !identical(cached$run_signature, fingerprint$run_signature) ||
      !identical(as.integer(cached$projection_n_dyn), projection_n_iter) ||
      !identical(as.integer(cached$projection_years), projection_years) ||
      !identical(as.integer(cached$projection_dynamics_years),
        projection_years) ||
      !identical(as.integer(cached$projection_iters), as.integer(iters)) ||
      !isTRUE(all.equal(
        as.matrix(cached$projection_post),
        as.matrix(post),
        tolerance = 0,
        check.attributes = FALSE
      )) ||
      !identical(cached$projection_rdev_dynamics_y, rdev_y) ||
      !identical(cached$proj_sel_fya, sel_fya) ||
      !isTRUE(all.equal(
        cached$projection_removal_multiplier_f,
        removal_multiplier_f,
        tolerance = 0,
        check.attributes = TRUE
      )) ||
      !identical(
        as.integer(length(cached$fixed_catch_years)),
        as.integer(inputs$fixed_catch_n_years)
      ) ||
      !isTRUE(all.equal(
        cached$projection_config$fixed_projection_tac,
        inputs$fixed_projection_tac,
        tolerance = 0,
        check.attributes = TRUE
      )) ||
      !isTRUE(all.equal(
        cached$projection_config$projection_tac_split_yf,
        inputs$projection_tac_split_yf,
        tolerance = 0,
        check.attributes = TRUE
      )) ||
      !identical(cached$ctp_schedule, ctp_schedule)) {
    return(NULL)
  }

  n_updates <- sum(cached$ctp_schedule$ctp_update %in% TRUE)
  expected_snapshot_names <- as.character(
    cached$ctp_schedule$year[cached$ctp_schedule$ctp_update %in% TRUE]
  )
  catch_totals_match <- isTRUE(all.equal(
    apply(cached$ctp_nominal_catch_iysf, c(1L, 2L), sum),
    cached$ctp_total_tac_iy,
    tolerance = 1e-10,
    check.attributes = FALSE
  )) && isTRUE(all.equal(
    apply(cached$ctp_catch_iysf, c(1L, 2L), sum),
    cached$ctp_total_removal_iy,
    tolerance = 1e-10,
    check.attributes = FALSE
  ))
  if (!is.data.frame(cached$ctp_status) ||
      nrow(cached$ctp_status) != projection_n_iter * n_updates ||
      !all(cached$ctp_status$status == "ok") ||
      !identical(names(cached$ctp_input_snapshots), expected_snapshot_names) ||
      !catch_totals_match) {
    return(NULL)
  }

  state_contract <- biological_state_contract()
  if (length(cached$proj_dyn_base) != projection_n_iter ||
      length(cached$proj_dyn) != projection_n_iter ||
      !all(vapply(cached$proj_dyn_base, is.list, logical(1L))) ||
      !all(vapply(cached$proj_dyn, is.list, logical(1L))) ||
      !is.numeric(fingerprint$maximum_raw_harvest) ||
      length(fingerprint$maximum_raw_harvest) != 1L ||
      !is.finite(fingerprint$maximum_raw_harvest) ||
      fingerprint$maximum_raw_harvest >
        state_contract$hrate_limit + state_contract$hrate_tolerance ||
      !identical(fingerprint$maximum_harvest_penalty, 0)) {
    return(NULL)
  }

  cached$timing$loaded_from_cache <- TRUE
  attr(cached, "esc31_reporting_cache_status") <-
    "audited_format_8_catch_staging_equivalent"
  cached
}

run_projection_scenario <- function(
    inputs, cache_file, cache_signature, arm, post,
    fit = base_fit,
    object = obj_projection,
    mcmc = projection_mcmc,
    iters = projection_iters,
    rdev_y = projection_rdev_dynamics_y,
    sel_fya = proj_sel_fya,
    compute = run_projection_computation,
    removal_multiplier_f = base_projection_removal_multiplier_f,
    legacy_cache_signatures = character()
  ) {
  legacy_reporting_cache <- load_legacy_projection_reporting_cache(
    inputs = inputs,
    cache_file = cache_file,
    arm = arm,
    fit = fit,
    post = post,
    iters = iters,
    rdev_y = rdev_y,
    sel_fya = sel_fya,
    removal_multiplier_f = removal_multiplier_f
  )
  if (!is.null(legacy_reporting_cache) &&
      !rebuild_projection_dynamics && !compute) {
    message(
      "Loaded audited format-8 reporting cache: ", basename(cache_file)
    )
    return(legacy_reporting_cache)
  }

  run_projections(
    data = fit$data,
    object = object,
    mcmc = mcmc,
    iters = iters,
    first_yr = projection_first_yr,
    last_yr = projection_last_yr,
    seed = projection_seed,
    cores = projection_cores,
    fixed_catch_n_years = inputs$fixed_catch_n_years,
    fixed_projection_tac = inputs$fixed_projection_tac,
    ctp_tac_schedule = ctp_tac_schedule,
    ctp_tac_calculation_lag = ctp_tac_calculation_lag,
    ctp_catch_data_lag = ctp_catch_data_lag,
    ctp_cpue_data_lag = ctp_cpue_data_lag,
    ctp_gt_data_lag = ctp_gt_data_lag,
    ctp_ckmr_data_lag = ctp_ckmr_data_lag,
    projection_tac_split_yf = inputs$projection_tac_split_yf,
    projection_removal_multiplier_f = removal_multiplier_f,
    gt_skip_years = gt_skip_years,
    gt_projection_nrel = gt_projection_nrel,
    gt_projection_nsam = gt_projection_nsam,
    pop_projection_n_juvenile = pop_projection_n_juvenile,
    pop_projection_n_adult = pop_projection_n_adult,
    pop_projection_min_cohort = pop_projection_min_cohort,
    pop_projection_adult_ages = pop_projection_adult_ages,
    hsp_projection_nC = hsp_projection_nC,
    rdev_y = rdev_y,
    sel_fya = sel_fya,
    cache_file = cache_file,
    cache_signature = cache_signature,
    legacy_cache_signatures = legacy_cache_signatures,
    cache_only = !compute,
    overwrite = rebuild_projection_dynamics,
    verbose = TRUE
  )
}

projection_run <- run_projection_scenario(
  fixed_projection_tac_inputs,
  cache_file = projection_run_file,
  cache_signature = paste(projection_input_signature, active_fixed_projection_tac_scenario, sep = "|"),
  arm = "base",
  post = projection_post,
  legacy_cache_signatures = legacy_projection_cache_signature("base")
)
projection_run_scenario_2 <- run_projection_scenario(
  fixed_projection_tac_scenario_2_inputs,
  cache_file = projection_run_scenario_2_file,
  cache_signature = paste(projection_input_signature, "scenario_2", sep = "|"),
  arm = "scenario_2",
  post = projection_post,
  legacy_cache_signatures =
    legacy_projection_cache_signature("scenario_2")
)
no_uam_projection_run <- run_projection_scenario(
  fixed_projection_tac_inputs,
  cache_file = no_uam_projection_run_file,
  cache_signature = paste(
    no_uam_projection_input_signature,
    "scenario_1",
    sep = "|"
  ),
  arm = "no_uam",
  post = no_uam_projection_post,
  fit = no_uam_fit,
  object = no_uam_obj_projection,
  mcmc = no_uam_mcmc,
  iters = no_uam_projection_iters,
  rdev_y = no_uam_projection_rdev_dynamics_y,
  sel_fya = no_uam_proj_sel_fya,
  compute = run_no_uam_projection_computation,
  removal_multiplier_f = no_uam_projection_removal_multiplier_f,
  legacy_cache_signatures = legacy_projection_cache_signature("no_uam")
)
mle_grid_projection_run <- run_projection_scenario(
  fixed_projection_tac_inputs,
  cache_file = mle_grid_projection_run_file,
  cache_signature = paste(
    mle_grid_projection_input_signature,
    "scenario_1",
    sep = "|"
  ),
  arm = "mle_grid",
  post = mle_grid_projection_post,
  fit = list(data = mle_grid_projection_data),
  object = mle_grid_projection_object,
  mcmc = mle_grid_projection_fit,
  iters = mle_grid_projection_iters,
  rdev_y = mle_grid_projection_rdev_dynamics_y,
  sel_fya = mle_grid_proj_sel_fya,
  compute = run_mle_grid_projection_computation,
  removal_multiplier_f = base_projection_removal_multiplier_f,
  legacy_cache_signatures =
    legacy_projection_cache_signature("mle_grid")
)
projection_scenario_runs <- setNames(
  list(projection_run, projection_run_scenario_2),
  unname(projection_scenario_labels)
)
if (!isTRUE(all.equal(
      projection_run$projection_removal_multiplier_f,
      base_projection_removal_multiplier_f,
      tolerance = 0
    )) ||
    !isTRUE(all.equal(
      projection_run_scenario_2$projection_removal_multiplier_f,
      base_projection_removal_multiplier_f,
      tolerance = 0
    )) ||
    !isTRUE(all.equal(
      no_uam_projection_run$projection_removal_multiplier_f,
      no_uam_projection_removal_multiplier_f,
      tolerance = 0
    )) ||
    !isTRUE(all.equal(
      mle_grid_projection_run$projection_removal_multiplier_f,
      base_projection_removal_multiplier_f,
      tolerance = 0
    ))) {
  stop("A projection cache does not match its required removal multipliers.",
       call. = FALSE)
}
fixed_tac_indices <- match(fixed_catch_years, projection_years)
fixed_total_difference <- max(abs(
  projection_run$ctp_total_tac_iy[, fixed_tac_indices, drop = FALSE] -
    projection_run_scenario_2$ctp_total_tac_iy[
      , fixed_tac_indices, drop = FALSE
    ]
))
if (!is.finite(fixed_total_difference) ||
    fixed_total_difference > 1e-8) {
  stop(
    "The two projection scenarios do not have identical fixed total TAC.",
    call. = FALSE
  )
}
no_uam_fixed_total_difference <- max(abs(
  projection_run$ctp_total_tac_iy[, fixed_tac_indices, drop = FALSE] -
    no_uam_projection_run$ctp_total_tac_iy[
      , fixed_tac_indices, drop = FALSE
    ]
))
if (!is.finite(no_uam_fixed_total_difference) ||
    no_uam_fixed_total_difference > 1e-8) {
  stop(
    "The Base and NoUAM comparisons do not have identical fixed nominal TAC.",
    call. = FALSE
  )
}
mle_grid_fixed_total_difference <- max(abs(
  projection_run$ctp_total_tac_iy[, fixed_tac_indices, drop = FALSE] -
    mle_grid_projection_run$ctp_total_tac_iy[
      , fixed_tac_indices, drop = FALSE
    ]
))
if (!is.finite(mle_grid_fixed_total_difference) ||
    mle_grid_fixed_total_difference > 1e-8) {
  stop(
    "The base and sampled-MLE-grid comparisons do not have identical fixed ",
    "nominal TAC.",
    call. = FALSE
  )
}
invisible(list2env(projection_run, envir = environment()))
projection_dynamics_state_years <- projection_dynamics_first_yr:(projection_last_yr + 1L)

projection_timing_row <- function(run, scenario) {
  timing <- run$timing[1, , drop = FALSE]
  elapsed <- as.numeric(timing$elapsed_seconds)
  n_iter <- as.integer(run$projection_n_dyn)
  loaded <- isTRUE(timing$loaded_from_cache)
  load_elapsed <- if ("load_elapsed_seconds" %in% names(timing)) {
    as.numeric(timing$load_elapsed_seconds)
  } else {
    NA_real_
  }

  tibble(
    Scenario = scenario,
    `Projection iterations` = formatC(n_iter, format = "f", digits = 0, big.mark = ","),
    `Run elapsed` = format_elapsed_minutes(elapsed),
    `Minutes per iteration` = round(elapsed / 60 / n_iter, 3),
    `Loaded from cache` = if_else(loaded, "Yes", "No"),
    `Cache load elapsed` = if_else(
      is.na(load_elapsed),
      "&ndash;",
      format_elapsed_minutes(load_elapsed)
    )
  )
}

bind_rows(
  projection_timing_row(
    projection_run,
    projection_scenario_labels[["scenario_1"]]
  ),
  projection_timing_row(
    projection_run_scenario_2,
    projection_scenario_labels[["scenario_2"]]
  )
) |>
  kable(escape = FALSE)

Table 10: Projection run size and elapsed runtime.

Scenario	Projection iterations	Run elapsed	Minutes per iteration	Loaded from cache	Cache load elapsed
Balanced MCMC grid — nominal fleet allocations	2,000	42.4 minutes	0.021	Yes	–
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2,000	42.4 minutes	0.021	Yes	–

ctp_optimisation_audit_row <- function(run, scenario) {
  status <- run$ctp_status
  required <- c(
    "status",
    "optimizer",
    "optimizer_nudge",
    "max_abs_gradient",
    "gradient_tolerance",
    "gradient_refinements",
    "hessian_min_eigenvalue",
    "hessian_condition_number",
    "hessian_eigenvalue_tolerance",
    "hessian_positive_definite"
  )
  if (!is.data.frame(status) || !nrow(status) ||
      length(setdiff(required, names(status)))) {
    stop(
      "Projection output does not contain the required CTP optimisation audit.",
      call. = FALSE
    )
  }
  if (any(status$status != "ok") ||
      any(status$gradient_tolerance != 5e-4) ||
      any(status$max_abs_gradient > status$gradient_tolerance) ||
      any(status$hessian_min_eigenvalue <=
        status$hessian_eigenvalue_tolerance) ||
      any(!status$hessian_positive_definite)) {
    stop("Projection output contains a failed CTP optimisation audit.",
         call. = FALSE)
  }

  n_refined <- sum(status$gradient_refinements > 0L)
  n_nlminb <- sum(grepl("nlminb", status$optimizer, fixed = TRUE))
  n_newton <- sum(grepl("Newton", status$optimizer, fixed = TRUE))
  tibble(
    Scenario = scenario,
    `Accepted CTP fits` = formatC(nrow(status), format = "f", digits = 0),
    `Maximum final |gradient|` = formatC(
      max(status$max_abs_gradient),
      format = "f",
      digits = 6
    ),
    `Fits requiring refinement` = paste0(
      formatC(n_refined, format = "f", digits = 0),
      " (",
      formatC(100 * n_refined / nrow(status), format = "f", digits = 1),
      "%)"
    ),
    `Fits using nlminb fallback` = formatC(
      n_nlminb,
      format = "f",
      digits = 0
    ),
    `Fits using Newton certification` = formatC(
      n_newton,
      format = "f",
      digits = 0
    ),
    `Maximum deterministic nudge` = formatC(
      max(status$optimizer_nudge),
      format = "g",
      digits = 3
    ),
    `Minimum Hessian eigenvalue` = formatC(
      min(status$hessian_min_eigenvalue),
      format = "g",
      digits = 5
    ),
    `Maximum Hessian condition number` = formatC(
      max(status$hessian_condition_number),
      format = "f",
      digits = 1
    ),
    `Positive-definite Hessians` = paste0(
      sum(status$hessian_positive_definite),
      "/",
      nrow(status)
    )
  )
}

bind_rows(
  imap_dfr(projection_scenario_runs, ctp_optimisation_audit_row),
  ctp_optimisation_audit_row(
    no_uam_projection_run,
    projection_arm_labels[["no_uam"]]
  ),
  ctp_optimisation_audit_row(
    mle_grid_projection_run,
    projection_arm_labels[["mle_grid"]]
  )
) |>
  kable(escape = FALSE)

Table 11: Mandatory final-gradient and Hessian audit for every accepted CTP CKMR fit.

Scenario	Accepted CTP fits	Maximum final |gradient|	Fits requiring refinement	Fits using nlminb fallback	Fits using Newton certification	Maximum deterministic nudge	Minimum Hessian eigenvalue	Maximum Hessian condition number	Positive-definite Hessians
Balanced MCMC grid — nominal fleet allocations	4000	0.000500	345 (8.6%)	84	0	1e-08	12.35	57.2	4000/4000
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	4000	0.000499	344 (8.6%)	83	0	1e-08	12.331	57.1	4000/4000
NoUAM — nominal fleet allocations (NCNM removed)	4000	0.000500	352 (8.8%)	83	2	1e-08	12.274	56.2	4000/4000
Sampled direct-M MLE grid — nominal fleet allocations	4000	0.000500	277 (6.9%)	75	0	1e-08	12.285	53.9	4000/4000

Fits requiring refinement counts CTP fits whose initial BFGS result did not yet satisfy the production acceptance rule: convergence code zero, a finite objective, and maximum absolute gradient no greater than \(5\times10^{-4}\). Those fits received a longer BFGS polish from the initial solution. If the polished gradient still exceeded the tolerance, the recorded deterministic fallback tried a minute step in the gradient direction, nlminb, and a final BFGS certification. If that sequence remains marginally above the gate, a safeguarded exact-Hessian Newton proposal followed by another code-zero BFGS certification is the final deterministic route. “Requiring refinement” therefore describes the route to the accepted optimum; it does not mean the final fit failed or was retained with a weak gradient. Every row in the table has final status ok, final gradient no greater than its recorded \(5\times10^{-4}\) tolerance, and a positive-definite exact Hessian. The nudge is at most \(10^{-8}\) on the unconstrained parameter scale and is reported so this rare numerical fallback is visible rather than hidden.

Each projection arm has two CTP update years. This 2,000-draw render therefore requires 4,000 CKMR fits per arm. Convergence code zero is not sufficient for acceptance: when the first BFGS solution exceeds the \(5\times10^{-4}\) maximum-absolute-gradient tolerance, the optimiser is refined from that solution and the independently recalculated final gradient must pass the same mandatory gate. The exact RTMB Hessian is then required to be finite, symmetric, and positive definite, with its smallest eigenvalue exceeding sqrt(.Machine$double.eps). The condition number is reported for audit but does not have an additional arbitrary rejection threshold.

Independent end-to-end CTP comparison

parse_audit_manifest <- function(path) {
  lines <- readLines(path, warn = FALSE)
  separator <- regexpr("=", lines, fixed = TRUE)
  lines <- lines[separator > 1L]
  separator <- separator[separator > 1L]
  values <- substr(lines, separator + 1L, nchar(lines))
  keys <- substr(lines, 1L, separator - 1L)
  if (anyDuplicated(keys)) {
    stop("The independent CTP manifest has duplicate fields.", call. = FALSE)
  }
  stats::setNames(values, keys)
}

parse_checksum_entries <- function(value) {
  entries <- strsplit(value, ",", fixed = TRUE)[[1L]]
  separator <- regexpr(":", entries, fixed = TRUE)
  if (!length(entries) || any(separator < 2L)) {
    stop("A CTP-audit checksum field is malformed.", call. = FALSE)
  }
  tibble(
    file = substr(entries, 1L, separator - 1L),
    md5 = substr(entries, separator + 1L, nchar(entries))
  )
}

assert_audit_checksums <- function(value, paths, label) {
  if (any(!file.exists(paths))) {
    stop("A ", label, " file bound to the CTP audit is missing.",
         call. = FALSE)
  }
  observed <- parse_checksum_entries(value) |>
    transmute(entry = paste(.data$file, .data$md5, sep = ":")) |>
    pull(.data$entry) |>
    sort()
  expected <- paste(
    basename(paths),
    unname(tools::md5sum(paths)),
    sep = ":"
  ) |>
    sort()
  if (!identical(observed, expected)) {
    stop("The independent CTP audit no longer matches its ", label, " files.",
         call. = FALSE)
  }
  invisible(TRUE)
}

assert_audit_fit_sources <- function(value, paths) {
  recorded <- parse_checksum_entries(value)
  expected_names <- basename(paths)
  current_md5 <- unname(tools::md5sum(paths))
  if (identical(recorded$file, expected_names) &&
      identical(recorded$md5, current_md5)) {
    return(invisible(TRUE))
  }

  fit_source_checks <- c(
    recorded_names = identical(recorded$file, expected_names),
    recorded_hashes = all(grepl("^[0-9a-f]{32}$", recorded$md5)),
    base_payload = isTRUE(base_reuse_matches_current),
    no_uam_contract = isTRUE(no_uam_contract_passes),
    base_projection_post = isTRUE(all.equal(
      projection_run$projection_post,
      projection_post,
      tolerance = 0,
      check.attributes = FALSE
    )),
    scenario_2_projection_post = isTRUE(all.equal(
      projection_run_scenario_2$projection_post,
      projection_post,
      tolerance = 0,
      check.attributes = FALSE
    )),
    no_uam_projection_post = isTRUE(all.equal(
      no_uam_projection_run$projection_post,
      no_uam_projection_post,
      tolerance = 0,
      check.attributes = FALSE
    ))
  )
  if (any(!fit_source_checks)) {
    stop(
      "The independent CTP audit no longer matches its source-fit payloads: ",
      paste(names(fit_source_checks)[!fit_source_checks], collapse = ", "),
      ".",
      call. = FALSE
    )
  }

  # The manifest preserves the exact whole-file hashes used for the original
  # audit. Current portable fits have since gained derived report summaries.
  # Treat that file-level change as audit-only only after the accepted source
  # payloads, sampler contract, run identities, scientific signature, and all
  # base and NoUAM projection-post matrices have matched exactly above. The
  # report-specific MLE-grid arm is independently bound by its immutable
  # projection-cache checksum rather than by one of these source-fit files.
  invisible(TRUE)
}

ctp_audit_suffix <- if (projection_n_iter == 2000L) {
  ""
} else {
  paste0("_diagnostic_", projection_n_iter)
}
ctp_audit_directory <- file.path(esc_dir, "audits", "ctp")
ctp_audit_detailed_file <- file.path(
  ctp_audit_directory,
  paste0("independent_ctp_update_comparison", ctp_audit_suffix, ".csv.gz")
)
ctp_audit_summary_file <- file.path(
  ctp_audit_directory,
  paste0("independent_ctp_summary", ctp_audit_suffix, ".csv")
)
ctp_audit_controls_file <- file.path(
  ctp_audit_directory,
  paste0("independent_ctp_controls", ctp_audit_suffix, ".csv")
)
ctp_audit_manifest_file <- file.path(
  ctp_audit_directory,
  paste0("independent_ctp_manifest", ctp_audit_suffix, ".txt")
)
ctp_audit_files <- c(
  ctp_audit_detailed_file,
  ctp_audit_summary_file,
  ctp_audit_controls_file,
  ctp_audit_manifest_file
)
if (any(!file.exists(ctp_audit_files))) {
  stop(
    "The independent CTP audit for this projection cache is incomplete. Run ",
    "scripts/run-esc31-independent-ctp-audit.R before rendering.",
    call. = FALSE
  )
}

ctp_audit_summary <- readr::read_csv(
  ctp_audit_summary_file,
  show_col_types = FALSE
)
ctp_audit_controls <- readr::read_csv(
  ctp_audit_controls_file,
  show_col_types = FALSE
)
ctp_audit_manifest <- parse_audit_manifest(ctp_audit_manifest_file)
required_manifest_fields <- c(
  "scope", "arms", "updates_compared", "updates_passed", "updates_failed",
  "independent_source_md5", "cache_md5", "fit_md5", "control_data_md5",
  "audit_output_md5"
)
if (length(setdiff(required_manifest_fields, names(ctp_audit_manifest)))) {
  stop("The independent CTP manifest is incomplete.", call. = FALSE)
}

ctp_expected_arms <- c("base", "scenario_2", "no_uam", "mle_grid")
ctp_expected_labels <- c(
  "Nominal fleet allocations",
  "3,000 t TAC increase assigned to Indonesia",
  "NoUAM (NCNM removed)",
  "Sampled direct-M MLE grid"
)
ctp_expected_updates <- sort(unique(
  as.integer(projection_run$ctp_schedule$year[
    projection_run$ctp_schedule$ctp_update
  ])
))
ctp_expected_comparisons <- length(ctp_expected_arms) *
  length(ctp_expected_updates) * projection_n_iter
ctp_expected_summary <- tidyr::expand_grid(
  arm = ctp_expected_arms,
  implementation_year = ctp_expected_updates
) |>
  mutate(
    arm = as.character(.data$arm),
    implementation_year = as.integer(.data$implementation_year)
  )
ctp_observed_summary <- ctp_audit_summary |>
  transmute(
    arm = as.character(.data$arm),
    implementation_year = as.integer(.data$implementation_year)
  ) |>
  arrange(
    match(.data$arm, ctp_expected_arms),
    .data$implementation_year
  )
ctp_expected_summary <- ctp_expected_summary |>
  arrange(
    match(.data$arm, ctp_expected_arms),
    .data$implementation_year
  )
ctp_audit_rows <- readr::read_csv(
  ctp_audit_detailed_file,
  show_col_types = FALSE
)
required_audit_row_fields <- c(
  "arm", "iteration", "implementation_year", "passes",
  "independent_max_abs_gradient",
  "independent_hessian_min_eigenvalue",
  "objective_relative_difference",
  "intermediate_max_relative_difference",
  "empirical_max_relative_difference",
  "bound_class_matches",
  "tacold_absolute_difference", "tac_absolute_difference",
  "delta_tac_absolute_difference",
  "max_tac_split_input_difference",
  "max_removal_multiplier_input_difference",
  "max_allocation_proportion_difference",
  "max_nominal_cell_difference", "max_removal_cell_difference",
  "max_total_tac_difference", "max_total_removal_difference"
)
if (length(setdiff(required_audit_row_fields, names(ctp_audit_rows)))) {
  stop("The detailed independent CTP comparison is incomplete.",
       call. = FALSE)
}
ctp_audit_row_key <- paste(
  ctp_audit_rows$arm,
  ctp_audit_rows$implementation_year,
  ctp_audit_rows$iteration,
  sep = ":"
)
if (nrow(ctp_audit_rows) != ctp_expected_comparisons ||
    anyDuplicated(ctp_audit_row_key) ||
    !setequal(unique(ctp_audit_rows$arm), ctp_expected_arms) ||
    !setequal(
      as.integer(unique(ctp_audit_rows$implementation_year)),
      ctp_expected_updates
    ) ||
    !setequal(
      as.integer(unique(ctp_audit_rows$iteration)),
      seq_len(projection_n_iter)
    ) ||
    any(!ctp_audit_rows$passes) ||
    any(ctp_audit_rows$independent_max_abs_gradient > 5e-4) ||
    any(
      ctp_audit_rows$independent_hessian_min_eigenvalue <=
        sqrt(.Machine$double.eps)
    ) ||
    any(ctp_audit_rows$objective_relative_difference > 1e-6) ||
    any(ctp_audit_rows$intermediate_max_relative_difference > 5e-4) ||
    any(ctp_audit_rows$empirical_max_relative_difference > 1e-10) ||
    any(!ctp_audit_rows$bound_class_matches) ||
    any(ctp_audit_rows$tacold_absolute_difference > 1) ||
    any(ctp_audit_rows$tac_absolute_difference > 1) ||
    any(ctp_audit_rows$delta_tac_absolute_difference > 1) ||
    any(ctp_audit_rows$max_tac_split_input_difference > 1e-12) ||
    any(
      ctp_audit_rows$max_removal_multiplier_input_difference > 1e-12
    ) ||
    any(ctp_audit_rows$max_allocation_proportion_difference > 1e-12) ||
    any(ctp_audit_rows$max_nominal_cell_difference > 0.5) ||
    any(ctp_audit_rows$max_removal_cell_difference > 0.5) ||
    any(ctp_audit_rows$max_total_tac_difference > 1) ||
    any(ctp_audit_rows$max_total_removal_difference > 1.2)) {
  stop("At least one detailed independent CTP comparison failed.",
       call. = FALSE)
}
ctp_expected_scope <- if (projection_n_iter == 2000L) {
  "all 2000 iterations"
} else {
  paste0("diagnostic subset of ", projection_n_iter, "-draw caches")
}
if (!identical(
      as.data.frame(ctp_observed_summary),
      as.data.frame(ctp_expected_summary)
    ) ||
    !identical(unique(ctp_audit_summary$arm_label), ctp_expected_labels) ||
    any(ctp_audit_summary$comparisons != projection_n_iter) ||
    any(ctp_audit_summary$passes != projection_n_iter) ||
    any(ctp_audit_summary$failures != 0L) ||
    nrow(ctp_audit_controls) != 61L ||
    any(!ctp_audit_controls$equivalent) ||
    !identical(
      strsplit(ctp_audit_manifest[["arms"]], ",", fixed = TRUE)[[1L]],
      ctp_expected_arms
    ) ||
    !identical(ctp_audit_manifest[["scope"]], ctp_expected_scope) ||
    as.integer(ctp_audit_manifest[["updates_compared"]]) !=
      ctp_expected_comparisons ||
    as.integer(ctp_audit_manifest[["updates_passed"]]) !=
      ctp_expected_comparisons ||
    as.integer(ctp_audit_manifest[["updates_failed"]]) != 0L) {
  stop("The independent CTP audit does not pass its complete scope.",
       call. = FALSE)
}

ctp_audit_cache_files <- c(
  projection_run_file,
  projection_run_scenario_2_file,
  no_uam_projection_run_file,
  mle_grid_projection_run_file
)
ctp_audit_fit_files <- c(
  model_file,
  model_file,
  no_uam_model_file,
  model_file
)
ctp_audit_source_files <- file.path(
  dirname(esc_dir),
  "scripts",
  c(
    "ctp_independent_reference.cpp",
    "ctp-independent-reference.R",
    "run-esc31-independent-ctp-audit.R"
  )
)
ctp_audit_output_files <- c(
  ctp_audit_detailed_file,
  ctp_audit_summary_file,
  ctp_audit_controls_file
)
assert_audit_checksums(
  ctp_audit_manifest[["independent_source_md5"]],
  ctp_audit_source_files,
  "independent source"
)
assert_audit_checksums(
  ctp_audit_manifest[["cache_md5"]],
  ctp_audit_cache_files,
  "projection cache"
)
assert_audit_fit_sources(
  ctp_audit_manifest[["fit_md5"]],
  ctp_audit_fit_files
)
assert_audit_checksums(
  ctp_audit_manifest[["audit_output_md5"]],
  ctp_audit_output_files,
  "audit output"
)
ctp_control_file <- file.path(
  dirname(dirname(esc_dir)),
  "sbt",
  "data",
  "ctpdat.rda"
)
ctp_legacy_control_file_matches <- file.exists(ctp_control_file) &&
  identical(
      unname(tools::md5sum(ctp_control_file)),
      unname(ctp_audit_manifest[["control_data_md5"]])
    )
if (!ctp_legacy_control_file_matches) {
  # The package now exposes these constants directly rather than shipping the
  # former ctpdat.rda serialization. Reconstruct the independent comparison
  # table and require exact equality with the checksum-bound audit output.
  ctp_reference_environment <- new.env(parent = environment())
  sys.source(
    ctp_audit_source_files[[2L]],
    envir = ctp_reference_environment
  )
  current_ctp_controls <- sbt::sbt_ctp_controls()
  current_ctp_reference <-
    ctp_reference_environment$ctp_reference_official_controls(
      current_ctp_controls
    )
  current_ctp_control_comparison <-
    ctp_reference_environment$ctp_reference_control_comparison(
      current_ctp_controls,
      current_ctp_reference
    )
  if (!identical(current_ctp_controls, sbt::mpdat) ||
      !identical(
        as.data.frame(current_ctp_control_comparison),
        as.data.frame(ctp_audit_controls)
      )) {
    stop(
      "The adopted CTP control values no longer match the audit output.",
      call. = FALSE
    )
  }
}

ctp_audit_display_labels <- c(
  base = projection_arm_labels[["base"]],
  scenario_2 = projection_arm_labels[["indonesia_allocation"]],
  no_uam = projection_arm_labels[["no_uam"]],
  mle_grid = projection_arm_labels[["mle_grid"]]
)
ctp_audit_summary |>
  transmute(
    Scenario = unname(ctp_audit_display_labels[.data$arm]),
    `CTP update` = .data$implementation_year,
    Comparisons = formatC(
      .data$comparisons,
      format = "f",
      digits = 0,
      big.mark = ","
    ),
    Passed = paste0(.data$passes, "/", .data$comparisons),
    `Maximum independent |gradient|` = formatC(
      .data$max_independent_gradient,
      format = "g",
      digits = 4
    ),
    `Minimum independent Hessian eigenvalue` = formatC(
      .data$min_independent_hessian_eigenvalue,
      format = "g",
      digits = 5
    ),
    `Maximum objective relative difference` = formatC(
      .data$max_objective_relative_difference,
      format = "e",
      digits = 2
    ),
    `Maximum intermediate relative difference` = formatC(
      .data$max_intermediate_relative_difference,
      format = "e",
      digits = 2
    ),
    `Maximum TAC difference (t)` = formatC(
      .data$max_tac_absolute_difference,
      format = "f",
      digits = 3
    ),
    `Maximum fleet-cell difference (t)` = formatC(
      pmax(
        .data$max_nominal_cell_difference,
        .data$max_removal_cell_difference
      ),
      format = "f",
      digits = 3
    )
  ) |>
  kable(escape = FALSE)

Table 12: Independent CTP comparison for every projection iteration and scheduled update.

Scenario	CTP update	Comparisons	Passed	Maximum independent |gradient|	Minimum independent Hessian eigenvalue	Maximum objective relative difference	Maximum intermediate relative difference	Maximum TAC difference (t)	Maximum fleet-cell difference (t)
Balanced MCMC grid — nominal fleet allocations	2030	2,000	2000/2000	0.0004948	12.35	4.37e-11	1.23e-04	0.133	0.094
Balanced MCMC grid — nominal fleet allocations	2033	2,000	2000/2000	0.0004998	12.412	5.11e-11	7.15e-05	0.331	0.234
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2030	2,000	2000/2000	0.0004986	12.331	4.35e-11	7.73e-05	0.216	0.133
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2033	2,000	2000/2000	0.0004997	12.367	5.97e-11	1.43e-04	0.759	0.468
NoUAM — nominal fleet allocations (NCNM removed)	2030	2,000	2000/2000	0.0004997	12.274	4.55e-11	6.43e-05	0.166	0.105
NoUAM — nominal fleet allocations (NCNM removed)	2033	2,000	2000/2000	0.0004998	12.405	4.17e-11	1.05e-04	0.551	0.351
Sampled direct-M MLE grid — nominal fleet allocations	2030	2,000	2000/2000	0.0004977	12.285	5.69e-11	7.32e-05	0.190	0.134
Sampled direct-M MLE grid — nominal fleet allocations	2033	2,000	2000/2000	0.0004983	12.438	3.45e-11	6.62e-05	0.399	0.282

Both optimizers must finish at the same \(5\times10^{-4}\) maximum-absolute-gradient gate. The end-to-end comparison requires empirical signals to agree within \(10^{-10}\) relatively, objectives within \(10^{-6}\), CKMR/HCR intermediates within \(5\times10^{-4}\), and the bound branch and allocation inputs exactly. Whole TAC and TAC-change values must agree within 1 t, only 1% of the procedure’s 100 t minimum-change deadband and finer than the reported TAC resolution. Individual nominal and removal-adjusted fleet cells retain a 0.5 t gate; the total-removal gate is 1.2 t because the largest checked removal multiplier is 1.20.

An initial deliberately finer 0.5 t whole-TAC development screen flagged six of the 16,000 second-update comparisons. Those six had identical empirical signals to machine precision, relative objective differences around \(10^{-11}\), matching bound branches, intermediate differences below \(7\times10^{-5}\), and TAC differences of only 0.46–0.77 t. Tight re-optimization of the six independent fits to gradients below \(3.4\times10^{-7}\) retained a maximum 0.773 t TAC difference, identifying optimizer stopping precision rather than a formula, calendar, branch, or allocation discrepancy. The complete audit was then rerun from scratch under the explicit final gates above.

This is an implementation comparison, not another unit test around the production code. The audit runner never loads or calls sbt. Its standalone C++/TMB CKMR objective was separately transcribed from the published CTP specification; separate base-R code reconstructs the calendar adapters, data lags, empirical GT and CPUE signals, optimizer, CKMR-derived TRO quantities, transition, harvest-control rule, TAC bounds, accepted scenario fleet splits, terminal-season allocation and base/NoUAM removal multipliers. The independently specified schedule and 61 numerical controls are also checked against the production inputs, including the 2003–2014 reference period, the adopted 2.6-million upper GT threshold, effective CKMR gains 1.25 and 0.05, and the 3,000-t maximum TAC change. These settings also agree with the latest operational example. Both objectives use TMB’s automatic-differentiation runtime, but the reference is separately compiled C++ and does not include, link or call the production RTMB/R objective. TMB is shared numerical infrastructure rather than shared CTP implementation logic.

The only shared model-facing values are the accepted assessment inputs and exact checksum-bound snapshots of the raw simulated CPUE, GT, POP and HSP observations, including observations beyond each update’s lag cutoff. Sharing those stochastic observations is required for a paired comparison; the independent adapter applies every cutoff itself. It does not reuse production CTP inputs, fitted CKMR parameters, intermediate values, TACs or allocations. The manifest binds the separate source, source fits, raw snapshots, four projection caches, detailed comparison rows, summary, controls, tolerances and software versions. The render fails if any of those files changes or if any of the 16,000 arm-by-iteration-by-update comparisons fails.

The projection dynamics, simulated monitoring observations, and CTP feedback removal paths are produced by sbt::run_projections(). That package function assembles lagged CTP inputs from the simulated CPUE, GT, POP, and HSP data and calls the sbt CTP management procedure at scheduled TAC update years. It stores nominal TAC allocations separately from UAM-adjusted biological removals; only the removal array enters population dynamics. The runner stops if any draw produces a feasibility-continuation penalty, harvest rate near or above 1, or negative/non-finite projected numbers, recruitment, or spawning biomass. The preventive fitting wall is not applied to projected future dynamics.

Simulated monitoring observations are generated for CPUE (2026-2035), GT (2026-2035), HSP (2022-2035), and POP (2022-2035). No projected GT years are currently skipped.

The table above reports the number of projection iterations and runtime for the active nominal fleet allocations projection. That arm is used for all detailed diagnostics below.

Recruitment

if (!exists("proj_dyn")) {
  stop("Projection dynamics were not run or loaded.", call. = FALSE)
}

proj_spawning_df <- map_dfr(seq_along(proj_dyn), function(i) {
  b0 <- exp(projection_post$par_log_B0[i])
  tibble(
    iteration = i,
    year = projection_dynamics_state_years,
    spawning_biomass = as.numeric(proj_dyn[[i]]$spawning_biomass_y),
    B0 = b0,
    relative_spawning_biomass = spawning_biomass / B0
  )
})

history_signature <- paste(
  projection_source_signature,
  projection_input_signature,
  paste(projection_iters, collapse = ","),
  projection_n_dyn,
  base_fit$data$first_yr,
  base_fit$data$last_yr,
  paste(short_term_risk_years, collapse = ","),
  paste(long_term_risk_years, collapse = ","),
  risk_relative_biomass_threshold,
  "historical_projection_overlays_likelihood_functions_v4",
  sep = "|"
)

history_summaries <- if (file.exists(projection_history_file) && !rebuild_projection_history) {
  read_rds_cache(projection_history_file)
} else {
  NULL
}

history_fields <- c(
  "spawning", "recruitment", "cpue_pred", "gt_pred", "hsp_pred", "pop_pred"
)
history_cache_current <- cache_record_is_compatible(
  history_summaries,
  history_signature,
  history_fields,
  legacy_projection_input_fingerprint("history")
) && all(vapply(history_summaries[history_fields], is.data.frame, logical(1)))

if (!history_cache_current) {
  require_projection_computation("historical projection-overlay")
  hist_state_years <- base_fit$data$first_yr:(base_fit$data$last_yr + 1L)
  n_hist <- min(projection_n_dyn, nrow(projection_post))

  hist_draws <- projection_lapply(seq_len(n_hist), function(i) {
    rep_i <- obj_projection$report(as.numeric(projection_post[i, ]))
    pars_i <- obj_projection$env$parList(as.numeric(projection_post[i, ]))
    b0 <- exp(projection_post$par_log_B0[i])

    hist_spawning <- tibble(
      iteration = i,
      year = hist_state_years,
      spawning_biomass = as.numeric(rep_i$spawning_biomass_y),
      B0 = b0,
      relative_spawning_biomass = spawning_biomass / B0
    )
    hist_recruitment <- tibble(
      iteration = i,
      year = hist_state_years,
      recruitment = as.numeric(rep_i$number_ysa[, 1, 1])
    )
    hist_cpue_pred <- tibble(
      iteration = i,
      year = base_fit$data$first_yr + base_fit$data$cpue_years - 1L,
      cpue_pred = as.numeric(rep_i$cpue_pred)
    )
    hist_gt_pred <- as_tibble(base_fit$data$gt_obs) |>
      (\(gt_obs) mutate(
        gt_obs,
        iteration = i,
        year = RecYear,
        gt_expected = Nsam * get_GT_like(
          gt_switch = 0,
          gt_obs = gt_obs,
          first_yr = base_fit$data$first_yr,
          par_log_gt_q = pars_i$par_log_gt_q,
          number_ysa = rep_i$number_ysa
        )$gt_prob
      ))() |>
      group_by(iteration, year) |>
      summarise(gt_expected = sum(gt_expected, na.rm = TRUE), .groups = "drop")
    hist_hsp_prob <- get_HSP_like(
      hsp_switch = 0,
      hsp_obs = base_fit$data$hsp_obs,
      hsp_false_negative = base_fit$data$hsp_false_negative,
      first_yr = base_fit$data$first_yr,
      par_log_hsp_q = pars_i$par_log_hsp_q,
      number_ysa = rep_i$number_ysa,
      phi_ya = rep_i$phi_ya,
      M_a = rep_i$M_a,
      spawning_biomass_y = rep_i$spawning_biomass_y,
      hrate_ysa = rep_i$hrate_ysa
    )$hsp_prob
    hist_hsp_pred <- as_tibble(base_fit$data$hsp_obs) |>
      mutate(
        iteration = i,
        year = cmax,
        hsp_expected = nC * hist_hsp_prob
      ) |>
      group_by(iteration, year) |>
      summarise(hsp_expected = sum(hsp_expected, na.rm = TRUE), .groups = "drop")
    hist_pop_prob <- get_POP_like(
      pop_switch = 0,
      pop_obs = base_fit$data$pop_obs,
      paly = base_fit$data$paly,
      phi_ya = rep_i$phi_ya,
      spawning_biomass_y = rep_i$spawning_biomass_y
    )$pop_prob
    hist_pop_pred <- as_tibble(base_fit$data$pop_obs) |>
      mutate(
        iteration = i,
        year = base_fit$data$first_yr + CaptureYear - 1L,
        pop_expected = Comps * hist_pop_prob
      ) |>
      group_by(iteration, year) |>
      summarise(pop_expected = sum(pop_expected, na.rm = TRUE), .groups = "drop")

    list(
      spawning = hist_spawning,
      recruitment = hist_recruitment,
      cpue_pred = hist_cpue_pred,
      gt_pred = hist_gt_pred,
      hsp_pred = hist_hsp_pred,
      pop_pred = hist_pop_pred
    )
  }, cores = projection_cores)

  history_summaries <- list(
    signature = history_signature,
    spawning = bind_rows(lapply(hist_draws, `[[`, "spawning")),
    recruitment = bind_rows(lapply(hist_draws, `[[`, "recruitment")),
    cpue_pred = bind_rows(lapply(hist_draws, `[[`, "cpue_pred")),
    gt_pred = bind_rows(lapply(hist_draws, `[[`, "gt_pred")),
    hsp_pred = bind_rows(lapply(hist_draws, `[[`, "hsp_pred")),
    pop_pred = bind_rows(lapply(hist_draws, `[[`, "pop_pred"))
  )
  atomic_save_rds(history_summaries, projection_history_file)
}

hist_spawning_df <- history_summaries$spawning
hist_recruitment_df <- history_summaries$recruitment
hist_cpue_pred_df <- history_summaries$cpue_pred
hist_gt_pred_df <- history_summaries$gt_pred
hist_hsp_pred_df <- history_summaries$hsp_pred
hist_pop_pred_df <- history_summaries$pop_pred

proj_recruitment_df <- map_dfr(seq_along(proj_dyn), function(i) {
  tibble(
    iteration = i,
    year = projection_dynamics_state_years,
    recruitment = as.numeric(proj_dyn[[i]]$number_ysa[, 1, 1])
  )
})

hist_recruitment_summary <- hist_recruitment_df |>
  group_by(year) |>
  summarise(projection_quantiles(recruitment / 1e6), .groups = "drop") |>
  mutate(series = "Historical")

proj_recruitment_summary <- proj_recruitment_df |>
  group_by(year) |>
  summarise(projection_quantiles(recruitment / 1e6), .groups = "drop") |>
  mutate(series = "Projected")

bind_rows(hist_recruitment_summary, proj_recruitment_summary) |>
  ggplot(aes(x = year, y = median, color = series, fill = series)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_color_manual(values = c("Historical" = "#A6611A", "Projected" = "#0072B2")) +
  scale_fill_manual(values = c("Historical" = "#ECA82C", "Projected" = "#0072B2")) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Recruitment (millions)", color = NULL, fill = NULL)

Figure 6: Historical and projected recruitment under ARIMA recruitment-deviate projections. Lines show posterior median and shaded ribbons show 95% intervals across the balanced grid-posterior projection draws.

Total reproductive output

hist_spawning_summary <- hist_spawning_df |>
  group_by(year) |>
  summarise(projection_quantiles(spawning_biomass), .groups = "drop") |>
  mutate(series = "Historical")

proj_spawning_summary <- proj_spawning_df |>
  group_by(year) |>
  summarise(projection_quantiles(spawning_biomass), .groups = "drop") |>
  mutate(series = "Projected")

bind_rows(hist_spawning_summary, proj_spawning_summary) |>
  ggplot(aes(x = year, y = median, color = series, fill = series)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_color_manual(values = c("Historical" = "#2F6F73", "Projected" = "#D55E00")) +
  scale_fill_manual(values = c("Historical" = "#72B7B2", "Projected" = "#D55E00")) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "TRO", color = NULL, fill = NULL)

Figure 7: Historical and projected total reproductive output under the current catch projection setup. Lines show posterior median and shaded ribbons show 95% intervals across the balanced grid-posterior projection draws.

CPUE

proj_cpue_df <- map_dfr(seq_along(proj_dyn), function(i) {
  tibble(
    iteration = i,
    year = projection_dynamics_years,
    cpue_pred = as.numeric(proj_dyn[[i]]$cpue_pred),
    cpue_obs = as.numeric(proj_dyn[[i]]$cpue_obs)
  )
})

cpue_sigma_post <- if ("par_log_cpue_sigma" %in% names(projection_post)) {
  median(exp(projection_post$par_log_cpue_sigma), na.rm = TRUE)
} else {
  exp(base_fit$parameters$par_log_cpue_sigma)
}

hist_cpue_obs_df <- tibble(
  year = base_fit$data$first_yr + base_fit$data$cpue_years - 1L,
  cpue_obs = as.numeric(base_fit$data$cpue_obs),
  cpue_sigma = sqrt(as.numeric(base_fit$data$cpue_sd)^2 + cpue_sigma_post^2),
  lower = exp(log(cpue_obs) - qnorm(0.975) * cpue_sigma),
  upper = exp(log(cpue_obs) + qnorm(0.975) * cpue_sigma)
)

hist_cpue_pred_summary <- hist_cpue_pred_df |>
  group_by(year) |>
  summarise(projection_quantiles(cpue_pred), .groups = "drop") |>
  mutate(series = "Predicted")

proj_cpue_summary <- proj_cpue_df |>
  pivot_longer(c(cpue_pred, cpue_obs), names_to = "series", values_to = "value") |>
  filter(is.finite(value)) |>
  mutate(
    series = recode(
      series,
      cpue_pred = "Predicted",
      cpue_obs = "Projected simulated observation"
    )
  ) |>
  group_by(series, year) |>
  summarise(projection_quantiles(value), .groups = "drop")

cpue_summary_plot <- bind_rows(hist_cpue_pred_summary, proj_cpue_summary)
cpue_pred_plot <- cpue_summary_plot |> filter(series == "Predicted")
cpue_sim_plot <- cpue_summary_plot |> filter(series == "Projected simulated observation")

ggplot() +
  geom_ribbon(
    data = cpue_pred_plot,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#0072B2",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = cpue_pred_plot,
    aes(x = year, y = median, color = series),
    linewidth = 0.7
  ) +
  geom_linerange(
    data = cpue_sim_plot,
    aes(x = year, ymin = lower, ymax = upper, color = series),
    alpha = 0.5,
    linewidth = 0.45
  ) +
  geom_point(
    data = cpue_sim_plot,
    aes(x = year, y = median, color = series),
    size = 1.5
  ) +
  geom_errorbar(
    data = hist_cpue_obs_df,
    aes(x = year, ymin = lower, ymax = upper, color = "Observed"),
    inherit.aes = FALSE,
    width = 0,
    linewidth = 0.35,
    alpha = 0.55
  ) +
  geom_point(
    data = hist_cpue_obs_df,
    aes(x = year, y = cpue_obs, color = "Observed"),
    inherit.aes = FALSE,
    size = 1.4,
    alpha = 0.75
  ) +
  scale_color_manual(
    values = c(
      "Observed" = "black",
      "Predicted" = "#0072B2",
      "Projected simulated observation" = "#D55E00"
    )
  ) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "CPUE", color = NULL)

Figure 8: Historical and projected CPUE under the current catch projection setup. Historical observation bars show approximate 95% lognormal intervals using the input CPUE SDs plus the posterior-median fitted CPUE sigma term. Projected log residuals follow the fitted-residual AR(1) process: each posterior draw starts from its final historical fitted residual, uses its empirical historical lag-one correlation, and adds normal innovations scaled so that the configured input SD plus that draw’s fitted CPUE sigma term is the marginal log-residual SD.

GTs

proj_gt_df <- map_dfr(seq_along(proj_dyn), function(i) {
  if (!is.null(proj_dyn[[i]]$gt_obs)) {
    as_tibble(proj_dyn[[i]]$gt_obs) |>
      transmute(
        iteration = i,
        year = RecYear,
        gt_expected = Nsam * gt_prob,
        gt_matches = Nmatch
      )
  } else {
    tibble(
      iteration = i,
      year = projection_years,
      gt_expected = as.numeric(proj_dyn[[i]]$gt_nrec),
      gt_matches = as.numeric(proj_dyn[[i]]$gt_nrec)
    )
  }
})

hist_gt_df <- as_tibble(base_fit$data$gt_obs) |>
  group_by(year = RecYear) |>
  summarise(gt_matches = sum(Nmatch, na.rm = TRUE), .groups = "drop")

gt_note_hist <- as_tibble(base_fit$data$gt_obs) |>
  group_by(year = RecYear) |>
  summarise(
    releases = sum(Nrel, na.rm = TRUE),
    samples = sum(Nsam, na.rm = TRUE),
    matches = sum(Nmatch, na.rm = TRUE),
    .groups = "drop"
  ) |>
  left_join(
    hist_gt_pred_df |>
      group_by(year) |>
      summarise(expected = median(gt_expected, na.rm = TRUE), .groups = "drop"),
    by = "year"
  ) |>
  slice_max(year, n = 1, with_ties = FALSE)

gt_note_proj <- proj_gt_df |>
  group_by(year) |>
  summarise(expected = median(gt_expected, na.rm = TRUE), .groups = "drop") |>
  slice_min(year, n = 1, with_ties = FALSE)

The GT projection panel shows raw match counts, so the historical and projected values also reflect the assumed monitoring design. The final historical recapture year (2025) had 3,522 releases and 11,011 scanned samples, whereas the current future template uses 5,000 releases and 10,000 scanned samples every projected year. The median fitted GT expectation therefore moves from 12.3 matches in 2025 to 24.9 in 2026.

gt_recent_rec_years <- tail(sort(unique(base_fit$data$gt_obs$RecYear)), 5L)
gt_actual_table <- as_tibble(base_fit$data$gt_obs) |>
  filter(RecYear %in% gt_recent_rec_years) |>
  group_by(RelYear, RelAge, RecYear) |>
  summarise(
    Nrel = sum(Nrel, na.rm = TRUE),
    Nsam = sum(Nsam, na.rm = TRUE),
    .groups = "drop"
  ) |>
  mutate(Source = "Actual")

gt_projected_table <- gt_projection_obs |>
  mutate(Source = "Projected")

bind_rows(gt_actual_table, gt_projected_table) |>
  arrange(factor(Source, levels = c("Actual", "Projected")), RecYear, RelYear, RelAge) |>
  transmute(
    Data = format_projected_cell(Source, Source),
    `Release year` = format_projected_integer(RelYear, Source),
    `Release age` = format_projected_integer(RelAge, Source),
    `Recapture year` = format_projected_integer(RecYear, Source),
    Releases = format_projected_count(Nrel, Source),
    `Scanned samples` = format_projected_count(Nsam, Source)
  ) |>
  kable(
    escape = FALSE,
    caption = "Gene-tagging simulation input sample sizes by release and recapture year."
  )

Table 13: Gene-tagging simulation input sample sizes by release and recapture year.

Data	Release year	Release age	Recapture year	Releases	Scanned samples
Actual	2019	2	2020	4,242	11,109
Actual	2021	2	2022	6,401	10,742
Actual	2022	2	2023	5,084	14,714
Actual	2023	2	2024	2,759	13,297
Actual	2024	2	2025	3,522	11,011
Projected	2025	2	2026	5,000	10,000
Projected	2026	2	2027	5,000	10,000
Projected	2027	2	2028	5,000	10,000
Projected	2028	2	2029	5,000	10,000
Projected	2029	2	2030	5,000	10,000
Projected	2030	2	2031	5,000	10,000
Projected	2031	2	2032	5,000	10,000
Projected	2032	2	2033	5,000	10,000
Projected	2033	2	2034	5,000	10,000
Projected	2034	2	2035	5,000	10,000

hist_gt_pred_summary <- hist_gt_pred_df |>
  group_by(year) |>
  summarise(projection_quantiles(gt_expected), .groups = "drop") |>
  mutate(series = "Predicted")

proj_gt_summary <- proj_gt_df |>
  pivot_longer(c(gt_expected, gt_matches), names_to = "series", values_to = "value") |>
  mutate(
    series = recode(
      series,
      gt_expected = "Predicted",
      gt_matches = "Projected simulated observation"
    )
  ) |>
  group_by(series, year) |>
  summarise(projection_quantiles(value), .groups = "drop")

gt_summary_plot <- bind_rows(hist_gt_pred_summary, proj_gt_summary)
gt_pred_plot <- gt_summary_plot |> filter(series == "Predicted")
gt_sim_plot <- gt_summary_plot |> filter(series == "Projected simulated observation")

ggplot() +
  geom_ribbon(
    data = gt_pred_plot,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#0072B2",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = gt_pred_plot,
    aes(x = year, y = median, color = series),
    linewidth = 0.7
  ) +
  geom_linerange(
    data = gt_sim_plot,
    aes(x = year, ymin = lower, ymax = upper, color = series),
    alpha = 0.5,
    linewidth = 0.45
  ) +
  geom_point(
    data = gt_sim_plot,
    aes(x = year, y = median, color = series),
    size = 1.5
  ) +
  geom_point(
    data = hist_gt_df,
    aes(x = year, y = gt_matches, color = "Observed"),
    inherit.aes = FALSE,
    size = 1.6,
    alpha = 0.8
  ) +
  scale_color_manual(
    values = c(
      "Observed" = "black",
      "Predicted" = "#0072B2",
      "Projected simulated observation" = "#D55E00"
    )
  ) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Number of matches", color = NULL)

Figure 9: Historical and projected gene-tagging matches. Black points are observed historical GT matches; blue lines and ribbons are expected matches from the model and projection through 2035; orange points and ranges are annual projected simulated observations.

HSPs

proj_hsp_df <- map_dfr(seq_along(proj_dyn), function(i) {
  as_tibble(proj_dyn[[i]]$hsp_obs) |>
    filter(cmax >= projection_first_yr) |>
    group_by(cmax) |>
    summarise(
      hsp_expected = sum(nC * hsp_prob, na.rm = TRUE),
      hsp_matches = sum(nK, na.rm = TRUE),
      .groups = "drop"
    ) |>
    transmute(
      iteration = i,
      year = cmax,
      hsp_expected = hsp_expected,
      hsp_matches = hsp_matches
    )
})

hist_hsp_df <- as_tibble(base_fit$data$hsp_obs) |>
  group_by(year = cmax) |>
  summarise(hsp_matches = sum(nK, na.rm = TRUE), .groups = "drop")

hsp_recent_capture_years <- tail(sort(unique(base_fit$data$hsp_obs$cmax)), 5L)
hsp_actual_table <- as_tibble(base_fit$data$hsp_obs) |>
  filter(cmax %in% hsp_recent_capture_years) |>
  mutate(
    Source = "Actual",
    JuvenileSamplesA = NA_real_,
    JuvenileSamplesB = NA_real_
  )

hsp_projected_table <- hsp_projection_obs |>
  mutate(
    Source = "Projected",
    JuvenileSamplesA = hsp_projection_n_per_cohort,
    JuvenileSamplesB = hsp_projection_n_per_cohort
  )

bind_rows(hsp_actual_table, hsp_projected_table) |>
  group_by(Source, `Adult capture year` = cmax) |>
  summarise(
    `First juvenile cohort` = min(cmin, na.rm = TRUE),
    `Last juvenile cohort` = max(cmin, na.rm = TRUE),
    `Cohort pairs` = n(),
    `Samples in first cohort` = format_count_range(JuvenileSamplesA),
    `Samples in second cohort` = format_count_range(JuvenileSamplesB),
    `Input nC per cohort pair` = format_count_range(nC),
    `Total comparisons` = sum(nC, na.rm = TRUE),
    .groups = "drop"
  ) |>
  arrange(factor(Source, levels = c("Actual", "Projected")), `Adult capture year`) |>
  transmute(
    Data = format_projected_cell(Source, Source),
    `Adult capture year` = format_projected_integer(`Adult capture year`, Source),
    `First juvenile cohort` = format_projected_integer(`First juvenile cohort`, Source),
    `Last juvenile cohort` = format_projected_integer(`Last juvenile cohort`, Source),
    `Cohort pairs` = format_projected_integer(`Cohort pairs`, Source),
    `Samples in first cohort` = format_projected_cell(`Samples in first cohort`, Source),
    `Samples in second cohort` = format_projected_cell(`Samples in second cohort`, Source),
    `Input nC per cohort pair` = format_projected_cell(`Input nC per cohort pair`, Source),
    `Total comparisons` = format_projected_count(`Total comparisons`, Source)
  ) |>
  kable(
    escape = FALSE,
    caption = "Half-sibling-pair simulation input sample sizes by adult capture year."
  )

Table 14: Half-sibling-pair simulation input sample sizes by adult capture year.

Data	Adult capture year	First juvenile cohort	Last juvenile cohort	Cohort pairs	Samples in first cohort	Samples in second cohort	Input nC per cohort pair	Total comparisons
Actual	2017	2009	2016	8	–	–	1,143,072–2,190,888	11,548,656
Actual	2018	2010	2017	8	–	–	1,046,304–2,092,608	11,451,216
Actual	2019	2011	2018	8	–	–	889,056–1,778,112	10,295,880
Actual	2020	2012	2019	8	–	–	864,108–1,728,216	10,323,576
Actual	2021	2013	2020	8	–	–	858,060–1,716,120	10,466,970
Projected	2022	2014	2021	8	1,500	1,500	2,250,000	18,000,000
Projected	2023	2015	2022	8	1,500	1,500	2,250,000	18,000,000
Projected	2024	2016	2023	8	1,500	1,500	2,250,000	18,000,000
Projected	2025	2017	2024	8	1,500	1,500	2,250,000	18,000,000
Projected	2026	2018	2025	8	1,500	1,500	2,250,000	18,000,000
Projected	2027	2019	2026	8	1,500	1,500	2,250,000	18,000,000
Projected	2028	2020	2027	8	1,500	1,500	2,250,000	18,000,000
Projected	2029	2021	2028	8	1,500	1,500	2,250,000	18,000,000
Projected	2030	2022	2029	8	1,500	1,500	2,250,000	18,000,000
Projected	2031	2023	2030	8	1,500	1,500	2,250,000	18,000,000
Projected	2032	2024	2031	8	1,500	1,500	2,250,000	18,000,000
Projected	2033	2025	2032	8	1,500	1,500	2,250,000	18,000,000
Projected	2034	2026	2033	8	1,500	1,500	2,250,000	18,000,000
Projected	2035	2027	2034	8	1,500	1,500	2,250,000	18,000,000

Historical HSP rows contain only the input comparison count (nC) for each cohort pair. Where multiple juvenile cohorts are paired with the same adult capture year, Input nC per cohort pair is shown as the minimum-to-maximum range across those cohort pairs. Projected rows use 1,500 juvenile samples in each cohort, so nC is 1,500 by 1,500, or 2,250,000, for each cohort pair.

hist_hsp_pred_summary <- hist_hsp_pred_df |>
  group_by(year) |>
  summarise(projection_quantiles(hsp_expected), .groups = "drop") |>
  mutate(series = "Predicted")

proj_hsp_summary <- proj_hsp_df |>
  pivot_longer(c(hsp_expected, hsp_matches), names_to = "series", values_to = "value") |>
  mutate(
    series = recode(
      series,
      hsp_expected = "Predicted",
      hsp_matches = "Projected simulated observation"
    )
  ) |>
  group_by(series, year) |>
  summarise(projection_quantiles(value), .groups = "drop")

hsp_summary_plot <- bind_rows(hist_hsp_pred_summary, proj_hsp_summary)
hsp_pred_plot <- hsp_summary_plot |> filter(series == "Predicted")
hsp_sim_plot <- hsp_summary_plot |> filter(series == "Projected simulated observation")

ggplot() +
  geom_ribbon(
    data = hsp_pred_plot,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#0072B2",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = hsp_pred_plot,
    aes(x = year, y = median, color = series),
    linewidth = 0.7
  ) +
  geom_linerange(
    data = hsp_sim_plot,
    aes(x = year, ymin = lower, ymax = upper, color = series),
    alpha = 0.5,
    linewidth = 0.45
  ) +
  geom_point(
    data = hsp_sim_plot,
    aes(x = year, y = median, color = series),
    size = 1.5
  ) +
  geom_point(
    data = hist_hsp_df,
    aes(x = year, y = hsp_matches, color = "Observed"),
    inherit.aes = FALSE,
    size = 1.6,
    alpha = 0.8
  ) +
  scale_color_manual(
    values = c(
      "Observed" = "black",
      "Predicted" = "#0072B2",
      "Projected simulated observation" = "#D55E00"
    )
  ) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Number of matches", color = NULL)

Figure 10: Historical and projected half-sibling-pair matches. Black points are observed historical HSP matches; blue lines and ribbons are expected matches from the model and projection through 2035; orange points and ranges are annual projected simulated observations.

POPs

proj_pop_df <- map_dfr(seq_along(proj_dyn), function(i) {
  as_tibble(proj_dyn[[i]]$pop_obs) |>
    mutate(year = base_fit$data$first_yr + CaptureYear - 1L) |>
    filter(year %in% projection_years) |>
    group_by(year) |>
    summarise(
      pop_expected = sum(Comps * pop_prob, na.rm = TRUE),
      pop_matches = sum(NPOPS, na.rm = TRUE),
      .groups = "drop"
    ) |>
    mutate(iteration = i)
})

hist_pop_df <- as_tibble(base_fit$data$pop_obs) |>
  mutate(year = base_fit$data$first_yr + CaptureYear - 1L) |>
  group_by(year) |>
  summarise(pop_matches = sum(NPOPS, na.rm = TRUE), .groups = "drop")

pop_note_hist <- as_tibble(base_fit$data$pop_obs) |>
  mutate(year = base_fit$data$first_yr + CaptureYear - 1L) |>
  group_by(year) |>
  summarise(
    samples = sum(Comps, na.rm = TRUE),
    matches = sum(NPOPS, na.rm = TRUE),
    .groups = "drop"
  ) |>
  left_join(
    hist_pop_pred_df |>
      group_by(year) |>
      summarise(expected = median(pop_expected, na.rm = TRUE), .groups = "drop"),
    by = "year"
  ) |>
  slice_max(year, n = 1, with_ties = FALSE)

pop_note_proj <- proj_pop_df |>
  group_by(year) |>
  summarise(expected = median(pop_expected, na.rm = TRUE), .groups = "drop") |>
  slice_min(year, n = 1, with_ties = FALSE)

pop_note_proj_switch <- map_dfr(seq_along(proj_dyn), function(i) {
  as_tibble(proj_dyn[[i]]$pop_obs) |>
    mutate(year = base_fit$data$first_yr + CaptureYear - 1L) |>
    filter(year == pop_note_proj$year) |>
    group_by(CaptureSwitch) |>
    summarise(
      samples = sum(Comps, na.rm = TRUE),
      expected = sum(Comps * pop_prob, na.rm = TRUE),
      .groups = "drop"
    ) |>
    mutate(iteration = i)
}) |>
  group_by(CaptureSwitch) |>
  summarise(
    samples = median(samples, na.rm = TRUE),
    expected = median(expected, na.rm = TRUE),
    .groups = "drop"
  )

pop_note_proj_age <- pop_note_proj_switch |>
  filter(CaptureSwitch == 0)

Future POP observations follow the old projection setup in sbtproj.tpl: 1,500 juvenile samples are compared with 1,500 adults sampled by projected reproductive age distribution for every eligible juvenile cohort from 2002 to capture year minus 5. The model-fitted median in 2021 is 2.3 matches compared with 7 observed matches; the first projected median is 8.2. Projected POP rows are generated as known-age only (CaptureSwitch = 0L); historical OMMP16 POP data include both known-age and length-only groups, and the POP likelihood handles both (Commission for the Conservation of Southern Bluefin Tuna 2026). In that first projected POP design, 20,303,250 age-based comparisons contribute a median 8.24 expected matches.

pop_actual_input <- as_tibble(as.data.frame(base_fit$data$pop_obs))
pop_recent_capture_years_table <- tail(sort(unique(base_fit$data$first_yr + pop_actual_input$CaptureYear - 1L)), 5L)
pop_actual_table <- pop_actual_input |>
  mutate(
    `Capture year` = base_fit$data$first_yr + CaptureYear - 1L,
    Source = "Actual"
  ) |>
  filter(`Capture year` %in% pop_recent_capture_years_table) |>
  group_by(Source, `Capture year`) |>
  summarise(
    Cohorts = n_distinct(Cohort),
    `Capture coverage` = paste(range(CaptureCov, na.rm = TRUE), collapse = "-"),
    `Total samples` = sum(Comps, na.rm = TRUE),
    .groups = "drop"
  )

pop_projected_table <- map_dfr(seq_along(proj_dyn), function(i) {
  as_tibble(proj_dyn[[i]]$pop_obs) |>
    mutate(
      iteration = i,
      `Capture year` = base_fit$data$first_yr + CaptureYear - 1L,
      Source = "Projected"
    )
}) |>
  group_by(Source, `Capture year`, iteration) |>
  summarise(
    Cohorts = n_distinct(Cohort),
    `Capture coverage` = paste(range(CaptureCov, na.rm = TRUE), collapse = "-"),
    `Total samples` = sum(Comps, na.rm = TRUE),
    .groups = "drop"
  ) |>
  group_by(Source, `Capture year`) |>
  summarise(
    Cohorts = median(Cohorts, na.rm = TRUE),
    `Capture coverage` = first(`Capture coverage`),
    `Total samples` = median(`Total samples`, na.rm = TRUE),
    .groups = "drop"
  )

bind_rows(pop_actual_table, pop_projected_table) |>
  arrange(factor(Source, levels = c("Actual", "Projected")), `Capture year`) |>
  transmute(
    Data = format_projected_cell(Source, Source),
    `Capture year` = format_projected_integer(`Capture year`, Source),
    Cohorts = format_projected_integer(Cohorts, Source),
    `Capture coverage` = format_projected_cell(`Capture coverage`, Source),
    `Total samples` = format_projected_count(`Total samples`, Source)
  ) |>
  kable(
    escape = FALSE,
    caption = "Parent-offspring-pair simulation input sample sizes by capture year."
  )

Table 15: Parent-offspring-pair simulation input sample sizes by capture year.

Data	Capture year	Cohorts	Capture coverage	Total samples
Actual	2017	14	5-30	14,787,909
Actual	2018	15	5-30	11,741,100
Actual	2019	16	1-30	42,342,124
Actual	2020	17	4-30	25,712,890
Actual	2021	18	4-30	12,899,880
Projected	2022	16	6-30	20,303,250
Projected	2023	17	6-30	21,181,500
Projected	2024	18	6-30	21,660,750
Projected	2025	19	6-30	22,244,250
Projected	2026	20	6-30	22,765,500
Projected	2027	21	6-30	22,917,000
Projected	2028	22	6-30	23,017,500
Projected	2029	23	6-30	23,244,750
Projected	2030	24	6-30	23,577,000
Projected	2031	25	6-30	23,984,250
Projected	2032	25	6-30	24,447,000
Projected	2033	25	6-30	24,874,500
Projected	2034	25	6-30	25,251,750
Projected	2035	25	6-30	25,521,750

hist_pop_pred_summary <- hist_pop_pred_df |>
  group_by(year) |>
  summarise(projection_quantiles(pop_expected), .groups = "drop") |>
  mutate(series = "Predicted")

proj_pop_summary <- proj_pop_df |>
  pivot_longer(c(pop_expected, pop_matches), names_to = "series", values_to = "value") |>
  mutate(
    series = recode(
      series,
      pop_expected = "Predicted",
      pop_matches = "Projected simulated observation"
    )
  ) |>
  group_by(series, year) |>
  summarise(projection_quantiles(value), .groups = "drop")

pop_summary_plot <- bind_rows(hist_pop_pred_summary, proj_pop_summary)
pop_pred_plot <- pop_summary_plot |> filter(series == "Predicted")
pop_sim_plot <- pop_summary_plot |> filter(series == "Projected simulated observation")

ggplot() +
  geom_ribbon(
    data = pop_pred_plot,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#0072B2",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = pop_pred_plot,
    aes(x = year, y = median, color = series),
    linewidth = 0.7
  ) +
  geom_linerange(
    data = pop_sim_plot,
    aes(x = year, ymin = lower, ymax = upper, color = series),
    alpha = 0.5,
    linewidth = 0.45
  ) +
  geom_point(
    data = pop_sim_plot,
    aes(x = year, y = median, color = series),
    size = 1.5
  ) +
  geom_point(
    data = hist_pop_df,
    aes(x = year, y = pop_matches, color = "Observed"),
    inherit.aes = FALSE,
    size = 1.6,
    alpha = 0.8
  ) +
  scale_color_manual(
    values = c(
      "Observed" = "black",
      "Predicted" = "#0072B2",
      "Projected simulated observation" = "#D55E00"
    )
  ) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Number of matches", color = NULL)

Figure 11: Historical and projected parent-offspring-pair matches. Black points are observed historical POP matches; blue lines and ribbons are expected matches from the model and projection through 2035; orange points and ranges are annual projected simulated observations using the old projection POP sampling design from sbtproj.tpl.

Catch and TAC

hist_catch_df <- as_tibble(
  as.data.frame.table(base_fit$data$catch_obs_ysf, responseName = "catch")
) |>
  transmute(
    year = as.integer(as.character(Year)),
    season = factor(Season),
    fishery = factor(fishery_names[as.integer(Fishery)], levels = fishery_names),
    catch = catch
  )

proj_catch_output_df <- map_dfr(seq_along(proj_dyn), function(i) {
  as_tibble(as.data.frame.table(proj_dyn[[i]]$catch_pred_ysf, responseName = "catch")) |>
    filter(as.integer(Var1) <= length(projection_dynamics_years)) |>
    transmute(
      iteration = i,
      year = projection_dynamics_years[as.integer(Var1)],
      season = factor(dimnames(ctp_catch_iysf)$Season[as.integer(Var2)], levels = dimnames(ctp_catch_iysf)$Season),
      fishery = factor(dimnames(ctp_catch_iysf)$Fishery[as.integer(Var3)], levels = dimnames(ctp_catch_iysf)$Fishery),
      catch = catch
    )
})

hist_catch_summary <- hist_catch_df |>
  group_by(year, fishery) |>
  summarise(catch = sum(catch, na.rm = TRUE), .groups = "drop")

catch_check_seasons <- dimnames(ctp_catch_iysf)$Season
catch_check_fisheries <- dimnames(ctp_catch_iysf)$Fishery

proj_catch_input_full <- as_tibble(
  as.data.frame.table(ctp_catch_iysf, responseName = "catch_input")
) |>
  transmute(
    iteration = as.integer(as.character(iteration)),
    year = as.integer(as.character(Year)),
    season = factor(as.character(Season), levels = catch_check_seasons),
    fishery = factor(as.character(Fishery), levels = catch_check_fisheries),
    catch_input = catch_input
  )

proj_catch_input_summary <- proj_catch_input_full |>
  group_by(iteration, year, fishery) |>
  summarise(catch = sum(catch_input, na.rm = TRUE), .groups = "drop") |>
  group_by(year, fishery) |>
  summarise(projection_quantiles(catch), .groups = "drop")

catch_output_comparison <- proj_catch_output_df |>
  left_join(proj_catch_input_full, by = join_by(iteration, year, season, fishery)) |>
  mutate(abs_diff = abs(catch - catch_input))

if (any(is.na(catch_output_comparison$catch_input))) {
  stop("Projection catch diagnostic could not match all output catch rows to input catch rows.", call. = FALSE)
}

max_projection_catch_abs_diff <- max(catch_output_comparison$abs_diff, na.rm = TRUE)
if (is.finite(max_projection_catch_abs_diff) && max_projection_catch_abs_diff > 1e-5) {
  stop(
    "Projection output catch differs from input catch. Maximum absolute difference: ",
    signif(max_projection_catch_abs_diff, 8),
    call. = FALSE
  )
}

proj_catch_summary <- catch_output_comparison |>
  group_by(iteration, year, fishery) |>
  summarise(catch = sum(catch, na.rm = TRUE), .groups = "drop") |>
  group_by(year, fishery) |>
  summarise(projection_quantiles(catch), .groups = "drop")

proj_total_tac_summary <- imap_dfr(
  projection_scenario_runs,
  function(run, scenario) {
    as_tibble(
      as.data.frame.table(run$ctp_total_tac_iy, responseName = "tac")
    ) |>
      transmute(
        iteration = as.integer(as.character(iteration)),
        year = as.integer(as.character(Year)),
        tac = tac
      ) |>
      filter(year %in% run$projection_future_years) |>
      group_by(year) |>
      summarise(projection_quantiles(tac), .groups = "drop") |>
      mutate(
        scenario = scenario,
        series = if_else(
          year %in% run$fixed_catch_years,
          "Fixed total TAC",
          "CTP total TAC"
        )
      )
  }
)

recent_tac_context <- hist_catch_df |>
  filter(.data$year %in% tail(projection_actual_catch_years, 5L)) |>
  group_by(year) |>
  summarise(median = sum(catch, na.rm = TRUE), .groups = "drop") |>
  mutate(
    lower = median,
    upper = median,
    series = "Observed context"
  )

ggplot(proj_total_tac_summary, aes(x = year, y = median)) +
  geom_ribbon(
    data = filter(proj_total_tac_summary, .data$series == "CTP total TAC"),
    aes(ymin = lower, ymax = upper),
    fill = "#D55E00",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = recent_tac_context,
    aes(x = year, y = median, color = series, group = 1),
    inherit.aes = FALSE,
    linewidth = 0.65
  ) +
  geom_line(
    data = filter(
      proj_total_tac_summary,
      .data$series == "Fixed total TAC"
    ),
    aes(color = series, group = scenario),
    linewidth = 0.75
  ) +
  geom_line(
    data = filter(proj_total_tac_summary, .data$series == "CTP total TAC"),
    aes(color = series),
    linewidth = 0.75
  ) +
  geom_point(
    data = proj_total_tac_summary,
    aes(color = series),
    size = 1.7
  ) +
  geom_point(
    data = recent_tac_context,
    aes(x = year, y = median, color = series),
    inherit.aes = FALSE,
    size = 1.8
  ) +
  facet_wrap(vars(scenario), ncol = 1) +
  scale_color_manual(values = c(
    "Observed context" = "black",
    "Fixed total TAC" = "#6A3D9A",
    "CTP total TAC" = "#D55E00"
  )) +
  scale_y_continuous(labels = comma, expand = expansion(mult = c(0.05, 0.08))) +
  scale_x_continuous(
    breaks = c(tail(projection_actual_catch_years, 5L), projection_future_years),
    minor_breaks = NULL
  ) +
  labs(x = "Year", y = "Total TAC", color = NULL)

Figure 12: Total TAC by allocation scenario. A black line connects recent conditioned catch totals, and a purple line connects the fixed total TAC for 2026-2029. The fixed total is identical in both scenarios; only its fleet allocation differs. Orange lines and 95% ribbons show scenario-specific CTP feedback from 2030 onward. Points mark every annual value or projection median.

The black segment is the sum of conditioned fishery catches in the five years before the future TAC period. The purple segment is prescribed and therefore has no uncertainty ribbon. From 2030, the orange median and interval summarize the CTP response across projection draws. The two orange trajectories may diverge because the alternative 2027–2029 fleet allocations generate different biological removals and hence different monitoring feedback, even though their fixed total TAC is identical.

ctp_status_prob <- ctp_status |>
  mutate(
    implementation_year = factor(implementation_year),
    direction = case_when(
      delta_tac > 0 ~ "Increase",
      delta_tac < 0 ~ "Decrease",
      TRUE ~ "No change"
    ),
    direction = factor(direction, levels = c("Decrease", "No change", "Increase"))
  ) |>
  count(implementation_year, direction, name = "draws") |>
  complete(implementation_year, direction, fill = list(draws = 0L)) |>
  group_by(implementation_year) |>
  mutate(probability = draws / sum(draws)) |>
  ungroup()

ggplot(ctp_status_prob, aes(x = "", y = probability, fill = direction)) +
  geom_col(width = 1, color = "white", linewidth = 0.35) +
  geom_text(
    data = filter(ctp_status_prob, probability >= 0.04),
    aes(label = percent(probability, accuracy = 1)),
    position = position_stack(vjust = 0.5),
    color = "white",
    size = 3
  ) +
  coord_polar(theta = "y") +
  facet_wrap(vars(implementation_year), nrow = 1) +
  scale_fill_manual(
    values = c("Decrease" = "#0072B2", "No change" = "grey55", "Increase" = "#D55E00")
  ) +
  labs(fill = NULL) +
  theme_void() +
  theme(
    legend.position = "bottom",
    strip.text = element_text(face = "bold", margin = margin(b = 8))
  )

Figure 13: CTP management-procedure TAC adjustment probabilities by implementation year. Pie slices show the probability across projection draws that the CTP TAC decreases, does not change, or increases.

ggplot() +
  geom_ribbon(
    data = proj_catch_input_summary,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#D55E00",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = proj_catch_input_summary,
    aes(x = year, y = median),
    color = "#D55E00",
    linewidth = 0.8,
    alpha = 0.75
  ) +
  geom_point(
    data = proj_catch_input_summary,
    aes(x = year, y = median),
    color = "#D55E00",
    size = 1.15,
    alpha = 0.8
  ) +
  geom_line(
    data = hist_catch_summary,
    aes(x = year, y = catch),
    color = "black",
    linewidth = 0.4,
    alpha = 0.65
  ) +
  geom_point(
    data = hist_catch_summary,
    aes(x = year, y = catch),
    color = "black",
    size = 1,
    alpha = 0.75
  ) +
  facet_wrap(vars(fishery), scales = "free_y") +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "Biological removals")

Figure 14: Conditioned historical removals and configured future biological removals by fishery. Points mark every historical value and annual projection median. Future values include the base projection UAM multipliers. Dynamics output is checked against the removal input and the report stops if they differ beyond numerical tolerance.

catch_zoom_start_year <- 2015L

hist_catch_zoom <- hist_catch_summary |>
  filter(year >= catch_zoom_start_year) |>
  mutate(series = "Observed context")

proj_catch_input_zoom <- proj_catch_input_summary |>
  filter(
    year >= catch_zoom_start_year,
    year %in% projection_future_years
  ) |>
  mutate(series = "Configured biological removals")

proj_catch_fixed_zoom <- proj_catch_input_zoom |>
  filter(year %in% fixed_catch_years) |>
  mutate(series = "Fixed nominal-TAC years")

proj_catch_output_zoom <- proj_catch_summary |>
  filter(
    year >= catch_zoom_start_year,
    year %in% projection_future_years
  ) |>
  mutate(series = "Dynamics output removals")

ggplot() +
  geom_line(
    data = hist_catch_zoom,
    aes(x = year, y = catch, color = series),
    linewidth = 0.6,
    alpha = 0.8
  ) +
  geom_point(
    data = hist_catch_zoom,
    aes(x = year, y = catch, color = series),
    size = 1.5,
    alpha = 0.8
  ) +
  geom_ribbon(
    data = proj_catch_input_zoom,
    aes(x = year, ymin = lower, ymax = upper),
    fill = "#D55E00",
    alpha = 0.15,
    color = NA
  ) +
  geom_line(
    data = proj_catch_input_zoom,
    aes(x = year, y = median, color = series),
    linewidth = 0.7
  ) +
  geom_point(
    data = proj_catch_input_zoom,
    aes(x = year, y = median, color = series),
    size = 1.25,
    alpha = 0.8
  ) +
  geom_line(
    data = proj_catch_fixed_zoom,
    aes(x = year, y = median, color = series),
    linewidth = 0.8,
    alpha = 0.9
  ) +
  geom_point(
    data = proj_catch_fixed_zoom,
    aes(x = year, y = median, color = series),
    size = 1.6,
    alpha = 0.85
  ) +
  geom_line(
    data = proj_catch_output_zoom,
    aes(x = year, y = median, color = series),
    linewidth = 0.45,
    alpha = 0.9
  ) +
  scale_color_manual(
    values = c(
      "Observed context" = "black",
      "Configured biological removals" = "#D55E00",
      "Fixed nominal-TAC years" = "#6A3D9A",
      "Dynamics output removals" = "#0072B2"
    )
  ) +
  facet_wrap(vars(fishery), scales = "free_y") +
  scale_y_zero(labels = comma) +
  scale_x_continuous(
    breaks = seq(catch_zoom_start_year + 1L, projection_last_yr, by = 4L),
    expand = expansion(mult = c(0.02, 0.03))
  ) +
  labs(x = "Year", y = "Biological removals", color = NULL)

Figure 15: Zoomed biological-removal inputs and checked dynamics output from 2015 onward. Black points and lines show conditioned removals through 2025. Purple points and lines show the fixed nominal-TAC years 2026-2029 after applying fishery-specific removal multipliers. Orange points, lines, and 95% ribbons show CTP-controlled removals from 2030 onward; the underlying configured-removal series begins in 2026. Blue lines are a diagnostic dynamics-output overlay that must match the configured removals.

The 2022–2025 orange points in the previous version were not future TAC decisions. They appeared because the projection dynamics array begins in 2022 and carries the already conditioned catch inputs through 2025. Those years are now shown once, in black, as historical conditioned removals. Purple then marks the fixed 2026–2029 period, while orange carries the configured future-removal trajectory into the CTP-controlled years. The blue overlay remains visible only as a numerical accounting diagnostic.

Risk

The risk summary is a projection diagnostic for the CTP period. The short-term window is the first three-year CTP block (2030-2032); the long-term window contains the subsequent CTP and final state years (2033-2036). For each window the table separates two quantities: the mean of the annual probabilities across draws, and the probability that a draw crosses the threshold at least once. Relative TRO is calculated as spawning_biomass / B0, and the configured threshold is 0.2. This follows the risk-style interpretation used in CCSBT assessment and management-procedure evaluations, where projected TRO outcomes are compared against relative TRO reference points rather than fitted as an additional data stream (Hillary et al. 2023; Commission for the Conservation of Southern Bluefin Tuna 2026).

risk_windows <- bind_rows(
  tibble(window = "Short term", year = short_term_risk_years),
  tibble(window = "Long term", year = long_term_risk_years)
) |>
  mutate(window = factor(window, levels = c("Short term", "Long term")))

risk_data <- proj_spawning_df |>
  inner_join(risk_windows, by = "year") |>
  mutate(below_threshold = relative_spawning_biomass < risk_relative_biomass_threshold)

risk_annual <- risk_data |>
  group_by(window, year) |>
  summarise(annual_probability = mean(below_threshold, na.rm = TRUE), .groups = "drop") |>
  group_by(window) |>
  summarise(
    years = paste(range(year), collapse = "-"),
    `Mean annual probability below threshold` = mean(annual_probability),
    .groups = "drop"
  )

risk_ever <- risk_data |>
  group_by(window, iteration) |>
  summarise(crossed_at_least_once = any(below_threshold, na.rm = TRUE), .groups = "drop") |>
  group_by(window) |>
  summarise(
    `Probability a draw crosses at least once` = mean(crossed_at_least_once),
    .groups = "drop"
  )

risk_summary <- risk_annual |>
  left_join(risk_ever, by = "window") |>
  mutate(threshold = risk_relative_biomass_threshold)

risk_summary |>
  mutate(across(where(is.numeric), ~ formatC(.x, format = "f", digits = 3))) |>
  kable(
    caption = "CTP-period risk summary showing mean annual and per-draw ever-crossing probabilities."
  )

Table 16: CTP-period risk summary showing mean annual and per-draw ever-crossing probabilities.

window	years	Mean annual probability below threshold	Probability a draw crosses at least once	threshold
Short term	2030-2032	0.001	0.001	0.200
Long term	2033-2036	0.009	0.018	0.200

Assessment Projection Summary

The following table is the current counterpart of the projection summary in the 2023 assessment. It contains the accepted 2,000-draw outputs for all four production arms. The files retain projection cache format 8. The format-9 change stages not-yet-decided future catches as zero during intermediate CTP updates and records realized catch separately. All four completed format-8 paths had zero harvest penalties and maximum raw harvest below 0.57, so the fail-closed runs had no catch shortfall. Because later placeholder catches cannot affect observations available to an earlier CTP update, the stored results remain scientifically valid and are read through exact signature and payload fingerprints; they are not relabelled as format 9. The current reporting contract projects dynamics through 2035 and reports the resulting state through 2036, so a 2036 state column replaces the former 2040 column. Mean TAC is the mean annual nominal TAC over 2026–2035 for each draw, summarized across draws. This assessment-format table deliberately reports a 2035 point-in-time stock-status probability; it is distinct from the mean-annual and per-draw ever-crossing window diagnostics used in the risk comparison tables.

projection_summary_file <- file.path(
  esc_dir,
  "assessment_doc", "data", "table6_projection_summary.csv"
)
if (!file.exists(projection_summary_file)) {
  stop(
    "The assessment summary tables are missing. Run ",
    "`Rscript scripts/build-esc31-assessment-tables.R .`.",
    call. = FALSE
  )
}
assessment_projection_summary <- readr::read_csv(
  projection_summary_file,
  show_col_types = FALSE
)

projection_table_interval <- function(median, lower, upper, digits = 3L) {
  paste0(
    format_decimal(median, digits),
    " (", format_decimal(lower, digits),
    "-", format_decimal(upper, digits), ")"
  )
}

assessment_projection_summary |>
  transmute(
    Scenario = .data$Scenario,
    `P(TRO_2035 > 0.2 TRO_0)` = scales::percent(
      .data$Probability_relative_TRO_2035_above_0_2,
      accuracy = 0.1
    ),
    `P(TRO_2035 > 0.3 TRO_0)` = scales::percent(
      .data$Probability_relative_TRO_2035_above_0_3,
      accuracy = 0.1
    ),
    `TRO_2025/TRO_0` = projection_table_interval(
      .data$Relative_TRO_2025_median,
      .data$Relative_TRO_2025_lower,
      .data$Relative_TRO_2025_upper
    ),
    `TRO_2035/TRO_0` = projection_table_interval(
      .data$Relative_TRO_2035_median,
      .data$Relative_TRO_2035_lower,
      .data$Relative_TRO_2035_upper
    ),
    `TRO_2036/TRO_0` = projection_table_interval(
      .data$Relative_TRO_2036_median,
      .data$Relative_TRO_2036_lower,
      .data$Relative_TRO_2036_upper
    ),
    `Mean nominal TAC, 2026-2035 (t)` = projection_table_interval(
      .data$Mean_2026_2035_nominal_TAC_t_median,
      .data$Mean_2026_2035_nominal_TAC_t_lower,
      .data$Mean_2026_2035_nominal_TAC_t_upper,
      0L
    )
  ) |>
  kable(align = c("l", rep("r", 6L)))

Table 17: Projection summary from the accepted format-8 production runs. Exact cache fingerprints, complete CTP-input audits, zero harvest penalties, and maximum raw harvest below 0.57 support reporting compatibility with the format-9 catch-staging change. Probabilities refer to the 2035 state. State and TAC entries are medians and equal-tailed 95% intervals across 2,000 draws. The sampled direct-M MLE-grid intervals describe weighted grid-resampling spread and are not posterior credible intervals.

Scenario	P(TRO_2035 > 0.2 TRO_0)	P(TRO_2035 > 0.3 TRO_0)	TRO_2025/TRO_0	TRO_2035/TRO_0	TRO_2036/TRO_0	Mean nominal TAC, 2026-2035 (t)
Balanced MCMC grid — nominal fleet allocations	99.0%	71.2%	0.271 (0.201-0.359)	0.339 (0.218-0.542)	0.338 (0.209-0.563)	25,399 (25,349-26,249)
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	99.2%	71.6%	0.271 (0.201-0.359)	0.341 (0.220-0.544)	0.341 (0.211-0.566)	25,400 (25,349-26,249)
NoUAM — nominal fleet allocations (NCNM removed)	100.0%	86.9%	0.279 (0.223-0.348)	0.364 (0.259-0.541)	0.364 (0.251-0.560)	25,399 (25,349-26,249)
Sampled direct-M MLE grid — nominal fleet allocations	100.0%	73.0%	0.268 (0.219-0.316)	0.337 (0.245-0.490)	0.332 (0.236-0.505)	25,392 (25,349-26,249)

Scenarios

Allocation Scenario Comparison

All detailed plots and tables above use the nominal fleet allocations arm. The 3,000 t TAC increase assigned to Indonesia arm uses the same posterior draws, recruitment, selectivity, UAM factors, observation simulation settings, CTP schedule, and fixed total TAC. The only intended difference is the pre-CTP nominal fleet allocation in 2027-2029 and the CTP fleet-allocation proportions derived from those non-carryover pre-CTP years. Total TAC is identical in both arms for every fixed year, 2026-2029.

These two scenarios implement the OMMP16 Alloc full-grid sensitivity. The scenario definitions are complete. The current canonical render shows the 2,000-draw production comparison. All report gates and the visual review for this production comparison have passed. The 108-fit direct-get_M MLE grid is the separate Og comparison rather than an input to these projections.

fixed_projection_tac_scenarios |>
  transmute(
    Scenario = unname(projection_scenario_labels[.data$scenario]),
    Year = .data$year,
    LL1 = format_count(.data$LL1),
    LL2 = format_count(.data$LL2),
    Indonesia = format_count(.data$Indonesia),
    Australia = format_count(.data$Australia),
    Total = format_count(.data$total)
  ) |>
  kable()

Table 18: Fixed pre-CTP nominal TAC allocations. Total TAC is identical in both allocation scenarios.

Scenario	Year	LL1	LL2	Indonesia	Australia	Total
Balanced MCMC grid — nominal fleet allocations	2026	14,326	1,534	1,248	5,562	22,671
Balanced MCMC grid — nominal fleet allocations	2027	15,063	1,653	1,370	5,561	23,647
Balanced MCMC grid — nominal fleet allocations	2028	15,063	1,653	1,370	5,561	23,647
Balanced MCMC grid — nominal fleet allocations	2029	15,063	1,653	1,370	5,561	23,647
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2026	14,326	1,534	1,248	5,562	22,671
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2027	13,151	1,443	4,197	4,856	23,647
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2028	13,151	1,443	4,197	4,856	23,647
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2029	13,151	1,443	4,197	4,856	23,647

This render uses 2,000 balanced posterior draws. Its recorded elapsed times are reported directly; any runtime rescaled from a different cache size is explicitly labelled as a linear estimate: Balanced MCMC grid — nominal fleet allocations: 42.4 minutes (observed); Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia: 42.4 minutes (observed); both scenarios sequentially: 84.7 minutes.

imap_dfr(projection_scenario_runs, projection_timing_row) |>
  kable(escape = FALSE)

Table 19: Projection run size and elapsed runtime by projection scenario.

Scenario	Projection iterations	Run elapsed	Minutes per iteration	Loaded from cache	Cache load elapsed
Balanced MCMC grid — nominal fleet allocations	2,000	42.4 minutes	0.021	Yes	–
Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	2,000	42.4 minutes	0.021	Yes	–

scenario_biomass_comparison_start_year <- 2010L

projection_scenario_spawning_summary <- imap_dfr(projection_scenario_runs, function(run, scenario) {
  state_years <- run$projection_dynamics_first_yr:(max(run$projection_dynamics_years) + 1L)

  map_dfr(seq_along(run$proj_dyn), function(i) {
    tibble(
      iteration = i,
      year = state_years,
      spawning_biomass = as.numeric(run$proj_dyn[[i]]$spawning_biomass_y)
    )
  }) |>
    group_by(year) |>
    summarise(projection_quantiles(spawning_biomass), .groups = "drop") |>
    mutate(scenario = scenario)
})

ggplot() +
  geom_line(
    data = filter(hist_spawning_summary, year >= scenario_biomass_comparison_start_year),
    aes(x = year, y = median),
    color = "grey25",
    linewidth = 0.55,
    alpha = 0.8
  ) +
  geom_ribbon(
    data = filter(projection_scenario_spawning_summary, year >= scenario_biomass_comparison_start_year),
    aes(x = year, ymin = lower, ymax = upper, fill = scenario),
    alpha = 0.14,
    color = NA
  ) +
  geom_line(
    data = filter(projection_scenario_spawning_summary, year >= scenario_biomass_comparison_start_year),
    aes(x = year, y = median, color = scenario),
    linewidth = 0.75
  ) +
  scale_color_manual(values = projection_scenario_colors) +
  scale_fill_manual(values = projection_scenario_colors) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "TRO", color = NULL, fill = NULL)

Figure 16: Projected total reproductive output for the nominal fleet allocations and the 3,000 t TAC increase assigned to Indonesia. Lines show posterior medians and shaded ribbons show 95% intervals; the historical median is shown for context.

projection_scenario_catch_summary <- imap_dfr(projection_scenario_runs, function(run, scenario) {
  as_tibble(
    as.data.frame.table(
      run$ctp_nominal_catch_iysf,
      responseName = "nominal_tac"
    )
  ) |>
    transmute(
      iteration = as.integer(as.character(iteration)),
      year = as.integer(as.character(Year)),
      fishery = factor(
        as.character(Fishery),
        levels = dimnames(run$ctp_nominal_catch_iysf)$Fishery
      ),
      nominal_tac = nominal_tac
    ) |>
    group_by(iteration, year, fishery) |>
    summarise(
      nominal_tac = sum(nominal_tac, na.rm = TRUE),
      .groups = "drop"
    ) |>
    group_by(year, fishery) |>
    summarise(projection_quantiles(nominal_tac), .groups = "drop") |>
    mutate(scenario = scenario)
}) |>
  filter(year %in% projection_future_years)

ggplot() +
  geom_ribbon(
    data = projection_scenario_catch_summary,
    aes(x = year, ymin = lower, ymax = upper, fill = scenario),
    alpha = 0.12,
    color = NA
  ) +
  geom_line(
    data = projection_scenario_catch_summary,
    aes(x = year, y = median, color = scenario),
    linewidth = 0.7
  ) +
  geom_point(
    data = projection_scenario_catch_summary,
    aes(x = year, y = median, color = scenario),
    size = 1.4,
    alpha = 0.85
  ) +
  facet_wrap(vars(fishery), scales = "free_y") +
  scale_color_manual(values = projection_scenario_colors) +
  scale_fill_manual(values = projection_scenario_colors) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(
    breaks = seq(
      min(projection_future_years),
      max(projection_future_years),
      by = 2L
    ),
    minor_breaks = projection_future_years,
    expand = expansion(mult = c(0.02, 0.03))
  ) +
  labs(
    x = "Year",
    y = "Nominal TAC allocation",
    color = NULL,
    fill = NULL
  ) +
  theme(legend.position = "bottom")

Figure 17: Nominal TAC allocations by fishery for the nominal fleet allocations and the 3,000 t TAC increase assigned to Indonesia. Points mark every annual median. Fixed total TAC is identical in 2026-2029; the panels show how that same total is allocated differently.

The fixed total TAC is equal on the nominal reporting scale. Because the LL1 and Australian removal multipliers exceed one, reallocating that TAC among fleets can nevertheless change total biological removals, as shown below.

projection_scenario_fixed_removals <- imap_dfr(
  projection_scenario_runs,
  function(run, scenario) {
    as_tibble(
      as.data.frame.table(
        run$ctp_total_removal_iy,
        responseName = "removals"
      )
    ) |>
      transmute(
        scenario = scenario,
        year = as.integer(as.character(Year)),
        removals = removals
      ) |>
      filter(year %in% run$fixed_catch_years) |>
      group_by(scenario, year) |>
      summarise(removals = median(removals), .groups = "drop")
  }
)
projection_scenario_fixed_removals |>
  pivot_wider(
    names_from = scenario,
    values_from = removals
  ) |>
  mutate(
    removal_difference =
      .data[[projection_scenario_labels[["scenario_2"]]]] -
        .data[[projection_scenario_labels[["scenario_1"]]]]
  ) |>
  transmute(
    Year = .data$year,
    `Identical nominal total TAC` = format_count(
      fixed_projection_tac[match(.data$year, fixed_catch_years)]
    ),
    !!projection_arm_labels[["base"]] := format_count(
      .data[[projection_scenario_labels[["scenario_1"]]]]
    ),
    !!projection_arm_labels[["indonesia_allocation"]] := format_count(
      .data[[projection_scenario_labels[["scenario_2"]]]]
    ),
    `Indonesia allocation minus nominal` =
      format_count(.data$removal_difference)
  ) |>
  kable()

Table 20: Fixed-period biological removals under the two allocation scenarios. Nominal total TAC is identical; removal totals differ only because the same TAC is allocated among fisheries with different removal multipliers.

Year	Identical nominal total TAC	Balanced MCMC grid — nominal fleet allocations	Balanced MCMC grid — 3,000 t TAC increase assigned to Indonesia	Indonesia allocation minus nominal
2026	22,671	25,359	25,359	0
2027	23,647	26,416	26,065	-351
2028	23,647	26,416	26,065	-351
2029	23,647	26,416	26,065	-351

projection_scenario_tac_summary <- imap_dfr(projection_scenario_runs, function(run, scenario) {
  as_tibble(as.data.frame.table(run$ctp_total_tac_iy, responseName = "tac")) |>
    transmute(
      scenario = scenario,
      iteration = as.integer(as.character(iteration)),
      year = as.integer(as.character(Year)),
      tac = tac,
      source = if_else(year %in% run$fixed_catch_years, "Pre-CTP", "CTP")
    ) |>
    filter(year %in% run$projection_future_years) |>
    group_by(scenario, source, year) |>
    summarise(projection_quantiles(tac), .groups = "drop")
})

projection_scenario_tac_wide <- projection_scenario_tac_summary |>
  select(scenario, source, year, median, lower, upper) |>
  pivot_wider(
    names_from = scenario,
    values_from = c(median, lower, upper),
    names_glue = "{scenario}_{.value}"
  )
nominal_label <- projection_scenario_labels[["scenario_1"]]
indonesia_label <- projection_scenario_labels[["scenario_2"]]

projection_scenario_tac_wide |>
  transmute(
    Source = source,
    Year = year,
    `Nominal-allocation median TAC` =
      format_count(.data[[paste0(nominal_label, "_median")]]),
    `Nominal-allocation 95% interval` = paste0(
      format_count(.data[[paste0(nominal_label, "_lower")]]),
      "-",
      format_count(.data[[paste0(nominal_label, "_upper")]])
    ),
    `Indonesia-allocation median TAC` =
      format_count(.data[[paste0(indonesia_label, "_median")]]),
    `Indonesia-allocation 95% interval` = paste0(
      format_count(.data[[paste0(indonesia_label, "_lower")]]),
      "-",
      format_count(.data[[paste0(indonesia_label, "_upper")]])
    ),
    `Indonesia minus nominal median` = format_count(
      .data[[paste0(indonesia_label, "_median")]] -
        .data[[paste0(nominal_label, "_median")]]
    )
  ) |>
  arrange(Year) |>
  kable()

Table 21: Total TAC comparison between the two meaning-labelled allocation scenarios.

Source	Year	Nominal-allocation median TAC	Nominal-allocation 95% interval	Indonesia-allocation median TAC	Indonesia-allocation 95% interval	Indonesia minus nominal median
Pre-CTP	2026	22,671	22,671-22,671	22,671	22,671-22,671	0
Pre-CTP	2027	23,647	23,647-23,647	23,647	23,647-23,647	0
Pre-CTP	2028	23,647	23,647-23,647	23,647	23,647-23,647	0
Pre-CTP	2029	23,647	23,647-23,647	23,647	23,647-23,647	0
CTP	2030	26,647	26,647-26,647	26,647	26,647-26,647	0
CTP	2031	26,647	26,647-26,647	26,647	26,647-26,647	0
CTP	2032	26,647	26,647-26,647	26,647	26,647-26,647	0
CTP	2033	26,814	26,647-29,647	26,816	26,647-29,647	2
CTP	2034	26,814	26,647-29,647	26,816	26,647-29,647	2
CTP	2035	26,814	26,647-29,647	26,816	26,647-29,647	2

NoUAM Projection Comparison

This comparison carries the OMMP16 NoUAM instruction through both stages. The base arm is conditioned with the recorded NCNM additions and applies the LL1 1.11 factor to future nominal TAC allocations. The NoUAM arm removes the NCNM catch additions in conditioning and uses an LL1 factor of 1.00 in projection. Both arms retain the separate Australian surface-fishery factor of 1.20 and use the nominal fleet-allocation scenario, the same fixed total TAC, monitoring programme, CTP lags, projection years, and projection seed.

bind_rows(
  projection_timing_row(
    projection_run,
    projection_arm_labels[["base"]]
  ),
  projection_timing_row(
    no_uam_projection_run,
    projection_arm_labels[["no_uam"]]
  )
) |>
  kable(escape = FALSE)

Table 22: Projection run size and elapsed runtime for the base-with-NCNM and NoUAM comparison arms.

Scenario	Projection iterations	Run elapsed	Minutes per iteration	Loaded from cache	Cache load elapsed
Balanced MCMC grid — nominal fleet allocations	2,000	42.4 minutes	0.021	Yes	–
NoUAM — nominal fleet allocations (NCNM removed)	2,000	42.9 minutes	0.021	Yes	–

no_uam_projection_comparison <- setNames(
  list(
    list(
      run = projection_run,
      post = projection_post
    ),
    list(
      run = no_uam_projection_run,
      post = no_uam_projection_post
    )
  ),
  unname(projection_arm_labels[c("base", "no_uam")])
)

no_uam_projection_spawning_draws <- imap_dfr(
  no_uam_projection_comparison,
  function(arm, model) {
    run <- arm$run
    post <- arm$post
    if (length(run$proj_dyn) != nrow(post)) {
      stop(
        "Projection dynamics and posterior rows differ for ", model, ".",
        call. = FALSE
      )
    }
    state_years <- run$projection_dynamics_first_yr:
      (max(run$projection_dynamics_years) + 1L)
    map_dfr(seq_along(run$proj_dyn), function(i) {
      b0 <- exp(post$par_log_B0[i])
      tibble(
        model = model,
        iteration = i,
        year = state_years,
        spawning_biomass =
          as.numeric(run$proj_dyn[[i]]$spawning_biomass_y),
        relative_spawning_biomass = spawning_biomass / b0
      )
    })
  }
)
no_uam_projection_spawning_summary <-
  no_uam_projection_spawning_draws |>
  group_by(.data$model, .data$year) |>
  summarise(
    projection_quantiles(.data$spawning_biomass),
    .groups = "drop"
  )
no_uam_projection_colors <- setNames(
  c("#D55E00", "#0072B2"),
  unname(projection_arm_labels[c("base", "no_uam")])
)

ggplot(
  filter(
    no_uam_projection_spawning_summary,
    year >= projection_first_yr
  ),
  aes(x = year, y = median, color = model, fill = model)
) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.14, color = NA) +
  geom_line(linewidth = 0.8) +
  scale_color_manual(
    values = no_uam_projection_colors
  ) +
  scale_fill_manual(
    values = no_uam_projection_colors
  ) +
  scale_y_zero(labels = comma) +
  scale_x_continuous(breaks = pretty_breaks()) +
  labs(x = "Year", y = "TRO", color = NULL, fill = NULL) +
  theme(
    plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 18)
  )

Figure 18: Projected total reproductive output when NCNM is included in both conditioning and projection versus removed from both. Both arms use the nominal fleet-allocation TAC scenario. Lines show posterior medians and ribbons show 95% intervals.

no_uam_projection_risk_data <- no_uam_projection_spawning_draws |>
  inner_join(risk_windows, by = "year") |>
  mutate(
    below_threshold =
      .data$relative_spawning_biomass < risk_relative_biomass_threshold
  )
no_uam_projection_risk_annual <- no_uam_projection_risk_data |>
  group_by(.data$model, .data$window, .data$year) |>
  summarise(
    annual_probability = mean(.data$below_threshold),
    .groups = "drop"
  ) |>
  group_by(.data$model, .data$window) |>
  summarise(
    Years = paste(range(.data$year), collapse = "-"),
    `Mean annual probability below threshold` =
      mean(.data$annual_probability),
    .groups = "drop"
  )
no_uam_projection_risk_ever <- no_uam_projection_risk_data |>
  group_by(.data$model, .data$window, .data$iteration) |>
  summarise(crossed = any(.data$below_threshold), .groups = "drop") |>
  group_by(.data$model, .data$window) |>
  summarise(
    `Probability a draw crosses at least once` = mean(.data$crossed),
    .groups = "drop"
  )
left_join(
  no_uam_projection_risk_annual,
  no_uam_projection_risk_ever,
  by = c("model", "window")
) |>
  transmute(
    Model = .data$model,
    Window = .data$window,
    Years,
    `Relative TRO threshold` =
      percent(risk_relative_biomass_threshold, accuracy = 0.1),
    `Mean annual probability` = percent(
      .data$`Mean annual probability below threshold`,
      accuracy = 0.1
    ),
    `Ever-crossing probability` = percent(
      .data$`Probability a draw crosses at least once`,
      accuracy = 0.1
    )
  ) |>
  kable()

Table 23: Relative-TRO risk comparison with NCNM included in both conditioning and projection versus removed from both stages. It uses the same mean-annual and per-draw ever-crossing definitions as the base and MLE-grid risk tables.

Model	Window	Years	Relative TRO threshold	Mean annual probability	Ever-crossing probability
Balanced MCMC grid — nominal fleet allocations	Short term	2030-2032	20.0%	0.1%	0.1%
Balanced MCMC grid — nominal fleet allocations	Long term	2033-2036	20.0%	0.9%	1.8%
NoUAM — nominal fleet allocations (NCNM removed)	Short term	2030-2032	20.0%	0.0%	0.0%
NoUAM — nominal fleet allocations (NCNM removed)	Long term	2033-2036	20.0%	0.1%	0.2%

imap_dfr(no_uam_projection_comparison, function(arm, model) {
  run <- arm$run
  as_tibble(
    as.data.frame.table(run$ctp_total_tac_iy, responseName = "tac")
  ) |>
    transmute(
      model = model,
      year = as.integer(as.character(Year)),
      tac = tac,
      source = if_else(
        year %in% run$fixed_catch_years,
        "Pre-CTP",
        "CTP"
      )
    ) |>
    filter(year %in% run$projection_future_years) |>
    group_by(.data$model, .data$source, .data$year) |>
    summarise(projection_quantiles(.data$tac), .groups = "drop")
}) |>
  transmute(
    Model = .data$model,
    Source = .data$source,
    Year = .data$year,
    `Median TAC` = format_count(.data$median),
    `95% interval` = paste0(
      format_count(.data$lower),
      "-",
      format_count(.data$upper)
    )
  ) |>
  arrange(.data$Year, .data$Model) |>
  kable()

Table 24: Nominal total TAC with NCNM included versus removed. Fixed TAC is the same; later CTP TAC may differ through population feedback.

Model	Source	Year	Median TAC	95% interval
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2026	22,671	22,671-22,671
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2026	22,671	22,671-22,671
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2027	23,647	23,647-23,647
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2027	23,647	23,647-23,647
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2028	23,647	23,647-23,647
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2028	23,647	23,647-23,647
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2029	23,647	23,647-23,647
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2029	23,647	23,647-23,647
Balanced MCMC grid — nominal fleet allocations	CTP	2030	26,647	26,647-26,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2030	26,647	26,647-26,647
Balanced MCMC grid — nominal fleet allocations	CTP	2031	26,647	26,647-26,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2031	26,647	26,647-26,647
Balanced MCMC grid — nominal fleet allocations	CTP	2032	26,647	26,647-26,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2032	26,647	26,647-26,647
Balanced MCMC grid — nominal fleet allocations	CTP	2033	26,814	26,647-29,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2033	26,812	26,647-29,647
Balanced MCMC grid — nominal fleet allocations	CTP	2034	26,814	26,647-29,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2034	26,812	26,647-29,647
Balanced MCMC grid — nominal fleet allocations	CTP	2035	26,814	26,647-29,647
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2035	26,812	26,647-29,647

imap_dfr(no_uam_projection_comparison, function(arm, model) {
  run <- arm$run
  as_tibble(
    as.data.frame.table(
      run$ctp_total_removal_iy,
      responseName = "removals"
    )
  ) |>
    transmute(
      model = model,
      year = as.integer(as.character(Year)),
      removals = removals,
      source = if_else(
        year %in% run$fixed_catch_years,
        "Pre-CTP",
        "CTP"
      )
    ) |>
    filter(year %in% run$projection_future_years) |>
    group_by(.data$model, .data$source, .data$year) |>
    summarise(
      projection_quantiles(.data$removals),
      .groups = "drop"
    )
}) |>
  transmute(
    Model = .data$model,
    Source = .data$source,
    Year = .data$year,
    `Median biological removals` = format_count(.data$median),
    `95% interval` = paste0(
      format_count(.data$lower),
      "-",
      format_count(.data$upper)
    )
  ) |>
  arrange(.data$Year, .data$Model) |>
  kable()

Table 25: Total biological removals with NCNM included in both conditioning and projection versus removed from both. These totals differ from nominal TAC where a removal multiplier exceeds one.

Model	Source	Year	Median biological removals	95% interval
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2026	25,359	25,359-25,359
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2026	23,783	23,783-23,783
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2027	26,416	26,416-26,416
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2027	24,759	24,759-24,759
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2028	26,416	26,416-26,416
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2028	24,759	24,759-24,759
Balanced MCMC grid — nominal fleet allocations	Pre-CTP	2029	26,416	26,416-26,416
NoUAM — nominal fleet allocations (NCNM removed)	Pre-CTP	2029	24,759	24,759-24,759
Balanced MCMC grid — nominal fleet allocations	CTP	2030	29,767	29,767-29,767
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2030	27,900	27,900-27,900
Balanced MCMC grid — nominal fleet allocations	CTP	2031	29,767	29,767-29,767
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2031	27,900	27,900-27,900
Balanced MCMC grid — nominal fleet allocations	CTP	2032	29,767	29,767-29,767
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2032	27,900	27,900-27,900
Balanced MCMC grid — nominal fleet allocations	CTP	2033	29,954	29,767-33,119
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2033	28,073	27,900-31,041
Balanced MCMC grid — nominal fleet allocations	CTP	2034	29,954	29,767-33,119
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2034	28,073	27,900-31,041
Balanced MCMC grid — nominal fleet allocations	CTP	2035	29,954	29,767-33,119
NoUAM — nominal fleet allocations (NCNM removed)	CTP	2035	28,073	27,900-31,041

Sampled MLE-grid Projection Comparison

This final section projects the full-objective-weighted direct-get_M MLE grid and compares it with the nominal base projection. It is an ESC31 report-specific calculation rather than a new public sbt interface. Both arms use the nominal fleet allocation, identical fixed TAC, monitoring programme, CTP schedule, projection seed, and component-specific stochastic streams.

The uncertainty bands have different meanings. The base arm carries within-cell posterior uncertainty as well as the accepted nine-cell h/psi grid. The MLE-grid arm resamples deterministic fits from the 108-cell direct-M grid using its full-objective weights; its spread represents grid and resampling variation only and is not a posterior credible interval.

mle_grid_projection_source_summary <- tibble(
  Model = unname(projection_arm_labels[c("base", "mle_grid")]),
  `Projection draws` = c(
    nrow(projection_post),
    nrow(mle_grid_projection_post)
  ),
  `Source grid cells represented` = c(
    n_distinct(projection_draw_metadata$cell),
    sum(mle_grid_projection_draw_plan$draws > 0L)
  ),
  `Mortality implementation` = c(
    "Length-scaled get_M_length; M10 and M30 estimated",
    "Direct get_M; M0 and M10 fixed by cell, M4 and M30 estimated"
  ),
  `Historical state audit` = c(
    projection_historical_state$summary$passes,
    mle_grid_projection_historical_state$summary$passes
  )
)
mle_grid_projection_source_summary |>
  mutate(
    `Projection draws` = format_count(.data$`Projection draws`),
    `Historical state audit` = if_else(
      .data$`Historical state audit`,
      "Pass",
      "Fail"
    )
  ) |>
  kable()

Table 26: Source and numerical audit for the base and sampled-MLE-grid projection arms.

Model	Projection draws	Source grid cells represented	Mortality implementation	Historical state audit
Balanced MCMC grid — nominal fleet allocations	2,000	9	Length-scaled get_M_length; M10 and M30 estimated	Pass
Sampled direct-M MLE grid — nominal fleet allocations	2,000	104	Direct get_M; M0 and M10 fixed by cell, M4 and M30 estimated	Pass

bind_rows(
  projection_timing_row(
    projection_run,
    projection_arm_labels[["base"]]
  ),
  projection_timing_row(
    mle_grid_projection_run,
    projection_arm_labels[["mle_grid"]]
  )
) |>
  kable(escape = FALSE)
bind_rows(
  ctp_optimisation_audit_row(
    projection_run,
    projection_arm_labels[["base"]]
  ),
  ctp_optimisation_audit_row(
    mle_grid_projection_run,
    projection_arm_labels[["mle_grid"]]
  )
) |>
  kable(escape = FALSE)

Table 27: Projection runtime and mandatory CTP optimisation audit for the base and sampled-MLE-grid arms.

Scenario	Projection iterations	Run elapsed	Minutes per iteration	Loaded from cache	Cache load elapsed
Balanced MCMC grid — nominal fleet allocations	2,000	42.4 minutes	0.021	Yes	–
Sampled direct-M MLE grid — nominal fleet allocations	2,000	41.7 minutes	0.021	Yes	–

Scenario	Accepted CTP fits	Maximum final |gradient|	Fits requiring refinement	Fits using nlminb fallback	Fits using Newton certification	Maximum deterministic nudge	Minimum Hessian eigenvalue	Maximum Hessian condition number	Positive-definite Hessians
Balanced MCMC grid — nominal fleet allocations	4000	0.000500	345 (8.6%)	84	0	1e-08	12.35	57.2	4000/4000
Sampled direct-M MLE grid — nominal fleet allocations	4000	0.000500	277 (6.9%)	75	0	1e-08	12.285	53.9	4000/4000

mle_grid_projection_comparison <- setNames(
  list(
    list(run = projection_run, post = projection_post),
    list(run = mle_grid_projection_run, post = mle_grid_projection_post)
  ),
  unname(projection_arm_labels[c("base", "mle_grid")])
)
mle_grid_projection_colors <- setNames(
  c("#D55E00", "#0072B2"),
  unname(projection_arm_labels[c("base", "mle_grid")])
)

mle_grid_projection_state_draws <- imap_dfr(
  mle_grid_projection_comparison,
  function(arm, model) {
    run <- arm$run
    post <- arm$post
    if (length(run$proj_dyn) != nrow(post)) {
      stop(
        "Projection dynamics and source rows differ for ",
        model,
        ".",
        call. = FALSE
      )
    }
    state_years <- run$projection_dynamics_first_yr:
      (max(run$projection_dynamics_years) + 1L)
    map_dfr(seq_along(run$proj_dyn), function(i) {
      b0 <- exp(post[i, "par_log_B0"])
      tibble(
        model = model,
        iteration = i,
        year = state_years,
        TRO = as.numeric(run$proj_dyn[[i]]$spawning_biomass_y),
        `Relative TRO` = TRO / b0
      )
    })
  }
)
mle_grid_projection_state_summary <-
  mle_grid_projection_state_draws |>
  pivot_longer(
    cols = c("TRO", "Relative TRO"),
    names_to = "quantity",
    values_to = "value"
  ) |>
  group_by(.data$model, .data$quantity, .data$year) |>
  summarise(projection_quantiles(.data$value), .groups = "drop")

ggplot(
  filter(
    mle_grid_projection_state_summary,
    .data$year >= projection_first_yr
  ),
  aes(
    x = .data$year,
    y = .data$median,
    color = .data$model,
    fill = .data$model
  )
) +
  geom_ribbon(
    aes(ymin = .data$lower, ymax = .data$upper),
    alpha = 0.14,
    color = NA
  ) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 1.35, alpha = 0.85) +
  facet_wrap(vars(.data$quantity), ncol = 1L, scales = "free_y") +
  scale_color_manual(values = mle_grid_projection_colors) +
  scale_fill_manual(values = mle_grid_projection_colors) +
  scale_x_continuous(
    breaks = seq(projection_first_yr, projection_last_yr + 1L, by = 2L),
    minor_breaks = projection_first_yr:(projection_last_yr + 1L)
  ) +
  scale_y_continuous(
    labels = label_comma(),
    expand = expansion(mult = c(0.02, 0.05))
  ) +
  labs(x = "Year", y = NULL, color = NULL, fill = NULL) +
  theme(legend.position = "bottom")

Figure 19: Projected total reproductive output and relative TRO for the balanced MCMC-grid base and the full-objective-weighted direct-M MLE grid. Lines show medians, ribbons show 95% intervals across the respective sampled rows, and points mark every projected year. MLE-grid intervals are not posterior credible intervals.

mle_grid_projection_risk_data <-
  mle_grid_projection_state_draws |>
  inner_join(risk_windows, by = "year") |>
  mutate(
    below_threshold =
      .data$`Relative TRO` < risk_relative_biomass_threshold
  )
mle_grid_projection_risk_annual <-
  mle_grid_projection_risk_data |>
  group_by(.data$model, .data$window, .data$year) |>
  summarise(
    annual_probability = mean(.data$below_threshold),
    .groups = "drop"
  ) |>
  group_by(.data$model, .data$window) |>
  summarise(
    Years = paste(range(.data$year), collapse = "-"),
    `Mean annual probability below threshold` =
      mean(.data$annual_probability),
    .groups = "drop"
  )
mle_grid_projection_risk_ever <-
  mle_grid_projection_risk_data |>
  group_by(.data$model, .data$window, .data$iteration) |>
  summarise(
    crossed = any(.data$below_threshold),
    .groups = "drop"
  ) |>
  group_by(.data$model, .data$window) |>
  summarise(
    `Probability a draw crosses at least once` = mean(.data$crossed),
    .groups = "drop"
  )
left_join(
  mle_grid_projection_risk_annual,
  mle_grid_projection_risk_ever,
  by = c("model", "window")
) |>
  transmute(
    Model = .data$model,
    Window = .data$window,
    Years,
    `Relative TRO threshold` =
      percent(risk_relative_biomass_threshold, accuracy = 0.1),
    `Mean annual probability` = percent(
      .data$`Mean annual probability below threshold`,
      accuracy = 0.1
    ),
    `Ever-crossing probability` = percent(
      .data$`Probability a draw crosses at least once`,
      accuracy = 0.1
    )
  ) |>
  kable()

Table 28: Relative-TRO risk comparison for the balanced MCMC-grid base and sampled direct-M MLE grid.

Model	Window	Years	Relative TRO threshold	Mean annual probability	Ever-crossing probability
Balanced MCMC grid — nominal fleet allocations	Short term	2030-2032	20.0%	0.1%	0.1%
Balanced MCMC grid — nominal fleet allocations	Long term	2033-2036	20.0%	0.9%	1.8%
Sampled direct-M MLE grid — nominal fleet allocations	Short term	2030-2032	20.0%	0.0%	0.0%
Sampled direct-M MLE grid — nominal fleet allocations	Long term	2033-2036	20.0%	0.0%	0.0%

mle_grid_projection_tac_summary <- imap_dfr(
  mle_grid_projection_comparison,
  function(arm, model) {
    as_tibble(
      as.data.frame.table(
        arm$run$ctp_total_tac_iy,
        responseName = "tac"
      )
    ) |>
      transmute(
        model = model,
        year = as.integer(as.character(Year)),
        tac = .data$tac
      ) |>
      filter(.data$year %in% arm$run$projection_future_years) |>
      group_by(.data$model, .data$year) |>
      summarise(projection_quantiles(.data$tac), .groups = "drop")
  }
)
mle_grid_projection_tac_agreement <-
  mle_grid_projection_tac_summary |>
  filter(.data$model == projection_arm_labels[["base"]]) |>
  select(
    .data$year,
    mcmc_lower = .data$lower,
    mcmc_median = .data$median,
    mcmc_upper = .data$upper
  ) |>
  inner_join(
    mle_grid_projection_tac_summary |>
      filter(.data$model == projection_arm_labels[["mle_grid"]]) |>
      select(
        .data$year,
        mle_lower = .data$lower,
        mle_median = .data$median,
        mle_upper = .data$upper
      ),
    by = "year"
  ) |>
  mutate(
    median_difference = .data$mle_median - .data$mcmc_median,
    absolute_relative_median_difference =
      abs(.data$median_difference) / .data$mcmc_median
  )

if (nrow(mle_grid_projection_tac_agreement) !=
    length(projection_future_years)) {
  stop(
    "The MCMC-grid and direct-M MLE-grid TAC summaries do not cover the same years.",
    call. = FALSE
  )
}

mle_grid_tac_tolerance <- 1e-8
mle_grid_tac_first_different_year <-
  mle_grid_projection_tac_agreement |>
  filter(abs(.data$median_difference) > mle_grid_tac_tolerance) |>
  summarise(year = min(.data$year)) |>
  pull("year")
mle_grid_tac_matching_median_last_year <-
  if (is.finite(mle_grid_tac_first_different_year)) {
    mle_grid_tac_first_different_year - 1L
  } else {
    max(mle_grid_projection_tac_agreement$year)
  }
mle_grid_tac_max_median_difference <- max(
  abs(mle_grid_projection_tac_agreement$median_difference)
)
mle_grid_tac_max_relative_median_difference <- max(
  mle_grid_projection_tac_agreement$absolute_relative_median_difference
)
mle_grid_tac_interval_endpoints_identical <- all(
  abs(
    mle_grid_projection_tac_agreement$mle_lower -
      mle_grid_projection_tac_agreement$mcmc_lower
  ) <= mle_grid_tac_tolerance &
    abs(
      mle_grid_projection_tac_agreement$mle_upper -
        mle_grid_projection_tac_agreement$mcmc_upper
    ) <= mle_grid_tac_tolerance
)

ggplot(
  mle_grid_projection_tac_summary,
  aes(
    x = .data$year,
    y = .data$median,
    color = .data$model,
    fill = .data$model
  )
) +
  geom_ribbon(
    aes(ymin = .data$lower, ymax = .data$upper),
    alpha = 0.14,
    color = NA
  ) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 1.4, alpha = 0.85) +
  scale_color_manual(values = mle_grid_projection_colors) +
  scale_fill_manual(values = mle_grid_projection_colors) +
  scale_x_continuous(
    breaks = projection_future_years,
    minor_breaks = NULL
  ) +
  scale_y_zero(labels = comma) +
  labs(x = "Year", y = "Nominal total TAC", color = NULL, fill = NULL) +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Figure 20: Nominal total TAC for the balanced MCMC-grid base and sampled direct-M MLE-grid projections under the same nominal fleet allocation. Lines show medians, ribbons show 95% intervals, and points mark every year.

The two TAC projections are effectively indistinguishable under the same nominal fleet allocation. Total TAC is exactly the same during the fixed-TAC period (2026–2029). Across the CTP feedback years, the median remains identical through 2032; the largest subsequent difference is only 23.8 tonnes (0.089%). The 95% interval endpoints are identical in every projected year.

References

Commission for the Conservation of Southern Bluefin Tuna. 2026. Report of the Sixteenth Operating Model and Management Procedure Technical Meeting. OMMP16 meeting report. Commission for the Conservation of Southern Bluefin Tuna.

Hillary, R. M., A. L. Preece, N. Takahashi, C. R. Davies, and T. Itoh. 2023. The Southern Bluefin Tuna Stock Assessment in 2023. CCSBT-ESC/2308/16. Commission for the Conservation of Southern Bluefin Tuna. https://www.ccsbt.org/system/files/2023-08/ESC28_16_stockAssessment2023.pdf.