ESC31 Grid

Show code

options(readr.show_col_types = FALSE)
invisible(knitr::knit_meta_add(list(rmarkdown::html_dependency_jquery())))

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

theme_set(theme_bw())

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

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))
  )
}

Grid acceptance status

This page implements the approved nine-cell grid for steepness (\(h\)) and \(\psi\), the power parameter for relative reproductive contribution by age. In the model, \(\psi\) scales mature fish length as \(L^{3\psi}\) when calculating relative reproductive output at age; it is not recruitment variability (Hillary et al. 2023). The accepted base posterior is reused as mid cell 5, and the other eight cells are fitted and sampled independently. All nine selected cells now pass every MLE, sampler, provenance, and biological-state gate. The balanced 2,000-draw posterior, direct-get_M 108-fit MLE grid, and its 2,000-draw resample have also been built and independently validated. These accepted grid products are ready for projection use; compatible diagnostic and production results are reported on the separate Projections page.

Configuration

Show code

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)
}

source(file.path(esc_dir, "esc31_workflow.R"))
source(file.path(esc_dir, "esc31_inputs.R"), local = TRUE)
source(file.path(esc_dir, "esc31_base_posterior.R"), local = TRUE)
source(file.path(esc_dir, "esc31_grid_contract.R"), local = TRUE)

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

model_file <- file.path(run_dir, "esc31_base.rds")
base_mcmc_file <- model_file
grid_root <- file.path(run_dir, "grid_mcmc")
mle_grid_root <- file.path(run_dir, "grid_mle")
projection_run_dir <- file.path(run_dir, "projections")
dir.create(projection_run_dir, recursive = TRUE, showWarnings = FALSE)

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")
}

env_grid_cells <- function(name, n_cells) {
  value <- trimws(Sys.getenv(name, ""))
  if (!nzchar(value)) return(seq_len(n_cells))

  cells <- suppressWarnings(as.integer(trimws(strsplit(value, ",", fixed = TRUE)[[1]])))
  if (anyNA(cells) || any(cells < 1L | cells > n_cells) || anyDuplicated(cells)) {
    stop(
      name, " must be a comma-separated list of unique cell numbers from 1 to ",
      n_cells, ".",
      call. = FALSE
    )
  }
  cells
}

env_positive_integer <- function(name, default) {
  value <- suppressWarnings(as.integer(trimws(Sys.getenv(name, as.character(default)))))
  if (length(value) != 1L || is.na(value) || value < 1L) {
    stop(name, " must be one positive integer.", call. = FALSE)
  }
  value
}

if (!exists("refit_model", inherits = FALSE)) refit_model <- FALSE
run_grid_mcmc <- env_flag("ESC31_RUN_GRID_MCMC")
run_grid_postprocessing <- env_flag("ESC31_RUN_GRID_POSTPROCESSING")
run_m10_mle_grid <- env_flag("ESC31_RUN_MLE_GRID")
render_failed_grid_diagnostics <-
  env_flag("ESC31_RENDER_FAILED_GRID_DIAGNOSTICS")
if (run_grid_mcmc || run_grid_postprocessing || run_m10_mle_grid) {
  esc31_require_grid_specification_approval(
    esc31_grid_specification,
    stage = "grid computation"
  )
}
if (run_m10_mle_grid) {
  esc31_require_direct_mle_grid_specification_approval(
    esc31_direct_mle_grid_specification,
    stage = "direct-M MLE-grid computation"
  )
}
detected_physical_cores <- parallel::detectCores(logical = FALSE)
if (!is.finite(detected_physical_cores) || detected_physical_cores < 1L) {
  detected_physical_cores <- 1L
}
grid_state_diagnostic_cores <- env_positive_integer(
  "ESC31_STATE_DIAGNOSTIC_CORES",
  min(16L, detected_physical_cores)
)
grid_fit_validation_mode <- local({
  configured <- tolower(trimws(Sys.getenv(
    "ESC31_FIT_VALIDATION_MODE", "auto"
  )))
  if (!configured %in% c("auto", "full", "skip")) {
    stop(
      "`ESC31_FIT_VALIDATION_MODE` must be `auto`, `full`, or `skip`.",
      call. = FALSE
    )
  }
  configured
})
# Plot-summary construction is deterministic and each cell writes only its
# own portable fit, so this performance control is not a scientific input.
grid_posterior_summary_cores <- env_positive_integer(
  "ESC31_GRID_POSTERIOR_SUMMARY_CORES",
  min(3L, detected_physical_cores)
)
# Each production MLE cell is independent. Cell-level caches make this
# computational parallelism resumable without changing the grid design,
# optimizer, or acceptance contract.
m10_mle_grid_cores <- min(
  env_positive_integer(
    "ESC31_MLE_GRID_CORES",
    min(12L, detected_physical_cores)
  ),
  detected_physical_cores
)
rebuild_m10_mle_grid <- FALSE
check_m10_mle_estimability <- TRUE
m10_mle_gradient_limit <-
  esc31_direct_mle_grid_specification$max_gradient
grid_mle_gradient_limit <- esc31_grid_specification$optimizer$max_gradient
rebuild_combined_mcmc <- FALSE
grid_default_run_id <-
  "esc31_grid_dense_prod_20260722_v8_base_only_reuse_base_cell5_150_750"
selected_grid_run <- Sys.getenv(
  "ESC31_SELECTED_GRID_RUN",
  esc31_grid_mixed_source_selection$output_run_id
)
grid_run_id <- trimws(Sys.getenv(
  "ESC31_GRID_RUN_ID",
  grid_default_run_id
))
if (!grepl("^esc31_grid_[A-Za-z0-9_-]+$", grid_run_id)) {
  stop(
    "ESC31_GRID_RUN_ID must start with 'esc31_grid_' and contain only letters, numbers, '_' or '-'.",
    call. = FALSE
  )
}
grid_output_dir <- file.path(grid_root, grid_run_id)
combined_mcmc_draws <- esc31_grid_specification$combined_posterior$draws
combined_mcmc_seed <- esc31_grid_specification$combined_posterior$seed
combined_mcmc_file <- file.path(
  projection_run_dir,
  paste0("esc31_projection_mcmc_", combined_mcmc_draws, ".rda")
)
m10_mle_grid_file <- file.path(
  mle_grid_root,
  "m0_m10_posterior_quantile_mle_grid_108.rds"
)

grid_values <- expand.grid(
  h = esc31_grid_specification$coordinates$h,
  psi = esc31_grid_specification$coordinates$psi
)
grid_base_cell <- esc31_grid_specification$coordinates$reused_base_cell
grid_sampled_cells <- esc31_grid_specification$coordinates$sampled_cells
if (!identical(grid_base_cell, 5L) ||
    !identical(grid_sampled_cells, setdiff(seq_len(nrow(grid_values)), 5L)) ||
    !isTRUE(all.equal(
      as.numeric(grid_values[grid_base_cell, c("h", "psi")]),
      c(
        esc31_grid_specification$base_cell_reuse$h,
        esc31_grid_specification$base_cell_reuse$psi
      ),
      tolerance = 0,
      check.attributes = FALSE
    ))) {
  stop("The approved base-cell reuse does not match grid cell 5.",
       call. = FALSE)
}
grid_mcmc_cells <- env_grid_cells("ESC31_GRID_CELLS", nrow(grid_values))
grid_mcmc_num_samples <- env_positive_integer(
  "ESC31_GRID_NUM_SAMPLES",
  esc31_grid_specification$sampler$num_samples
)
grid_failed_cell_rerun_active <-
  grid_mcmc_num_samples != esc31_grid_specification$sampler$num_samples
if (grid_failed_cell_rerun_active) {
  allowed_rerun_cells <- c(
    esc31_grid_failed_cell_rerun$failed_source_cells,
    grid_base_cell
  )
  if (!run_grid_mcmc ||
      !identical(
        as.integer(grid_mcmc_num_samples),
        as.integer(
          esc31_grid_failed_cell_rerun$sampler$num_samples
        )
      ) ||
      !identical(
        grid_run_id,
        esc31_grid_failed_cell_rerun$output_run_id
      ) ||
      any(!grid_mcmc_cells %in% allowed_rerun_cells)) {
    stop(
      "A non-default grid chain length is approved only for the recorded ",
      "1,500-draw reruns of failed cells 1, 2, 4, 6, 7, and 9 ",
      "(with cell 5 permitted only to materialize the reused base fit).",
      call. = FALSE
    )
  }
}
m10_grid_range_probs <-
  esc31_direct_mle_grid_specification$m10_posterior_probabilities
m10_grid_n <- esc31_direct_mle_grid_specification$m10_values
m0_grid_probs <-
  esc31_direct_mle_grid_specification$m0_posterior_probabilities
m0_grid_n <- esc31_direct_mle_grid_specification$m0_values
mle_grid_h_values <-
  esc31_direct_mle_grid_specification$h_values
mle_grid_psi_values <-
  esc31_direct_mle_grid_specification$psi_values
mle_grid_resample_n <-
  esc31_direct_mle_grid_specification$resample_draws
mle_grid_resample_seed <-
  esc31_direct_mle_grid_specification$resample_seed
m10_mle_tmbfit_file <- file.path(
  mle_grid_root,
  paste0(
    "m0_m10_posterior_quantile_mle_grid_",
    mle_grid_resample_n,
    "_tmbfit.rds"
  )
)

pair_plot_cells <- seq_len(nrow(grid_values))
rhat_acceptance_limit <- esc31_grid_specification$acceptance$max_rhat
ess_acceptance_limit <- esc31_grid_specification$acceptance$min_bulk_ess

Approved Grid Design

Show code

approved_grid_cells <- grid_values |>
  mutate(
    Cell = row_number(),
    Role = if_else(
      .data$Cell == grid_base_cell,
      "Accepted base posterior",
      "New MLE and MCMC"
    )
  ) |>
  select(.data$Cell, .data$h, .data$psi, .data$Role)

approved_grid_cells |>
  mutate(
    h = formatC(.data$h, format = "f", digits = 2L),
    psi = formatC(.data$psi, format = "f", digits = 2L)
  ) |>
  kable(align = c("r", "r", "r", "l")) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  ) |>
  kableExtra::row_spec(grid_base_cell, background = "#DDEBF7", bold = TRUE)

Table 1: Approved nine-cell MCMC grid. Cell 5 is the accepted base posterior; the other eight cells require a new MLE and MCMC.

Cell	h	psi	Role
1	0.60	1.50	New MLE and MCMC
2	0.70	1.50	New MLE and MCMC
3	0.80	1.50	New MLE and MCMC
4	0.60	1.75	New MLE and MCMC
5	0.70	1.75	Accepted base posterior
6	0.80	1.75	New MLE and MCMC
7	0.60	2.00	New MLE and MCMC
8	0.70	2.00	New MLE and MCMC
9	0.80	2.00	New MLE and MCMC

Show code

accepted_grid_samples_per_chain <-
  esc31_grid_mixed_source_selection$retained_draws_per_chain
accepted_grid_cell_ids <-
  suppressWarnings(as.integer(names(accepted_grid_samples_per_chain)))
if (length(accepted_grid_samples_per_chain) != nrow(grid_values) ||
    anyNA(accepted_grid_cell_ids) ||
    !setequal(accepted_grid_cell_ids, seq_len(nrow(grid_values)))) {
  stop(
    "The approved mixed-source sampler contract does not cover all grid cells.",
    call. = FALSE
  )
}

format_grid_draw_groups <- function(multiplier = 1L) {
  paste(
    vapply(
      sort(unique(as.integer(accepted_grid_samples_per_chain))),
      function(draws_per_chain) {
        cells <- accepted_grid_cell_ids[
          accepted_grid_samples_per_chain == draws_per_chain
        ]
        paste0(
          formatC(
            draws_per_chain * multiplier,
            format = "f", digits = 0L, big.mark = ","
          ),
          " (cells ", paste(cells, collapse = ", "), ")"
        )
      },
      character(1L)
    ),
    collapse = "; "
  )
}

tribble(
  ~Setting, ~Value,
  "Chains per cell", formatC(
    esc31_grid_specification$sampler$chains,
    format = "f", digits = 0L, big.mark = ","
  ),
  "Warmup per chain", formatC(
    esc31_grid_specification$sampler$num_warmup,
    format = "f", digits = 0L, big.mark = ","
  ),
  "Retained draws per chain", format_grid_draw_groups(),
  "Retained draws per cell", format_grid_draw_groups(
    esc31_grid_mixed_source_selection$chains
  ),
  "Total retained draws across grid", formatC(
    sum(
      accepted_grid_samples_per_chain *
        esc31_grid_mixed_source_selection$chains
    ),
    format = "f", digits = 0L, big.mark = ","
  ),
  "Metric", stringr::str_to_title(esc31_grid_specification$sampler$metric),
  "Adapt delta", formatC(
    esc31_grid_specification$sampler$adapt_delta,
    format = "f", digits = 3L
  ),
  "Maximum treedepth", formatC(
    esc31_grid_specification$sampler$max_treedepth,
    format = "f", digits = 0L
  ),
  "Chain initialization", "Exact cell MLE",
  "Balanced combined posterior", paste0(
    formatC(
      esc31_grid_specification$combined_posterior$draws,
      format = "f", digits = 0L, big.mark = ","
    ),
    " draws"
  ),
  "Dependent MLE grid", paste0(
    formatC(
      esc31_direct_mle_grid_specification$cells,
      format = "f", digits = 0L, big.mark = ","
    ),
    " cells"
  )
) |>
  kable(align = c("l", "l")) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE
  )

Table 2: Approved mixed-source grid sampling and downstream design.

Setting	Value
Chains per cell	4
Warmup per chain	150
Retained draws per chain	750 (cells 3, 5, 8); 1,500 (cells 1, 2, 4, 6, 7, 9)
Retained draws per cell	3,000 (cells 3, 5, 8); 6,000 (cells 1, 2, 4, 6, 7, 9)
Total retained draws across grid	45,000
Metric	Dense
Adapt delta	0.999
Maximum treedepth	13
Chain initialization	Exact cell MLE
Balanced combined posterior	2,000 draws
Dependent MLE grid	108 cells

Load Base Fit

Show code

if (run_grid_mcmc || run_grid_postprocessing || run_m10_mle_grid) {
  esc31_require_ll34_conditioning_approval(
    base_model_contract,
    stage = "grid 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 approved grid 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))
grid_harvest_wall_identity <- esc31_grid_harvest_wall_identity(
  base_model_contract
)
base_wall_data <- lapply(
  names(grid_harvest_wall_identity$data),
  function(name) base_fit$data[[name]]
)
names(base_wall_data) <- names(grid_harvest_wall_identity$data)
if (!identical(base_wall_data, grid_harvest_wall_identity$data) ||
    !identical(
      base_fit$provenance$metadata$run_identity$config$
        base_model_contract$harvest_wall,
      grid_harvest_wall_identity$decision
    )) {
  stop(
    "The current base fit does not carry the approved harvest-wall contract.",
    call. = FALSE
  )
}
base_grid_posterior_fit <- sbt_fit_read(
  base_mcmc_file,
  strict = TRUE,
  rebuild = FALSE
)
base_grid_posterior_contract <- esc31_base_mcmc_contract(
  esc_dir,
  base_fit
)
if (!esc31_mcmc_identity_equivalent(
      base_grid_posterior_fit$provenance$metadata$mcmc_run_identity,
      base_grid_posterior_contract$identity,
      fit = base_grid_posterior_fit,
      base_fit = base_fit
    )) {
  stop(
    "The current base posterior does not match its approved MCMC contract.",
    call. = FALSE
  )
}
grid_prerequisite_acceptance <- esc31_build_grid_prerequisite_acceptance(
  esc_dir = esc_dir,
  fits = list(base = base_grid_posterior_fit),
  files = c(base = base_mcmc_file),
  base_run_signature =
    base_fit$provenance$metadata$run_identity$signature,
  harvest_wall_identity = grid_harvest_wall_identity
)
esc31_atomic_save_grid_prerequisite_acceptance(
  grid_prerequisite_acceptance,
  esc31_grid_prerequisite_acceptance_file(esc_dir)
)
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 = grid_harvest_wall_identity
)
grid_prerequisite_validation <- grid_prerequisite_check$provenance

base_grid_posterior_payload_checksum <-
  sbt:::.grid_mcmc_posterior_payload_checksum(
    as_tmbfit(base_grid_posterior_fit)
  )
base_grid_reuse_identity <- list(
  schema_version = 1L,
  decision = esc31_grid_specification$base_cell_reuse,
  cell = grid_base_cell,
  grid_values = as.list(grid_values[grid_base_cell, , drop = FALSE]),
  source_file = basename(base_mcmc_file),
  source_file_md5 = unname(tools::md5sum(base_mcmc_file)),
  source_runtime_id = base_grid_posterior_fit$runtime_id,
  source_model_signature = base_grid_posterior_fit$model$signature,
  source_fit_scientific_signature =
    esc31_fit_scientific_signature(base_grid_posterior_fit),
  source_mle_run_signature =
    base_grid_posterior_fit$provenance$metadata$run_identity$signature,
  source_mcmc_run_signature =
    base_grid_posterior_fit$provenance$metadata$mcmc_run_identity$signature,
  source_posterior_payload_checksum =
    base_grid_posterior_payload_checksum,
  source_sampler = base_grid_posterior_fit$mcmc$settings$sampler,
  source_acceptance_thresholds =
    base_grid_posterior_fit$mcmc$settings$acceptance_thresholds,
  implementation = "esc31_grid_reuse_accepted_base_cell_v1"
)
if (!identical(
      base_grid_reuse_identity$source_file_md5,
      unname(grid_prerequisite_validation$artifact_checksums[["base"]])
    ) ||
    !esc31_mcmc_identity_equivalent(
      base_grid_posterior_fit$provenance$metadata$mcmc_run_identity,
      base_grid_posterior_contract$identity,
      fit = base_grid_posterior_fit,
      base_fit = base_fit
    ) ||
    !identical(
      base_grid_reuse_identity$source_sampler,
      base_grid_posterior_contract$sampler
    ) ||
    !identical(
      base_grid_reuse_identity$source_acceptance_thresholds,
      base_grid_posterior_contract$thresholds
    ) ||
    !identical(
      as.integer(base_grid_reuse_identity$source_sampler$seed),
      esc31_grid_specification$base_cell_reuse$source_sampler_seed
    ) ||
    !identical(base_grid_posterior_fit$data, base_fit$data) ||
    !identical(base_grid_posterior_fit$parameters, base_fit$parameters) ||
    !identical(base_grid_posterior_fit$map, base_fit$map)) {
  stop(
    "The accepted base posterior does not match the approved grid-cell reuse contract.",
    call. = FALSE
  )
}

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
)
target_par_names <- expand_parameter_names(names(obj_projection$par))
target_sample_names <- c(target_par_names, "lp__")

Helper Functions

Show code

grid_cell_pattern <- "^grid([0-9]+)(\\.sbt\\.rds|\\.rda)$"

cell_number <- function(path) {
  as.integer(sub(grid_cell_pattern, "\\1", basename(path), perl = TRUE))
}

validate_grid_dir <- function(path) {
  if (!dir.exists(path)) {
    stop("Grid run directory does not exist: ", path, call. = FALSE)
  }

  grid_files <- list.files(path, pattern = grid_cell_pattern, full.names = TRUE)
  grid_cells <- cell_number(grid_files)
  diagnostics_file <- file.path(path, "grid_mcmc_diagnostics.csv")
  complete_cells <- length(grid_files) == 9L &&
    !anyNA(grid_cells) &&
    !anyDuplicated(grid_cells) &&
    identical(sort(grid_cells), seq_len(9L))
  if (!complete_cells || !file.exists(diagnostics_file)) {
    stop("Grid run is not a completed nine-cell MCMC directory: ", path, call. = FALSE)
  }
  normalizePath(path, mustWork = TRUE)
}

grid_signature <- function(files) {
  esc31_object_md5(esc31_md5_manifest(files, basename(files)))
}

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))
}

atomic_save_rdata <- function(objects, path) {
  if (!is.list(objects) || is.null(names(objects)) || any(!nzchar(names(objects)))) {
    stop("objects must be a named list.", call. = FALSE)
  }
  cache_env <- list2env(objects, parent = emptyenv())
  atomic_cache_write(
    path,
    function(temporary) save(
      list = names(objects),
      envir = cache_env,
      file = temporary
    )
  )
}

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

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
}

format_diagnostic_cell <- function(x, fails, digits = 3) {
  out <- format_decimal(x, digits = digits)
  flagged <- !is.na(fails) & fails
  out[flagged] <- paste0(
    "<span style=\"display:block; margin:-8px; padding:8px; ",
    "background-color:#f8d7da; color:#842029; font-weight:600;\">",
    out[flagged],
    "</span>"
  )
  out
}

format_flagged_text <- function(x, fails) {
  out <- as.character(x)
  flagged <- is.na(fails) | fails
  out[flagged] <- paste0(
    "<span style=\"display:block; margin:-8px; padding:8px; ",
    "background-color:#f8d7da; color:#842029; font-weight:600;\">",
    out[flagged],
    "</span>"
  )
  out
}

replace_missing <- function(x, missing = "-") {
  x <- as.data.frame(x)
  x[] <- lapply(x, function(column) {
    column <- as.character(column)
    column[is.na(column) | column == "NA"] <- missing
    column
  })
  as_tibble(x)
}

fit_expected_color <- "#0072B2"
lev_vars <- c("h", "M0", "M10", "psi")

resolve_grid_file <- function(file, selected_grid_dir = NULL) {
  if (is.null(selected_grid_dir) || !grepl(grid_cell_pattern, basename(file))) {
    stop("Grid files must resolve within the selected grid run.", call. = FALSE)
  }
  selected_grid_dir <- normalizePath(selected_grid_dir, mustWork = TRUE)
  candidate <- file.path(selected_grid_dir, basename(file))
  if (!file.exists(candidate)) {
    stop("Could not find grid output file: ", file, call. = FALSE)
  }
  candidate <- normalizePath(candidate, mustWork = TRUE)
  if (!identical(dirname(candidate), selected_grid_dir)) {
    stop("Grid output escaped the selected grid run: ", file, call. = FALSE)
  }
  candidate
}

load_grid_output <- function(file, fallback_cell = NA_integer_) {
  sbt:::load_grid_mcmc_cell(
    resolve_grid_file(file, selected_grid_dir = grid_dir),
    grid_values = grid_values,
    fallback_cell = fallback_cell
  )
}

read_grid_run_metadata <- function(path) {
  if (!dir.exists(path)) return(NULL)

  metadata_file <- file.path(path, "run_metadata.rds")
  if (file.exists(metadata_file)) {
    metadata <- read_rds_cache(metadata_file)
    if (is.list(metadata)) return(metadata)
  }

  summary_file <- file.path(path, "grid_mcmc_summary.rds")
  if (file.exists(summary_file)) {
    summary <- read_rds_cache(summary_file)
    if (is.list(summary) && is.list(summary$run_metadata)) {
      return(summary$run_metadata)
    }
  }

  cell_files <- list.files(
    path,
    pattern = "^grid[0-9]+\\.sbt\\.rds$",
    full.names = TRUE
  )
  if (length(cell_files)) {
    cell <- sbt:::load_grid_mcmc_cell(cell_files[[1]], grid_values = grid_values)
    if (is.list(cell$run_metadata)) return(cell$run_metadata)
  }

  NULL
}

grid_input_signatures <- function(base_fit, grid_values) {
  lower <- base_fit$bounds$lower
  upper <- base_fit$bounds$upper
  shared_names <- base_fit$bounds$parameter %||% NULL
  lower_names <- sbt:::.grid_mcmc_bound_names(
    lower, shared_names, "lower"
  )
  upper_names <- sbt:::.grid_mcmc_bound_names(
    upper, shared_names, "upper"
  )
  if (!identical(lower_names, upper_names)) {
    stop("The lower and upper grid-bound parameter names differ.",
         call. = FALSE)
  }
  bound_signature_input <- list(
    parameter = lower_names %||% character(),
    lower = lower,
    upper = upper
  )
  par_map <- c(
    h = "par_log_h", psi = "par_log_psi",
    m0 = "par_log_m0", m10 = "par_log_m10"
  )
  cell_inputs <- lapply(seq_len(nrow(grid_values)), function(cell) {
    sbt:::.grid_mcmc_cell_inputs(
      parameters = base_fit$parameters,
      map = base_fit$map,
      grid_row = grid_values[cell, , drop = FALSE],
      par_map = par_map
    )
  })
  list(
    grid = sbt:::.grid_mcmc_signature(list(
      data = base_fit$data,
      parameters = base_fit$parameters,
      map = base_fit$map,
      random = as.character(base_fit$random),
      grid = grid_values,
      bounds = bound_signature_input
    )),
    cells = vapply(seq_len(nrow(grid_values)), function(cell) {
      sbt:::.grid_mcmc_signature(list(
        data = base_fit$data,
        parameters = cell_inputs[[cell]]$parameters,
        map = cell_inputs[[cell]]$map,
        random = as.character(base_fit$random),
        bounds = bound_signature_input
      ))
    }, character(1L)),
    inputs = cell_inputs
  )
}

materialize_reused_base_grid_cell <- function(
    source_fit, source_identity, base_fit, grid_values, cell,
    output_file, run_settings, run_signature, run_metadata,
    requested_cells) {
  signatures <- grid_input_signatures(base_fit, grid_values)
  expected_grid_values <- as.list(grid_values[cell, , drop = FALSE])
  expected_cell_inputs <- signatures$inputs[[cell]]
  validate_existing <- function() {
    existing <- sbt:::.grid_mcmc_existing_fit(
      output_file,
      cell = cell,
      grid_row = grid_values[cell, , drop = FALSE],
      input_signature = signatures$cells[[cell]],
      grid_input_signature = signatures$grid,
      run_signature = run_signature,
      require_run_signature = TRUE
    )
    saved <- existing$fit$provenance$metadata
    if (!identical(saved$grid_reuse, source_identity) ||
        !identical(saved$grid$reuse, source_identity) ||
        !identical(
          sbt:::.grid_mcmc_posterior_payload_checksum(existing$fit$mcmc),
          source_identity$source_posterior_payload_checksum
        )) {
      stop(
        "The existing reused base grid cell has different source provenance.",
        call. = FALSE
      )
    }
    existing
  }

  if (file.exists(output_file)) {
    existing <- validate_existing()
    message("Skipping validated reused base posterior in grid cell ", cell, ".")
    return(invisible(existing))
  }
  if (!identical(source_fit$data, base_fit$data) ||
      !identical(source_fit$parameters, expected_cell_inputs$parameters) ||
      !identical(source_fit$map, expected_cell_inputs$map) ||
      !identical(source_fit$random, base_fit$random)) {
    stop(
      "The accepted base posterior is not the exact model represented by grid cell ",
      cell, ".",
      call. = FALSE
    )
  }

  source_mcmc <- as_tmbfit(source_fit)
  source_object <- sbt_obj(source_fit, fresh = TRUE)
  source_object$par <- source_fit$fit$opt$par
  source_object$env$last.par.best <- source_fit$fit$opt$par
  invisible(source_object$fn(source_fit$fit$opt$par))
  source_mle_state <- source_fit$fit$diagnostics$biological_state_mle
  source_posterior_state <-
    source_fit$fit$diagnostics$biological_state_posterior
  source_nlminb_passes <-
    as.integer(source_fit$fit$diagnostics$optimizer$n_passes %||% NA_integer_)
  diagnostic <- sbt:::diagnose_grid_mcmc(
    mcmc = source_mcmc,
    cell = cell,
    grid_cell = grid_values[cell, , drop = FALSE],
    opt = source_fit$fit$opt,
    start_fn = source_fit$fit$opt$objective,
    b0_start_multiplier = 1,
    nlminb_passes = source_nlminb_passes,
    estimability = source_fit$fit$estimability,
    file = output_file,
    max_treedepth = source_fit$mcmc$max_treedepth,
    run_signature = run_signature,
    state_diagnostics = source_posterior_state
  )
  state_validation <- sbt:::.grid_mcmc_validate_biological_state(
    mle_state = source_mle_state,
    posterior_state = source_posterior_state,
    diagnostic = diagnostic,
    mcmc = source_fit$mcmc
  )
  if (!isTRUE(state_validation$passes)) {
    stop(
      "The accepted base posterior cannot be reused as grid cell ", cell,
      ": ", state_validation$message, ".",
      call. = FALSE
    )
  }

  reused_fit <- sbt_fit(
    data = source_fit$data,
    obj = source_object,
    bounds = source_fit$bounds,
    control = source_fit$control,
    opt = source_fit$fit$opt,
    estimability = source_fit$fit$estimability,
    diagnostics = utils::modifyList(
      source_fit$fit$diagnostics,
      list(grid_mcmc = diagnostic)
    ),
    metadata = utils::modifyList(
      source_fit$provenance$metadata,
      list(
        grid = list(
          cell = cell,
          values = expected_grid_values,
          input_signature = signatures$cells[[cell]],
          grid_input_signature = signatures$grid,
          reuse = source_identity
        ),
        start = list(
          objective = source_fit$fit$opt$objective,
          b0_multiplier = 1
        ),
        run = list(
          signature = run_signature,
          settings = run_settings,
          requested_cells = requested_cells
        ),
        run_metadata = run_metadata,
        grid_reuse = source_identity
      )
    ),
    mcmc = source_fit$mcmc,
    mcmc_settings = source_fit$mcmc$settings,
    makeadfun_args = source_fit$makeadfun_args,
    model = source_fit$model$name
  )
  reused_fit$validation <- source_fit$validation
  sbt_fit_validate(reused_fit)
  sbt:::.save_grid_sbt_fit(reused_fit, output_file)
  validate_existing()
  message(
    "Materialized accepted base posterior as canonical grid cell ", cell, "."
  )
  invisible(TRUE)
}

MCMC grid

Run grid

Decision esc31_2026_h_psi_dense_grid_v8_base_only_reuse_base_cell5_sparse_nuts_150_750 approves the configuration below after the current-identity base MCMC passes its numeric and biological-state gates. Sensitivity posteriors are comparisons and do not gate either grid. The base-only acceptance is rebuilt and machine-validated at launch. The accepted base posterior is exactly grid cell 5 (h = 0.7, psi = 1.75), so it is retained as that grid element with its original seed, sampler identity, and posterior payload checksum. Only the other eight cells are newly optimized and sampled. The stable output directory is selected with ESC31_GRID_RUN_ID; rerunning with the same ID validates and skips compatible completed cells. ESC31_GRID_CELLS can restrict a launch to comma-separated cell numbers while filling an incomplete run and preserving the original cell numbers (for example, ESC31_GRID_CELLS=1,4,8). Existing cell files are never overwritten. If any newly sampled cell fails a production diagnostic, a new run identity is required before changing the sampling settings. The first completed grid failed only the R-hat gate in cells 1, 2, 4, 6, 7, and 9. Decision esc31_2026_h_psi_grid_failed_cells_sparse_nuts_150_1500_v1 therefore authorizes 1,500 retained draws per chain for those six cells while preserving their seeds and all other controls. Cells 3, 5, and 8 are not resampled. The six longer reruns pass every numeric, sampler-pathology, and saved biological-state gate. Selection esc31_2026_grid_mixed_source_selection_failed_cells_1500_v1 therefore combines those six cells with original passing cells 3, 5, and 8. The selection copies every portable source sbt_fit unchanged and binds its source-run signature, posterior checksum, state checksum, retained-draw count, and source file checksum in a nine-cell manifest. Each source fit stores the complete MLE and posterior validation performed at fit completion. This page verifies that payload-bound record and the cross-file grid manifest; it does not normally rebuild nine objectives and rescan every retained state. A deliberate full audit remains available with ESC31_FIT_VALIDATION_MODE=full. The accepted base remains the central element only while its exact source posterior continues to pass the contract. New cells are saved as portable grid<N>.sbt.rds sbt_fit objects. Historical grid<N>.rda cells remain readable for audit and comparison only. Each cell is checked for a finite starting objective, successful nlminb() convergence, and estimability before MCMC sampling starts. The production starting settings use four chains with 150 SparseNUTS warmup and 750 retained draws per chain. The earlier 500-retained grid had a weakest tail ESS of 161, and the selected base’s 500-retained run narrowly failed maximum R-hat; consequently, a cell is not accepted unless the revised run clears every gate. The recorded seed defines independent chain RNG streams, while every chain starts at the cell MLE’s exact last.par.best. This replaces dispersed Gaussian starts that were finite but could be catastrophically biologically infeasible under the combined-harvest wall. Identical mode starts trade some initial-overdispersion diagnostic power for stable initialization; the strict rank-plot, R-hat, ESS, sampler-pathology, and biological-state gates therefore remain mandatory. These settings do not guarantee convergence; the diagnostic thresholds below still determine whether each cell can be combined. The sampler resume signature receives only scientific metadata; the complete workflow and sibling-source manifests remain in the run-directory provenance sidecar.

Grid postprocessing and the 108-cell MLE computation retain the base-only acceptance manifest and recorded grid provenance. Every dependent stage still stops if any of the nine grid cells fails its production gates.

All nine grid elements and all 108 dependent MLE cells use the approved preventive harvest wall in their fitted objective. The wall is recorded separately from the feasibility continuation: a positive wall contribution is expected and is not a state-gate failure, whereas every accepted cell must retain zero continuation penalty within tolerance.

Show code

grid_mcmc_config <- list(
  grid_values = grid_values,
  control = base_fit$control,
  n_passes = esc31_grid_specification$optimizer$n_passes,
  b0_start_step = esc31_grid_specification$optimizer$b0_start_step,
  b0_start_max = esc31_grid_specification$optimizer$b0_start_max,
  check = TRUE,
  stop_on_failure = TRUE,
  metric = esc31_grid_specification$sampler$metric,
  num_samples = grid_mcmc_num_samples,
  num_warmup = esc31_grid_specification$sampler$num_warmup,
  chains = esc31_grid_specification$sampler$chains,
  cores = esc31_grid_specification$sampler$cores,
  adapt_delta = esc31_grid_specification$sampler$adapt_delta,
  max_treedepth = esc31_grid_specification$sampler$max_treedepth,
  seed = esc31_grid_specification$sampler$seed,
  init = esc31_grid_specification$sampler$init,
  refresh = esc31_grid_specification$sampler$refresh
)
grid_mcmc_sampler_args <- list(
  metric = grid_mcmc_config$metric,
  num_samples = grid_mcmc_config$num_samples,
  num_warmup = grid_mcmc_config$num_warmup,
  chains = grid_mcmc_config$chains,
  cores = grid_mcmc_config$cores,
  seed = grid_mcmc_config$seed,
  control = list(
    adapt_delta = grid_mcmc_config$adapt_delta,
    max_treedepth = grid_mcmc_config$max_treedepth
  ),
  init = grid_mcmc_config$init,
  globals = sbt_globals(),
  refresh = grid_mcmc_config$refresh
)
grid_record_identity <- esc31_workflow_identity(
  esc_dir,
  workflow_files = "esc31_grid_contract.R",
  sbt_entry_points = esc31_sbt_entry_points("grid"),
  config = list(
    stage = "grid_mcmc",
    grid_specification = esc31_grid_specification,
    prerequisite_acceptance = grid_prerequisite_validation,
    base_fit_md5 = base_fit_md5,
    base_run_signature =
      base_fit$provenance$metadata$run_identity$signature,
    grid = grid_mcmc_config,
    failed_cell_rerun_decision = if (grid_failed_cell_rerun_active) {
      esc31_grid_failed_cell_rerun
    } else {
      NULL
    },
    fit_workflow_version =
      "esc31_grid_mcmc_v8_reuse_accepted_base_cell5",
    base_cell_reuse = base_grid_reuse_identity,
    biological_state_contract = biological_state_contract(),
    harvest_wall_contract = grid_harvest_wall_identity$executable,
    harvest_wall_decision = grid_harvest_wall_identity$decision
  )
)
grid_signature_inputs <- list(
  schema_version = 9L,
  stage = "grid_mcmc",
  grid_specification = esc31_grid_specification,
  prerequisite_acceptance = grid_prerequisite_validation,
  base_fit_md5 = base_fit_md5,
  base_run_signature = base_fit$provenance$metadata$run_identity$signature,
  raw_input_md5 = grid_record_identity$raw_input_md5,
  sbt_function_identity = grid_record_identity$sbt_function_identity,
  # Retain the normalized content fingerprint field used by accepted grid
  # metadata while the full workflow identity records the sbtdata release.
  sbt_data_md5 =
    esc31_workflow_scientific_data_identity(grid_record_identity),
  scientific_package_versions = grid_record_identity$scientific_package_versions,
  grid = grid_mcmc_config,
  failed_cell_rerun_decision = if (grid_failed_cell_rerun_active) {
    esc31_grid_failed_cell_rerun
  } else {
    NULL
  },
  base_cell_reuse = base_grid_reuse_identity,
  biological_state_contract = biological_state_contract(),
  harvest_wall_contract = grid_harvest_wall_identity$executable,
  harvest_wall_decision = grid_harvest_wall_identity$decision,
  harvest_wall_data = grid_harvest_wall_identity$data,
  implementation =
    "esc31_grid_driver_v10_reuse_accepted_base_cell5_harvest_wall_bound"
)
if (!grid_failed_cell_rerun_active) {
  # Preserve the exact original run signature when the optional rerun
  # provenance does not apply.
  grid_signature_inputs$failed_cell_rerun_decision <- NULL
}
grid_run_metadata <- list(
  schema_version = 9L,
  stage = "grid_mcmc",
  run_id = grid_run_id,
  signature_inputs = grid_signature_inputs,
  signature = esc31_object_md5(grid_signature_inputs),
  prerequisite_review_confirmation = list(
    confirmed = isTRUE(grid_prerequisite_check$base_posterior_accepted),
    scope = "accepted_base_mcmc_only",
    mechanism = "esc31_base_only_grid_launch_v1"
  ),
  failed_cell_rerun_decision = if (grid_failed_cell_rerun_active) {
    esc31_grid_failed_cell_rerun
  } else {
    NULL
  },
  recorded_workflow = grid_record_identity
)
grid_sampler_run_metadata <- grid_run_metadata[c(
  "schema_version", "stage", "run_id", "signature_inputs", "signature"
)]
expected_grid_run_settings <- sbt:::.grid_mcmc_run_settings(
  control = grid_mcmc_config$control,
  n_passes = grid_mcmc_config$n_passes,
  b0_start_step = grid_mcmc_config$b0_start_step,
  b0_start_max = grid_mcmc_config$b0_start_max,
  check = grid_mcmc_config$check,
  stop_on_failure = grid_mcmc_config$stop_on_failure,
  mcmc_args = grid_mcmc_sampler_args,
  run_metadata = grid_sampler_run_metadata
)
expected_grid_run_signature <- sbt:::.grid_mcmc_signature(
  expected_grid_run_settings
)

if (run_grid_mcmc) {
  existing_grid_metadata <- read_grid_run_metadata(grid_output_dir)
  if (!is.null(existing_grid_metadata) &&
      !identical(existing_grid_metadata$signature, grid_run_metadata$signature)) {
    stop(
      "ESC31_GRID_RUN_ID identifies a run with different scientific provenance. ",
      "Use the original code and controls to resume it, or choose a new run ID.",
      call. = FALSE
    )
  }
  dir.create(grid_output_dir, recursive = TRUE, showWarnings = FALSE)
  atomic_save_rds(
    grid_run_metadata,
    file.path(grid_output_dir, "run_metadata.rds")
  )
  materialize_reused_base_grid_cell(
    source_fit = base_grid_posterior_fit,
    source_identity = base_grid_reuse_identity,
    base_fit = base_fit,
    grid_values = grid_values,
    cell = grid_base_cell,
    output_file = file.path(
      grid_output_dir,
      paste0("grid", grid_base_cell, ".sbt.rds")
    ),
    run_settings = expected_grid_run_settings,
    run_signature = expected_grid_run_signature,
    run_metadata = grid_sampler_run_metadata,
    requested_cells = grid_mcmc_cells
  )

  grid_args <- list(
    data = base_fit$data,
    parameters = base_fit$parameters,
    grid = grid_values,
    bounds = base_fit$bounds,
    map = base_fit$map,
    grid_dir = grid_output_dir,
    file_prefix = "grid",
    control = grid_mcmc_config$control,
    n_passes = grid_mcmc_config$n_passes,
    b0_start_step = grid_mcmc_config$b0_start_step,
    b0_start_max = grid_mcmc_config$b0_start_max,
    check = grid_mcmc_config$check,
    stop_on_failure = grid_mcmc_config$stop_on_failure,
    overwrite = FALSE,
    cells = grid_mcmc_cells,
    resume = TRUE,
    run_metadata = grid_sampler_run_metadata,
    mcmc_args = grid_mcmc_sampler_args
  )
  grid_run <- do.call(sbt::run_grid_mcmc, grid_args)
  atomic_save_rds(
    grid_run_metadata,
    file.path(grid_output_dir, "run_metadata.rds")
  )
}

Grid acceptance results

Show code

selected_grid_dir <- if (run_grid_mcmc) {
  grid_output_dir
} else {
  file.path(grid_root, selected_grid_run)
}
grid_dir <- validate_grid_dir(selected_grid_dir)
grid_posterior_files <- file.path(
  grid_dir,
  paste0("grid", seq_len(nrow(grid_values)), ".sbt.rds")
)
prepare_grid_posterior_summary <- function(file) {
  fit <- sbt_fit_read(file, strict = TRUE, rebuild = FALSE)
  # This page does not draw the pooled conventional-tag panel. Reuse and
  # validate an existing complete plot summary instead of expanding every
  # portable grid fit solely to add that later optional derived component.
  if (!is.null(fit$posterior_summaries$entries$plots)) {
    get_posterior_summary(
      fit,
      name = "plots",
      probs = c(0.025, 0.975),
      transform = "plots"
    )
    return(invisible(TRUE))
  }
  summaries_before <- fit$posterior_summaries
  fit <- add_posterior_summaries(fit)
  if (!identical(summaries_before, fit$posterior_summaries)) {
    sbt_fit_save(fit, file, overwrite = TRUE)
  }
  invisible(TRUE)
}
grid_posterior_summary_results <- if (
  .Platform$OS.type != "windows" &&
    grid_posterior_summary_cores > 1L
) {
  parallel::mclapply(
    grid_posterior_files,
    prepare_grid_posterior_summary,
    mc.cores = min(
      grid_posterior_summary_cores,
      length(grid_posterior_files)
    ),
    mc.preschedule = FALSE
  )
} else {
  lapply(grid_posterior_files, prepare_grid_posterior_summary)
}
grid_posterior_summary_ok <- vapply(
  grid_posterior_summary_results,
  isTRUE,
  logical(1)
)
if (any(!grid_posterior_summary_ok)) {
  stop(
    "Posterior-summary preparation failed for: ",
    paste(
      basename(grid_posterior_files[!grid_posterior_summary_ok]),
      collapse = ", "
    ),
    call. = FALSE
  )
}
grid_metadata_file <- file.path(grid_dir, "run_metadata.rds")
selected_grid_metadata <- if (file.exists(grid_metadata_file)) {
  read_rds_cache(grid_metadata_file)
} else {
  NULL
}
grid_diagnostics <- read_csv(
  file.path(grid_dir, "grid_mcmc_diagnostics.csv"),
  show_col_types = FALSE
)
required_grid_diagnostics <- c(
  "Cell", "h", "psi", "convergence", "objective", "estimable", "max_rhat",
  "min_bulk_ess", "min_tail_ess", "divergences", "chains", "warmup",
  "samples_per_chain", "metric",
  "max_treedepth_hits", "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", "run_signature", "file"
)
if (!all(required_grid_diagnostics %in% names(grid_diagnostics)) ||
    nrow(grid_diagnostics) != 9L || anyDuplicated(grid_diagnostics$Cell) ||
    !identical(sort(as.integer(grid_diagnostics$Cell)), seq_len(9L)) ||
    anyNA(grid_diagnostics$run_signature) ||
    any(!nzchar(grid_diagnostics$run_signature))) {
  stop("Grid diagnostics are not a complete nine-cell table.", call. = FALSE)
}
grid_diagnostics <- grid_diagnostics |>
  arrange(.data$Cell) |>
  mutate(
    file = vapply(
      .data$file,
      resolve_grid_file,
      character(1),
      selected_grid_dir = grid_dir
    )
  )
expected_grid_values <- grid_values |>
  mutate(Cell = row_number(), .before = 1)
if (!isTRUE(all.equal(
      as.numeric(grid_diagnostics$h), as.numeric(expected_grid_values$h),
      tolerance = 1e-12
    )) ||
    !isTRUE(all.equal(
      as.numeric(grid_diagnostics$psi), as.numeric(expected_grid_values$psi),
      tolerance = 1e-12
    )) ||
    !identical(cell_number(grid_diagnostics$file), seq_len(9L))) {
  stop("Grid diagnostic cells do not match the expected h/psi grid.", call. = FALSE)
}

grid_base_reuse_identity_equivalent <- function(saved, expected, source_fit) {
  if (!is.list(saved) || !is.list(expected) ||
      !esc31_saved_fit_signature_matches(
        saved$source_fit_scientific_signature,
        source_fit
      ) ||
      !esc31_saved_fit_signature_matches(
        expected$source_fit_scientific_signature,
        source_fit
      )) {
    return(FALSE)
  }
  stable_fields <- c(
    "schema_version", "cell", "grid_values",
    "source_mle_run_signature", "source_mcmc_run_signature",
    "source_posterior_payload_checksum", "source_sampler",
    "source_acceptance_thresholds", "implementation"
  )
  normalize_decision <- function(decision) {
    if (is.list(decision)) decision$source_posterior_file <- NULL
    decision
  }
  identical(saved[stable_fields], expected[stable_fields]) &&
    identical(
      normalize_decision(saved$decision),
      normalize_decision(expected$decision)
    )
}

grid_prerequisite_identity <- function(x) {
  if (!is.list(x)) return(NULL)
  x[c(
    "harvest_wall_contract", "harvest_wall_decision", "harvest_wall_data",
    "harvest_wall_recomputation_tolerance",
    "visual_review_required_before_computation", "implementation"
  )]
}

grid_specification_identity <- function(x) {
  if (!is.list(x)) return(NULL)
  if (is.list(x$base_cell_reuse)) {
    x$base_cell_reuse$source_posterior_file <- NULL
  }
  x
}

grid_signature_inputs_equivalent <- function(saved, expected, source_fit) {
  if (!is.list(saved) || !is.list(expected) ||
      !grid_base_reuse_identity_equivalent(
        saved$base_cell_reuse,
        expected$base_cell_reuse,
        source_fit
      ) ||
      !identical(
        grid_prerequisite_identity(saved$prerequisite_acceptance),
        grid_prerequisite_identity(expected$prerequisite_acceptance)
      )) {
    return(FALSE)
  }
  if (is.null(saved$failed_cell_rerun_decision) &&
      is.null(expected$failed_cell_rerun_decision)) {
    saved$failed_cell_rerun_decision <- NULL
    expected$failed_cell_rerun_decision <- NULL
  }
  storage_or_implementation_fields <- c(
    "base_fit_md5", "sbt_function_identity", "base_cell_reuse",
    "prerequisite_acceptance"
  )
  saved_stable <- saved[
    setdiff(names(saved), storage_or_implementation_fields)
  ]
  expected_stable <- expected[
    setdiff(names(expected), storage_or_implementation_fields)
  ]
  saved_stable$grid_specification <-
    grid_specification_identity(saved_stable$grid_specification)
  expected_stable$grid_specification <-
    grid_specification_identity(expected_stable$grid_specification)
  identical(saved_stable, expected_stable)
}

grid_metadata_signature_valid <-
  is.list(selected_grid_metadata) &&
  is.character(selected_grid_metadata$signature) &&
  length(selected_grid_metadata$signature) == 1L &&
  !is.na(selected_grid_metadata$signature) &&
  nzchar(selected_grid_metadata$signature) &&
  identical(
    selected_grid_metadata$signature,
    esc31_object_md5(selected_grid_metadata$signature_inputs)
  )
grid_is_mixed_selection <- isTRUE(
  grid_metadata_signature_valid &&
    identical(
      selected_grid_metadata$stage,
      "grid_mcmc_mixed_source_selection"
    )
)

if (grid_is_mixed_selection) {
  grid_mixed_artifacts <- selected_grid_metadata$artifacts
  required_mixed_artifact_fields <- c(
    "Cell", "source_run_id", "source_run_metadata_signature",
    "source_grid_run_signature", "retained_draws_per_chain",
    "source_file_md5", "posterior_payload_checksum", "state_checksum"
  )
  mixed_selection_valid <-
    identical(selected_grid_metadata$schema_version, 1L) &&
    identical(
      selected_grid_metadata$run_id,
      esc31_grid_mixed_source_selection$output_run_id
    ) &&
    identical(
      selected_grid_metadata$selection,
      esc31_grid_mixed_source_selection
    ) &&
    identical(
      selected_grid_metadata$signature_inputs$selection,
      esc31_grid_mixed_source_selection
    ) &&
    identical(
      selected_grid_metadata$signature_inputs$grid_specification,
      esc31_grid_specification
    ) &&
    identical(
      selected_grid_metadata$signature_inputs$failed_cell_rerun,
      esc31_grid_failed_cell_rerun
    ) &&
    is.data.frame(grid_mixed_artifacts) &&
    nrow(grid_mixed_artifacts) == 9L &&
    all(required_mixed_artifact_fields %in% names(grid_mixed_artifacts)) &&
    identical(as.integer(grid_mixed_artifacts$Cell), seq_len(9L)) &&
    identical(
      as.character(grid_mixed_artifacts$source_run_id),
      unname(
        esc31_grid_mixed_source_selection$source_run_by_cell[
          as.character(seq_len(9L))
        ]
      )
    ) &&
    identical(
      as.integer(grid_mixed_artifacts$retained_draws_per_chain),
      unname(
        esc31_grid_mixed_source_selection$retained_draws_per_chain[
          as.character(seq_len(9L))
        ]
      )
    ) &&
    identical(
      as.character(grid_mixed_artifacts$source_grid_run_signature),
      as.character(grid_diagnostics$run_signature)
    ) &&
    identical(
      as.character(grid_mixed_artifacts$posterior_payload_checksum),
      as.character(grid_diagnostics$posterior_payload_checksum)
    ) &&
    identical(
      as.character(grid_mixed_artifacts$state_checksum),
      as.character(grid_diagnostics$state_checksum)
    ) &&
    all(grepl(
      "^[0-9a-f]{32}$",
      as.character(grid_mixed_artifacts$source_file_md5)
    )) &&
    all(grepl(
      "^[0-9a-f]{32}$",
      as.character(
        grid_mixed_artifacts$source_run_metadata_signature
      )
    ))
  if (!isTRUE(mixed_selection_valid)) {
    stop(
      "The selected mixed-source grid manifest is incomplete or does not ",
      "match the approved per-cell selection.",
      call. = FALSE
    )
  }
  grid_expected_samples_per_chain <-
    as.integer(grid_mixed_artifacts$retained_draws_per_chain)
  grid_expected_run_signature <-
    as.character(grid_mixed_artifacts$source_grid_run_signature)
  grid_expected_source_run_id <-
    as.character(grid_mixed_artifacts$source_run_id)
  grid_expected_source_metadata_signature <-
    as.character(
      grid_mixed_artifacts$source_run_metadata_signature
    )
  expected_grid_run_settings <- NULL
  expected_grid_run_signature <- grid_expected_run_signature
} else {
  if (!grid_metadata_signature_valid ||
      !identical(selected_grid_metadata$schema_version, 9L) ||
      !identical(selected_grid_metadata$stage, "grid_mcmc") ||
      !grid_signature_inputs_equivalent(
        selected_grid_metadata$signature_inputs,
        grid_signature_inputs,
        base_grid_posterior_fit
      )) {
    stop(
      "The selected grid has no matching current run metadata.",
      call. = FALSE
    )
  }
  selected_grid_sampler_run_metadata <- selected_grid_metadata[c(
    "schema_version", "stage", "run_id", "signature_inputs", "signature"
  )]
  expected_grid_run_settings <- sbt:::.grid_mcmc_run_settings(
    control = grid_mcmc_config$control,
    n_passes = grid_mcmc_config$n_passes,
    b0_start_step = grid_mcmc_config$b0_start_step,
    b0_start_max = grid_mcmc_config$b0_start_max,
    check = grid_mcmc_config$check,
    stop_on_failure = grid_mcmc_config$stop_on_failure,
    mcmc_args = grid_mcmc_sampler_args,
    run_metadata = selected_grid_sampler_run_metadata
  )
  expected_grid_run_signature <- sbt:::.grid_mcmc_signature(
    expected_grid_run_settings
  )
  grid_expected_samples_per_chain <- rep(
    as.integer(grid_mcmc_config$num_samples),
    nrow(grid_values)
  )
  grid_expected_run_signature <- rep(
    expected_grid_run_signature,
    nrow(grid_values)
  )
  grid_expected_source_run_id <- rep(
    selected_grid_metadata$run_id,
    nrow(grid_values)
  )
  grid_expected_source_metadata_signature <- rep(
    selected_grid_metadata$signature,
    nrow(grid_values)
  )
}
grid_expected_draws <- as.integer(
  grid_mcmc_config$chains * grid_expected_samples_per_chain
)
embedded_diagnostic_fields <- c(
  "Cell", "h", "psi", "convergence", "objective", "estimable",
  "chains", "warmup", "samples_per_chain", "metric",
  "max_rhat", "min_bulk_ess", "min_tail_ess", "divergences",
  "max_treedepth_hits", "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", "run_signature"
)
diagnostic_values_match <- function(saved, reported) {
  all(vapply(embedded_diagnostic_fields, function(field) {
    left <- saved[[field]]
    right <- reported[[field]]
    if (is.numeric(left) || is.numeric(right)) {
      return(isTRUE(all.equal(
        as.numeric(left), as.numeric(right),
        tolerance = 1e-10, check.attributes = FALSE
      )))
    }
    identical(as.character(left), as.character(right))
  }, logical(1)))
}

validate_grid_cell <- function(cell) {
  tryCatch({
    expected_samples_cell <-
      grid_expected_samples_per_chain[[cell]]
    expected_draws_cell <- grid_expected_draws[[cell]]
    expected_run_signature_cell <-
      grid_expected_run_signature[[cell]]
    expected_source_run_id_cell <-
      grid_expected_source_run_id[[cell]]
    expected_source_metadata_signature_cell <-
      grid_expected_source_metadata_signature[[cell]]
    loaded <- sbt:::load_grid_mcmc_cell(
      grid_diagnostics$file[[cell]],
      grid_values = grid_values,
      fallback_cell = cell,
      verify_validation = FALSE
    )
    if (!inherits(loaded$fit, "sbt_fit")) {
      stop("legacy or non-portable grid cell", call. = FALSE)
    }
    saved_reuse <- loaded$fit$provenance$metadata$grid_reuse
    embedded_reuse <- loaded$fit$provenance$metadata$grid$reuse
    reuse_contract_valid <- if (identical(cell, grid_base_cell)) {
      grid_base_reuse_identity_equivalent(
        saved_reuse,
        base_grid_reuse_identity,
        base_grid_posterior_fit
      ) &&
        grid_base_reuse_identity_equivalent(
          embedded_reuse,
          base_grid_reuse_identity,
          base_grid_posterior_fit
        ) &&
        identical(
          sbt:::.grid_mcmc_posterior_payload_checksum(loaded$fit$mcmc),
          base_grid_reuse_identity$source_posterior_payload_checksum
        ) &&
        identical(
          loaded$fit$mcmc$settings$sampler,
          base_grid_reuse_identity$source_sampler
        ) &&
        identical(
          loaded$fit$mcmc$settings$acceptance_thresholds,
          base_grid_reuse_identity$source_acceptance_thresholds
        )
    } else {
      is.null(saved_reuse) && is.null(embedded_reuse)
    }
    if (!isTRUE(reuse_contract_valid)) {
      stop("grid base-cell reuse provenance mismatch", call. = FALSE)
    }
    loaded_wall_data <- lapply(
      names(grid_harvest_wall_identity$data),
      function(name) loaded$fit$data[[name]]
    )
    names(loaded_wall_data) <- names(grid_harvest_wall_identity$data)
    if (!identical(loaded_wall_data, grid_harvest_wall_identity$data) ||
        !identical(
          loaded$run_metadata$signature_inputs$harvest_wall_contract,
          grid_harvest_wall_identity$executable
        ) ||
        !identical(
          loaded$run_metadata$signature_inputs$harvest_wall_decision,
          grid_harvest_wall_identity$decision
        )) {
      stop("grid cell harvest-wall contract mismatch", call. = FALSE)
    }
    embedded_diagnostic <- loaded$diagnostic
    portable_mcmc <- loaded$fit$mcmc
    sample_dimensions <- dim(portable_mcmc$samples)
    actual_sampler_layout_valid <-
      is.array(portable_mcmc$samples) &&
      length(sample_dimensions) == 3L &&
      identical(as.integer(sample_dimensions[[1L]]), as.integer(
        grid_mcmc_config$num_warmup + expected_samples_cell
      )) &&
      identical(as.integer(sample_dimensions[[2L]]),
                as.integer(grid_mcmc_config$chains)) &&
      identical(as.integer(portable_mcmc$chains),
                as.integer(grid_mcmc_config$chains)) &&
      identical(as.integer(portable_mcmc$warmup),
                as.integer(grid_mcmc_config$num_warmup)) &&
      identical(as.integer(portable_mcmc$iter),
                as.integer(expected_samples_cell)) &&
      identical(as.character(portable_mcmc$metric),
                as.character(grid_mcmc_config$metric)) &&
      identical(as.integer(portable_mcmc$max_treedepth),
                as.integer(grid_mcmc_config$max_treedepth))
    if (!isTRUE(actual_sampler_layout_valid)) {
      stop("actual sampler layout does not match the grid contract",
           call. = FALSE)
    }
    saved_grid_run_settings <- loaded$fit$provenance$metadata$run$settings
    saved_sampler_arguments <- saved_grid_run_settings$sampler$arguments
    saved_sampler_contract_valid <-
      is.list(saved_sampler_arguments) &&
      identical(as.integer(saved_sampler_arguments$chains),
                as.integer(grid_mcmc_config$chains)) &&
      identical(as.integer(saved_sampler_arguments$num_warmup),
                as.integer(grid_mcmc_config$num_warmup)) &&
      identical(as.integer(saved_sampler_arguments$num_samples),
                as.integer(expected_samples_cell)) &&
      identical(as.character(saved_sampler_arguments$metric),
                as.character(grid_mcmc_config$metric)) &&
      identical(as.integer(saved_sampler_arguments$seed),
                as.integer(grid_mcmc_config$seed)) &&
      identical(
        as.numeric(saved_sampler_arguments$control$adapt_delta),
        as.numeric(grid_mcmc_config$adapt_delta)
      ) &&
      identical(
        as.integer(saved_sampler_arguments$control$max_treedepth),
        as.integer(grid_mcmc_config$max_treedepth)
      ) &&
      identical(as.character(saved_sampler_arguments$init),
                as.character(grid_mcmc_config$init))
    recalculated_grid_run_signature <- tryCatch(
      sbt:::.grid_mcmc_signature(saved_grid_run_settings),
      error = function(error) NA_character_
    )
    saved_source_metadata_valid <-
      is.list(loaded$run_metadata) &&
      identical(
        as.character(loaded$run_metadata$run_id),
        as.character(expected_source_run_id_cell)
      ) &&
      identical(
        as.character(loaded$run_metadata$signature),
        as.character(expected_source_metadata_signature_cell)
      ) &&
      identical(
        as.character(recalculated_grid_run_signature),
        as.character(expected_run_signature_cell)
      ) &&
      identical(
        as.character(loaded$run_signature),
        as.character(expected_run_signature_cell)
      )
    if (!isTRUE(saved_sampler_contract_valid) ||
        !isTRUE(saved_source_metadata_valid) ||
        (!grid_is_mixed_selection &&
          !identical(saved_grid_run_settings, expected_grid_run_settings))) {
      stop(
        "saved sampler settings or recomputed run signature do not match ",
        "the grid contract",
        call. = FALSE
      )
    }
    expected_draws <- sample_dimensions[2L] *
      (sample_dimensions[1L] - portable_mcmc$warmup)
    if (!identical(as.integer(expected_draws),
                   as.integer(expected_draws_cell))) {
      stop("posterior draw count does not match the grid contract", call. = FALSE)
    }

    prior_validation <- sbt_fit_validation(
      loaded$fit,
      scope = "mcmc",
      verify = FALSE
    )
    loaded$fit <- check_mcmc(
      loaded$fit,
      cores = grid_state_diagnostic_cores,
      stop_on_failure = FALSE,
      mode = grid_fit_validation_mode
    )
    validation_record <- sbt_fit_validation(
      loaded$fit,
      scope = "mcmc",
      verify = !identical(grid_fit_validation_mode, "skip"),
      require_pass = TRUE
    )
    if (is.null(prior_validation) ||
        identical(grid_fit_validation_mode, "full")) {
      sbt_fit_save(
        loaded$fit,
        grid_diagnostics$file[[cell]],
        overwrite = TRUE
      )
    }
    loaded$mcmc <- as_tmbfit(loaded$fit)
    recalculated_mle_state <-
      loaded$fit$fit$diagnostics$biological_state_mle
    recalculated_posterior_state <-
      loaded$fit$fit$diagnostics$biological_state_posterior
    mle_max_gradient <-
      loaded$fit$fit$diagnostics$mle$max_gradient
    if (is.null(validation_record) ||
        !esc31_biological_state_diagnostics_match(
          loaded$state_diagnostics,
          recalculated_posterior_state
        )) {
      stop("embedded fit validation or biological-state diagnostics differ",
           call. = FALSE)
    }

    recalculated_diagnostic <- sbt:::diagnose_grid_mcmc(
      mcmc = loaded$mcmc,
      cell = cell,
      grid_cell = loaded$grid_cell,
      opt = loaded$fit$fit$opt,
      start_fn = embedded_diagnostic$start_fn,
      b0_start_multiplier = embedded_diagnostic$b0_start_multiplier,
      nlminb_passes = embedded_diagnostic$nlminb_passes,
      estimability = loaded$fit$fit$estimability,
      file = grid_diagnostics$file[[cell]],
      max_treedepth = grid_mcmc_config$max_treedepth,
      run_signature = loaded$run_signature,
      state_diagnostics = recalculated_posterior_state
    )

    valid <- identical(as.integer(loaded$cell), cell) &&
      is.data.frame(embedded_diagnostic) &&
      nrow(embedded_diagnostic) == 1L &&
      all(embedded_diagnostic_fields %in% names(embedded_diagnostic)) &&
      isTRUE(all.equal(
        as.numeric(loaded$grid_cell$h), expected_grid_values$h[[cell]],
        tolerance = 1e-12
      )) &&
      isTRUE(all.equal(
        as.numeric(loaded$grid_cell$psi), expected_grid_values$psi[[cell]],
        tolerance = 1e-12
      )) &&
      identical(
        loaded$run_metadata$signature,
        expected_source_metadata_signature_cell
      ) &&
      identical(
        loaded$run_signature,
        as.character(grid_diagnostics$run_signature[[cell]])
      ) &&
      identical(
        loaded$run_signature,
        as.character(expected_run_signature_cell)
      ) &&
      esc31_biological_state_diagnostics_pass(
        recalculated_mle_state, 1L
      ) &&
      esc31_biological_state_diagnostics_pass(
        recalculated_posterior_state, expected_draws
      ) &&
      diagnostic_values_match(
        recalculated_diagnostic,
        embedded_diagnostic
      ) &&
      diagnostic_values_match(
        recalculated_diagnostic,
        grid_diagnostics[cell, , drop = FALSE]
      )
    list(
      valid = isTRUE(valid),
      error = if (isTRUE(valid)) "" else "metadata or diagnostics mismatch",
      mle_state = recalculated_mle_state,
      posterior_state = recalculated_posterior_state,
      diagnostic = recalculated_diagnostic,
      mle_max_gradient = mle_max_gradient
    )
  }, error = function(error) {
    list(valid = FALSE, error = conditionMessage(error))
  })
}

grid_cell_validation <- lapply(seq_len(9L), validate_grid_cell)
cell_metadata_valid <- vapply(
  grid_cell_validation, `[[`, logical(1), "valid"
)
if (!all(cell_metadata_valid)) {
  failed <- which(!cell_metadata_valid)
  reasons <- vapply(grid_cell_validation[failed], `[[`, character(1), "error")
  stop(
    "Grid-cell validation failed: ",
    paste0("cell ", failed, " (", reasons, ")", collapse = "; "), ".",
    call. = FALSE
  )
}

grid_mle_gradient_summary <- tibble(
  Cell = seq_len(nrow(grid_diagnostics)),
  max_gradient = vapply(
    grid_cell_validation, `[[`, numeric(1), "mle_max_gradient"
  ),
  file_md5 = unname(tools::md5sum(grid_diagnostics$file))
)
grid_mle_gradient_checksum <- esc31_object_md5(list(
  limit = grid_mle_gradient_limit,
  cells = grid_mle_gradient_summary
))
grid_diagnostics <- grid_diagnostics |>
  mutate(
    mle_max_gradient = grid_mle_gradient_summary$max_gradient,
    .after = objective
  )

grid_diagnostic_checks <- grid_diagnostics |>
  mutate(
    expected_samples_per_chain = grid_expected_samples_per_chain,
    expected_state_draws = grid_expected_draws,
    diagnostics_pass =
      !is.na(convergence) & convergence == 0 &
      !is.na(estimable) & estimable &
      is.finite(objective) &
      is.finite(mle_max_gradient) &
        mle_max_gradient <= grid_mle_gradient_limit &
      !is.na(chains) & chains == grid_mcmc_config$chains &
      !is.na(warmup) & warmup == grid_mcmc_config$num_warmup &
      !is.na(samples_per_chain) &
        samples_per_chain == expected_samples_per_chain &
      !is.na(metric) & metric == grid_mcmc_config$metric &
      is.finite(max_rhat) & max_rhat < rhat_acceptance_limit &
      is.finite(min_bulk_ess) & min_bulk_ess >= ess_acceptance_limit &
      is.finite(min_tail_ess) & min_tail_ess >= ess_acceptance_limit &
      !is.na(divergences) & divergences == 0 &
      !is.na(max_treedepth_hits) & max_treedepth_hits == 0 &
      !is.na(state_passes) & state_passes &
      !is.na(state_draws_expected) &
        state_draws_expected == expected_state_draws &
      !is.na(state_draws_evaluated) &
        state_draws_evaluated == state_draws_expected &
      !is.na(state_invalid_draws) & state_invalid_draws == 0L &
      !is.na(state_non_finite_draws) & state_non_finite_draws == 0L &
      !is.na(state_invalid_cells) & state_invalid_cells == 0L &
      is.finite(state_max_raw_harvest) &
      state_max_raw_harvest <=
        biological_state_contract()$hrate_limit +
        biological_state_contract()$hrate_tolerance &
      is.finite(state_min_number) &
      state_min_number > biological_state_contract()$number_tolerance &
      is.finite(state_max_harvest_penalty) &
      abs(state_max_harvest_penalty) <=
        biological_state_contract()$penalty_tolerance &
      !is.na(state_checksum) & grepl("^[0-9a-f]{32}$", state_checksum) &
      !is.na(posterior_payload_checksum) &
        grepl("^[0-9a-f]{32}$", posterior_payload_checksum)
  )

grid_status <- grid_diagnostics |>
  summarise(
    cells = n(),
    converged = sum(convergence == 0, na.rm = TRUE),
    estimable = sum(estimable, na.rm = TRUE),
    max_mle_gradient = max(mle_max_gradient, na.rm = TRUE),
    total_divergences = sum(divergences, na.rm = TRUE),
    total_treedepth_hits = sum(max_treedepth_hits, na.rm = TRUE),
    state_draws_expected = sum(state_draws_expected, na.rm = TRUE),
    state_draws_evaluated = sum(state_draws_evaluated, na.rm = TRUE),
    invalid_state_draws = sum(state_invalid_draws, na.rm = TRUE),
    invalid_state_cells = sum(state_invalid_cells, na.rm = TRUE),
    max_raw_harvest = max(state_max_raw_harvest, na.rm = TRUE),
    min_number = min(state_min_number, na.rm = TRUE),
    max_harvest_penalty = max(state_max_harvest_penalty, na.rm = TRUE),
    max_rhat = max(max_rhat, na.rm = TRUE),
    min_bulk_ess = min(min_bulk_ess, na.rm = TRUE),
    min_tail_ess = min(min_tail_ess, na.rm = TRUE)
  )

Show code

grid_status |>
  mutate(
    across(
      c(
        cells,
        converged,
        estimable,
        total_divergences,
        total_treedepth_hits,
        state_draws_expected,
        state_draws_evaluated,
        invalid_state_draws,
        invalid_state_cells
      ),
      ~ format_decimal(.x, digits = 0)
    ),
    max_mle_gradient = format_decimal(max_mle_gradient, digits = 6),
    max_raw_harvest = format_decimal(max_raw_harvest, digits = 4),
    min_number = format_decimal(min_number, digits = 3),
    max_harvest_penalty =
      format_decimal(max_harvest_penalty, digits = 3),
    max_rhat = format_decimal(max_rhat, digits = 4),
    across(
      c(min_bulk_ess, min_tail_ess),
      ~ format_decimal(.x, digits = 3)
    )
  ) |>
  replace_missing() |>
  kable(caption = paste0(
    "Overall status for completed grid directory ",
    basename(grid_dir),
    "."
  ))

Table 3: Overall status for completed grid directory esc31_grid_dense_mixed_accepted_20260723_failed_cells_1500.

cells	converged	estimable	max_mle_gradient	total_divergences	total_treedepth_hits	state_draws_expected	state_draws_evaluated	invalid_state_draws	invalid_state_cells	max_raw_harvest	min_number	max_harvest_penalty	max_rhat	min_bulk_ess	min_tail_ess
9	9	9	0.004886	0	0	45,000	45,000	0	0	0.8821	291.315	0.000	1.0095	990.747	996.028

For production inference, the acceptance limits used in Table 4 are maximum \(\widehat{R} < 1.01\), minimum bulk effective sample size (ESS) at least 400, and minimum tail ESS at least 400. Monnahan’s fisheries-specific good-practice guidance recommends \(\widehat{R} < 1.01\), bulk ESS \(> 400\) for every parameter, and no NUTS divergences (Monnahan 2024). These limits follow the rank-normalized \(\widehat{R}\) and bulk/tail ESS guidance of Vehtari et al. (2021), who also recommend checking rank plots and whether Monte Carlo standard errors are sufficiently small for the quantities of interest. They are screening limits, not proof of convergence. A red cell identifies any failed gate: unsuccessful or missing MLE convergence, a non-finite or greater-than-0.01 validated MLE maximum gradient, false or missing estimability, excessive \(\widehat{R}\), insufficient ESS, nonzero or missing divergences, or nonzero or missing maximum-treedepth hits. It also identifies incomplete biological-state coverage, any invalid or non-finite retained draw, raw combined seasonal harvest above 0.9, non-positive abundance, or a continuation penalty above the numerical zero tolerance. Every retained draw is checked once when its portable cell is completed. The embedded payload-bound validation checksum must match the fit, and the saved grid diagnostics must match both the cell and CSV diagnostics. The formatting is calculated from the selected run rather than from a hard-coded description; a column without red cells therefore means that every selected cell passed that gate. Every failed cell must be investigated and a complete approved-grid rerun before any cell is included in the combined production posterior.

Show code

grid_diagnostic_checks |>
  transmute(
    Cell,
    `Overall gate` = diagnostics_pass,
    h,
    psi,
    Source = if_else(
      Cell %in% esc31_grid_failed_cell_rerun$failed_source_cells,
      "Long rerun",
      "Initial pass"
    ),
    `Retained / chain` = samples_per_chain,
    `MLE convergence` = convergence,
    Estimable = estimable,
    Objective = objective,
    `Max MLE gradient` = mle_max_gradient,
    `Max Rhat` = max_rhat,
    `Min bulk ESS` = min_bulk_ess,
    `Min tail ESS` = min_tail_ess,
    Divergences = divergences,
    `Max treedepth hits` = max_treedepth_hits,
    `State pass` = state_passes,
    `State draws expected` = state_draws_expected,
    `State draws checked` = state_draws_evaluated,
    `Invalid state draws` = state_invalid_draws,
    `Non-finite state draws` = state_non_finite_draws,
    `Invalid state cells` = state_invalid_cells,
    `Max raw harvest` = state_max_raw_harvest,
    `Min abundance` = state_min_number,
    `Max continuation penalty` = state_max_harvest_penalty
  ) |>
  mutate(
    Cell = format_decimal(Cell, digits = 0),
    h = format_decimal(h, digits = 1),
    psi = format_decimal(psi, digits = 2),
    `Overall gate` = format_flagged_text(
      if_else(
        is.na(`Overall gate`),
        "Missing",
        if_else(`Overall gate`, "Pass", "Fail")
      ),
      is.na(`Overall gate`) | !`Overall gate`
    ),
    `Retained / chain` = format_diagnostic_cell(
      `Retained / chain`,
      is.na(`Retained / chain`) |
        `Retained / chain` != grid_expected_samples_per_chain,
      digits = 0
    ),
    `MLE convergence` = format_diagnostic_cell(
      `MLE convergence`,
      is.na(`MLE convergence`) | `MLE convergence` != 0,
      digits = 0
    ),
    Estimable = format_flagged_text(Estimable, is.na(Estimable) | !Estimable),
    Objective = format_diagnostic_cell(
      Objective,
      !is.finite(Objective),
      digits = 3
    ),
    `Max MLE gradient` = format_diagnostic_cell(
      `Max MLE gradient`,
      !is.finite(`Max MLE gradient`) |
        `Max MLE gradient` > grid_mle_gradient_limit,
      digits = 4
    ),
    `Max Rhat` = format_diagnostic_cell(
      `Max Rhat`,
      !is.finite(`Max Rhat`) | `Max Rhat` >= rhat_acceptance_limit,
      digits = 4
    ),
    `Min bulk ESS` = format_diagnostic_cell(
      `Min bulk ESS`,
      !is.finite(`Min bulk ESS`) | `Min bulk ESS` < ess_acceptance_limit,
      digits = 3
    ),
    `Min tail ESS` = format_diagnostic_cell(
      `Min tail ESS`,
      !is.finite(`Min tail ESS`) | `Min tail ESS` < ess_acceptance_limit,
      digits = 3
    ),
    Divergences = format_diagnostic_cell(
      Divergences,
      is.na(Divergences) | Divergences != 0,
      digits = 0
    ),
    `Max treedepth hits` = format_diagnostic_cell(
      `Max treedepth hits`,
      is.na(`Max treedepth hits`) | `Max treedepth hits` != 0,
      digits = 0
    ),
    `State pass` = format_flagged_text(
      `State pass`, is.na(`State pass`) | !`State pass`
    ),
    `State draws expected` = format_diagnostic_cell(
      `State draws expected`,
      is.na(`State draws expected`) |
        `State draws expected` != grid_expected_draws,
      digits = 0
    ),
    `State draws checked` = format_diagnostic_cell(
      `State draws checked`,
      is.na(`State draws checked`) |
        `State draws checked` != grid_expected_draws,
      digits = 0
    ),
    `Invalid state draws` = format_diagnostic_cell(
      `Invalid state draws`,
      is.na(`Invalid state draws`) | `Invalid state draws` != 0,
      digits = 0
    ),
    `Non-finite state draws` = format_diagnostic_cell(
      `Non-finite state draws`,
      is.na(`Non-finite state draws`) | `Non-finite state draws` != 0,
      digits = 0
    ),
    `Invalid state cells` = format_diagnostic_cell(
      `Invalid state cells`,
      is.na(`Invalid state cells`) | `Invalid state cells` != 0,
      digits = 0
    ),
    `Max raw harvest` = format_diagnostic_cell(
      `Max raw harvest`,
      !is.finite(`Max raw harvest`) |
        `Max raw harvest` > biological_state_contract()$hrate_limit +
          biological_state_contract()$hrate_tolerance,
      digits = 4
    ),
    `Min abundance` = format_diagnostic_cell(
      `Min abundance`,
      !is.finite(`Min abundance`) |
        `Min abundance` <= biological_state_contract()$number_tolerance,
      digits = 3
    ),
    `Max continuation penalty` = format_diagnostic_cell(
      `Max continuation penalty`,
      !is.finite(`Max continuation penalty`) |
        abs(`Max continuation penalty`) >
          biological_state_contract()$penalty_tolerance,
      digits = 3
    ),
    across(where(is.numeric), ~ format_decimal(.x, digits = 3))
  ) |>
  replace_missing() |>
  kable(
    escape = FALSE,
    caption = "Cell-level MLE and MCMC diagnostics for the completed grid."
  )

Table 4: Cell-level MLE and MCMC diagnostics for the completed grid.

Cell	Overall gate	h	psi	Source	Retained / chain	Estimable	Objective	Max MLE gradient	Max Rhat	Min bulk ESS	Min tail ESS	State pass	State draws expected	State draws checked	Max raw harvest	Min abundance
1	Pass	0.6	1.50	Long rerun	1,500	TRUE	7,617.816	0.0000	1.0063	2,195.011	1,739.950	TRUE	6,000	6,000	0.8692	549.878
2	Pass	0.7	1.50	Long rerun	1,500	TRUE	7,617.761	0.0003	1.0052	1,842.748	2,771.568	TRUE	6,000	6,000	0.8791	512.787
3	Pass	0.8	1.50	Initial pass	750	TRUE	7,619.536	0.0000	1.0095	1,002.989	1,001.219	TRUE	3,000	3,000	0.8678	340.952
4	Pass	0.6	1.75	Long rerun	1,500	TRUE	7,616.730	0.0000	1.0048	1,711.965	2,437.581	TRUE	6,000	6,000	0.8785	594.976
5	Pass	0.7	1.75	Initial pass	750	TRUE	7,616.638	0.0000	1.0087	1,157.424	1,178.248	TRUE	3,000	3,000	0.8787	490.556
6	Pass	0.8	1.75	Long rerun	1,500	TRUE	7,618.383	0.0049	1.0059	1,788.697	2,561.164	TRUE	6,000	6,000	0.8748	363.130
7	Pass	0.6	2.00	Long rerun	1,500	TRUE	7,617.612	0.0000	1.0068	1,811.713	1,703.064	TRUE	6,000	6,000	0.8699	411.923
8	Pass	0.7	2.00	Initial pass	750	TRUE	7,617.483	0.0001	1.0071	990.747	996.028	TRUE	3,000	3,000	0.8622	541.746
9	Pass	0.8	2.00	Long rerun	1,500	TRUE	7,619.193	0.0004	1.0050	2,351.868	3,009.794	TRUE	6,000	6,000	0.8821	291.315

MCMC-grid TRO trajectories

Figures Figure 1 and Figure 2 both show the accepted nine-cell MCMC grid. The first separates the cells into their \(h\)/\(\psi\) panels; the second overlays their posterior medians and identifies cell 5 as the accepted base posterior. Every median and interval comes from the validated embedded all-draw summary: 3,000 retained draws in cells 3, 5, and 8, and 6,000 retained draws in cells 1, 2, 4, 6, 7, and 9.

Show code

grid_biomass_summary <- bind_rows(lapply(
  seq_len(nrow(grid_diagnostics)),
  function(i) {
    row <- grid_diagnostics[i, ]
    fit <- sbt_fit_read(row$file, strict = TRUE, rebuild = FALSE)
    summary <- get_posterior_summary(
      fit,
      name = "plots",
      probs = c(0.025, 0.975),
      transform = "plots"
    )
    relative_tro <- summary$derived$relative_tro
    probability_index <- match(
      c(0.025, 0.5, 0.975),
      as.numeric(summary$probabilities)
    )
    n_years <- prod(summary$report_dimensions$spawning_biomass_y)
    if (!is.matrix(relative_tro) || anyNA(probability_index) ||
        nrow(relative_tro) != length(summary$probabilities) ||
        ncol(relative_tro) != n_years ||
        !identical(
          as.integer(summary$draws),
          as.integer(row$state_draws_expected)
        )) {
      stop(
        "Grid cell ", row$Cell,
        " lacks a complete all-draw relative-TRO posterior summary.",
        call. = FALSE
      )
    }
    tibble(
      Cell = as.integer(row$Cell),
      h = as.numeric(row$h),
      psi = as.numeric(row$psi),
      Year = seq.int(fit$data$first_yr, length.out = n_years),
      lower = as.numeric(relative_tro[probability_index[[1L]], ]),
      median = as.numeric(relative_tro[probability_index[[2L]], ]),
      upper = as.numeric(relative_tro[probability_index[[3L]], ]),
      Retained_draws = as.integer(summary$draws)
    )
  }
)) |>
  mutate(
    h_label = paste0("h = ", h),
    psi_label = paste0("psi = ", psi)
  )

Show code

ggplot(grid_biomass_summary, aes(x = Year, y = median)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = fit_expected_color,
    alpha = 0.18,
    linewidth = 0
  ) +
  geom_line(color = fit_expected_color, linewidth = 0.5) +
  facet_grid(psi_label ~ h_label) +
  labs(x = "Year", y = "Relative TRO") +
  scale_x_continuous(breaks = pretty_breaks(n = 6)) +
  scale_y_zero()

Figure 1: Relative total reproductive output trajectories for each of the nine accepted MCMC grid cells. Lines show posterior medians and shaded ribbons show 95% credible intervals calculated from every retained MCMC draw within each cell.

Show code

grid_biomass_overlay <- grid_biomass_summary |>
  mutate(
    `Grid cell` = factor(
      paste0(
        "Cell ", Cell,
        if_else(Cell == 5L, " (base)", ""),
        ": h=", h, ", psi=", psi
      ),
      levels = unique(paste0(
        "Cell ", Cell,
        if_else(Cell == 5L, " (base)", ""),
        ": h=", h, ", psi=", psi
      ))
    )
  )

grid_overlay_levels <- levels(grid_biomass_overlay$`Grid cell`)
grid_overlay_colors <- setNames(
  viridisLite::viridis(length(grid_overlay_levels), option = "D", end = 0.92),
  grid_overlay_levels
)
grid_overlay_colors[grepl("\\(base\\)", names(grid_overlay_colors))] <-
  fit_expected_color

ggplot(grid_biomass_overlay, aes(x = Year, y = median)) +
  geom_ribbon(
    data = filter(grid_biomass_overlay, Cell == 5L),
    aes(x = Year, ymin = lower, ymax = upper),
    inherit.aes = FALSE,
    fill = fit_expected_color,
    alpha = 0.18,
    linewidth = 0
  ) +
  geom_line(aes(color = `Grid cell`), linewidth = 0.7) +
  scale_color_manual(values = grid_overlay_colors) +
  labs(
    x = "Year",
    y = "Relative TRO",
    color = NULL
  ) +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero() +
  theme(legend.position = "bottom")

Figure 2: Posterior-median relative total reproductive output trajectories overlaid for all nine accepted MCMC grid cells. The all-draw shaded 95% credible interval is shown only for mid cell 5, the accepted base posterior; the other cells are shown as medians without uncertainty bands.

Combine grid MCMCs

The projection posterior is built here so downstream projection runs all use the same balanced 2,000-draw sample from the nine-cell MCMC grid. Draws are taken approximately evenly across cells; the two remainder draws are assigned randomly using a fixed seed.

Show code

failed_grid_cells <- grid_diagnostic_checks |>
  filter(!diagnostics_pass) |>
  pull(Cell)
if (length(failed_grid_cells)) {
  stop(
    "Grid cells ", paste(failed_grid_cells, collapse = ", "),
    " fail the production MCMC thresholds and cannot be combined. ",
    "Run a complete approved grid that passes every threshold.",
    call. = FALSE
  )
}

grid_selected_provenance_matches <- if (grid_is_mixed_selection) {
  identical(
    selected_grid_metadata$selection,
    esc31_grid_mixed_source_selection
  )
} else {
  identical(
    selected_grid_metadata$signature,
    grid_run_metadata$signature
  )
}
grid_provenance_matches <- !is.null(selected_grid_metadata) &&
  isTRUE(grid_selected_provenance_matches) &&
  is.list(selected_grid_metadata$prerequisite_review_confirmation) &&
  isTRUE(selected_grid_metadata$prerequisite_review_confirmation$confirmed) &&
  identical(
    selected_grid_metadata$prerequisite_review_confirmation$scope,
    "accepted_base_mcmc_only"
  ) &&
  identical(
    selected_grid_metadata$prerequisite_review_confirmation$mechanism,
    "esc31_base_only_grid_launch_v1"
  )
if (!grid_provenance_matches) {
  stop(
    "The selected grid has no matching base-data/model/workflow provenance. ",
    "Run a complete approved grid before combining.",
    call. = FALSE
  )
}

grid_files <- grid_diagnostic_checks |>
  arrange(Cell) |>
  pull(file)

combined_grid_state_summary <- grid_diagnostic_checks |>
  arrange(Cell) |>
  select(
    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
  )
combined_grid_state_checksum <- esc31_object_md5(list(
  contract = biological_state_contract(),
  diagnostics = combined_grid_state_summary
))

combined_run_metadata <- list(
  schema_version = 9L,
  stage = "combined_grid_posterior",
  grid_specification = esc31_grid_specification,
  prerequisite_acceptance = grid_prerequisite_validation,
  prerequisite_review_confirmation =
    selected_grid_metadata$prerequisite_review_confirmation,
  base_fit_md5 = base_fit_md5,
  base_run_signature = base_fit$provenance$metadata$run_identity$signature,
  grid_run_signature = selected_grid_metadata$signature %||% NA_character_,
  grid_mixed_source_selection = if (grid_is_mixed_selection) {
    selected_grid_metadata$selection
  } else {
    NULL
  },
  grid_source_artifacts = if (grid_is_mixed_selection) {
    selected_grid_metadata$artifacts
  } else {
    NULL
  },
  grid_sampler_run_signature = expected_grid_run_signature,
  grid_sampler_contract = list(
    chains = grid_mcmc_config$chains,
    num_warmup = grid_mcmc_config$num_warmup,
    num_samples_by_cell = grid_expected_samples_per_chain,
    metric = grid_mcmc_config$metric,
    adapt_delta = grid_mcmc_config$adapt_delta,
    max_treedepth = grid_mcmc_config$max_treedepth,
    init = grid_mcmc_config$init
  ),
  base_cell_reuse = base_grid_reuse_identity,
  grid_file_signature = grid_signature(grid_files),
  diagnostics_md5 = unname(tools::md5sum(file.path(grid_dir, "grid_mcmc_diagnostics.csv"))),
  mle_gradient_limit = grid_mle_gradient_limit,
  mle_gradient_summary = grid_mle_gradient_summary,
  mle_gradient_checksum = grid_mle_gradient_checksum,
  biological_state_contract = biological_state_contract(),
  harvest_wall_contract = grid_harvest_wall_identity$executable,
  harvest_wall_decision = grid_harvest_wall_identity$decision,
  harvest_wall_data = grid_harvest_wall_identity$data,
  harvest_wall_application = "fitted_grid_objective_conditioning",
  biological_state_draws_per_cell = grid_expected_draws,
  biological_state_summary = combined_grid_state_summary,
  biological_state_checksum = combined_grid_state_checksum,
  target_sample_names = target_sample_names,
  draws = combined_mcmc_draws,
  seed = combined_mcmc_seed,
  implementation =
    "grid_mcmc_to_tmbfit_v11_mixed_source_reuse_accepted_base_cell5_harvest_wall_bound"
)
combined_metadata_signature <- esc31_object_md5(combined_run_metadata)
combined_payload_checksum <- function(posterior, metadata, plan) {
  esc31_object_md5(list(
    samples = posterior$samples,
    warmup = posterior$warmup,
    sample_names = posterior$sample_names,
    draw_metadata = as.data.frame(metadata),
    draw_plan = as.data.frame(plan)
  ))
}
combined_provenance_signature <- function(metadata_signature,
                                          payload_checksum) {
  esc31_object_md5(list(
    metadata_signature = metadata_signature,
    payload_checksum = payload_checksum
  ))
}

use_saved_combined <- file.exists(combined_mcmc_file) && !rebuild_combined_mcmc
if (use_saved_combined) {
  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"
  )
  saved_combined_metadata <- if (!is.null(combined_saved) &&
      exists(
        "combined_run_metadata",
        envir = combined_saved,
        inherits = FALSE
      )) {
    combined_saved$combined_run_metadata
  } else {
    NULL
  }
  # The portable base file can gain derived summaries without changing its
  # accepted runtime/scientific identity. Keep its exact file hash for audit,
  # but do not invalidate the payload-bound combined posterior for that alone.
  substantive_metadata_fields <- setdiff(
    names(combined_run_metadata),
    c("prerequisite_acceptance", "base_fit_md5", "base_cell_reuse")
  )
  combined_base_reuse_matches <- is.list(saved_combined_metadata) &&
    grid_base_reuse_identity_equivalent(
      saved_combined_metadata$base_cell_reuse,
      combined_run_metadata$base_cell_reuse,
      base_grid_posterior_fit
    )
  combined_prerequisite_acceptance_equivalent <- function(saved, current) {
    if (!isTRUE(combined_base_reuse_matches) || !is.list(saved) ||
        !is.list(current) ||
        !identical(names(saved$artifact_checksums), "base") ||
        !identical(names(current$artifact_checksums), "base")) {
      return(FALSE)
    }
    # The base-cell reuse gate above binds the exact accepted posterior payload,
    # sampler, thresholds, and run identities. This lets the prerequisite's
    # whole-file checksum remain audit-only when derived summaries are added.
    saved$artifact_checksums <- current$artifact_checksums
    esc31_grid_prerequisite_acceptance_equivalent(saved, current)
  }
  combined_metadata_mismatch_fields <- if (
    is.list(saved_combined_metadata)
  ) {
    substantive_metadata_fields[!vapply(
      substantive_metadata_fields,
      function(field) identical(
        combined_run_metadata[[field]],
        saved_combined_metadata[[field]]
      ),
      logical(1)
    )]
  } else {
    substantive_metadata_fields
  }
  combined_reuse_checks <- c(
    readable_cache = !is.null(combined_saved),
    required_objects = !is.null(combined_saved) && all(vapply(
        required_combined_objects,
        exists,
        logical(1),
        envir = combined_saved,
        inherits = FALSE
      )),
    metadata_schema = is.list(saved_combined_metadata) &&
      setequal(names(saved_combined_metadata), names(combined_run_metadata)),
    base_cell_reuse = isTRUE(combined_base_reuse_matches),
    metadata_fields = !length(combined_metadata_mismatch_fields),
    prerequisite_acceptance =
      combined_prerequisite_acceptance_equivalent(
        saved_combined_metadata$prerequisite_acceptance,
        combined_run_metadata$prerequisite_acceptance
      ),
    posterior_container = is.list(combined_saved$projection_mcmc),
    samples_array = is.array(combined_saved$projection_mcmc$samples),
    samples_dimensions =
      length(dim(combined_saved$projection_mcmc$samples)) == 3L,
    draw_metadata = is.data.frame(combined_saved$draw_metadata),
    draw_plan = is.data.frame(combined_saved$draw_plan)
  )
  use_saved_combined <- all(combined_reuse_checks)
  if (use_saved_combined) {
    combined_run_metadata <- saved_combined_metadata
    combined_metadata_signature <- esc31_object_md5(combined_run_metadata)
    recalculated_payload_checksum <- tryCatch(
      combined_payload_checksum(
        combined_saved$projection_mcmc,
        combined_saved$draw_metadata,
        combined_saved$draw_plan
      ),
      error = function(error) NA_character_
    )
    use_saved_combined <- is.character(recalculated_payload_checksum) &&
      length(recalculated_payload_checksum) == 1L &&
      !is.na(recalculated_payload_checksum) &&
      identical(
        combined_saved$combined_payload_checksum,
        recalculated_payload_checksum
      ) &&
      identical(
        combined_saved$saved_source_signature,
        combined_provenance_signature(
          combined_metadata_signature,
          recalculated_payload_checksum
        )
      )
  }
  if (use_saved_combined) {
    projection_mcmc <- combined_saved$projection_mcmc
    draw_metadata <- combined_saved$draw_metadata
    draw_plan <- combined_saved$draw_plan
    combined_payload_checksum_value <-
      combined_saved$combined_payload_checksum
    saved_source_signature <- combined_saved$saved_source_signature
  }
}

if (!use_saved_combined && !run_grid_postprocessing) {
  if (exists("combined_reuse_checks", inherits = FALSE) &&
      any(!combined_reuse_checks)) {
    stop(
      "The saved combined-grid posterior failed reuse checks: ",
      paste(names(combined_reuse_checks)[!combined_reuse_checks], collapse = ", "),
      ".",
      call. = FALSE
    )
  }
  if (exists("combined_base_reuse_matches", inherits = FALSE) &&
      !isTRUE(combined_base_reuse_matches)) {
    stop(
      "The saved combined-grid base-cell reuse identity is incompatible.",
      call. = FALSE
    )
  }
  if (exists("combined_metadata_mismatch_fields", inherits = FALSE) &&
      length(combined_metadata_mismatch_fields)) {
    stop(
      "Combined-grid metadata fields requiring review: ",
      paste(combined_metadata_mismatch_fields, collapse = ", "), ".",
      call. = FALSE
    )
  }
  cat(
    "::: {.callout-warning}\n",
    "## Accepted grid; downstream artifacts pending\n\n",
    "All nine selected MCMC grid cells pass the production gates. The ",
    "compatible balanced 2,000-draw posterior and production direct-M ",
    "108-fit MLE grid are currently missing, so neither grid-derived management output ",
    "nor projection results are shown. Rebuild and check those artifacts ",
    "before continuing the production sequence.\n",
    ":::\n",
    sep = ""
  )
  knitr::knit_exit()
} else {
if (!use_saved_combined) {
  if (!run_grid_postprocessing) {
    stop(
      "The compatible combined-grid posterior is unavailable. Set ",
      "ESC31_RUN_GRID_POSTPROCESSING=true only after accepting the approved grid.",
      call. = FALSE
    )
  }
  combined_grid <- grid_mcmc_to_tmbfit(
    grid_files = grid_files,
    grid_values = grid_values,
    target_sample_names = target_sample_names,
    target_par_names = target_par_names,
    allow_mixed_run_signatures = grid_is_mixed_selection,
    n_draws = combined_mcmc_draws,
    seed = combined_mcmc_seed
  )

  projection_mcmc <- combined_grid$fit
  draw_metadata <- combined_grid$draw_metadata
  draw_plan <- combined_grid$draw_plan
  combined_payload_checksum_value <- combined_payload_checksum(
    projection_mcmc, draw_metadata, draw_plan
  )
  saved_source_signature <- combined_provenance_signature(
    combined_metadata_signature,
    combined_payload_checksum_value
  )
  combined_grid_dir <- file.path("runs", "grid_mcmc", basename(grid_dir))

  atomic_save_rdata(
    list(
      projection_mcmc = projection_mcmc,
      draw_metadata = draw_metadata,
      draw_plan = draw_plan,
      combined_grid_dir = combined_grid_dir,
      combined_run_metadata = combined_run_metadata,
      combined_payload_checksum = combined_payload_checksum_value,
      saved_source_signature = saved_source_signature
    ),
    combined_mcmc_file
  )
}
combined_source_signature <- saved_source_signature

if (sum(draw_plan$draws) != combined_mcmc_draws) {
  stop("Balanced grid draw plan does not sum to ", combined_mcmc_draws, call. = FALSE)
}
}

Show code

draw_plan |>
  mutate(
    h = format_decimal(h, digits = 1),
    psi = format_decimal(psi, digits = 2)
  ) |>
  replace_missing() |>
  kable(
    caption = paste0(
      "Balanced grid posterior draw plan. Draw counts sum to ",
      format_decimal(sum(draw_plan$draws), digits = 0),
      " draws."
    )
  )

Table 5: Balanced grid posterior draw plan. Draw counts sum to 2,000 draws.

cell	h	psi	draws	file
1	0.6	1.50	222	grid1.sbt.rds
2	0.7	1.50	222	grid2.sbt.rds
3	0.8	1.50	222	grid3.sbt.rds
4	0.6	1.75	222	grid4.sbt.rds
5	0.7	1.75	222	grid5.sbt.rds
6	0.8	1.75	223	grid6.sbt.rds
7	0.6	2.00	223	grid7.sbt.rds
8	0.7	2.00	222	grid8.sbt.rds
9	0.8	2.00	222	grid9.sbt.rds

Show code

combined_grid_post <- extract_samples(projection_mcmc)
if (!"par_log_m10" %in% names(combined_grid_post)) {
  stop("The balanced grid posterior lacks par_log_m10.", call. = FALSE)
}
m10_posterior <- exp(combined_grid_post$par_log_m10)
if ("par_log_m0" %in% names(combined_grid_post)) {
  m0_posterior <- exp(combined_grid_post$par_log_m0)
} else {
  if (!identical(as.integer(base_fit$data$M_switch), 2L)) {
    stop(
      "par_log_m0 is absent outside the length-based mortality model.",
      call. = FALSE
    )
  }
  ages <- seq.int(base_fit$data$min_age, base_fit$data$max_age)
  age0 <- match(0L, ages)
  age10 <- match(10L, ages)
  length_at_age <- base_fit$data$length_mu_ysa[1L, 1L, ]
  mc <- as.numeric(base_fit$parameters$par_mc)
  if (anyNA(c(age0, age10)) ||
      !all(is.finite(length_at_age[c(age0, age10)])) ||
      any(length_at_age[c(age0, age10)] <= 0) ||
      length(mc) != 1L || !is.finite(mc)) {
    stop(
      "Could not derive M0 from the length-based mortality contract.",
      call. = FALSE
    )
  }
  m0_posterior <- m10_posterior *
    (length_at_age[[age0]] / length_at_age[[age10]])^mc
}

Show code

mcmc_lev_df <- draw_metadata |>
  arrange(draw) |>
  transmute(
    h = h,
    M0 = m0_posterior[seq_len(n())],
    M10 = m10_posterior[seq_len(n())],
    psi = psi
  )

plot_levs(
  mcmc_lev_df,
  lev_vars,
  diag_fill = "#0072B2",
  point_color = "#1F1F1F",
  point_alpha = 0.18,
  point_size = 0.25,
  discrete_jitter = 0.42,
  discrete_levels = list(
    h = grid_values$h,
    psi = grid_values$psi
  )
)

Figure 3: Balanced combined MCMC grid posterior for h, M0, M10, and psi. Diagonal panels show marginal counts or histograms; off-diagonal panels show jittered posterior draws from the 2,000-draw combined grid posterior.

Show code

combined_grid_state <- function(n_draws) {
  combined_file <- file.path(
    projection_run_dir,
    paste0("esc31_projection_mcmc_", n_draws, ".rda")
  )
  if (!file.exists(combined_file)) {
    stop("Could not find combined grid posterior file: ", combined_file, call. = FALSE)
  }

  cache_file <- file.path(
    grid_dir,
    paste0("grid_combined_state_", n_draws, "_draws.rds")
  )
  cache_identity <- list(
    contract = "esc31_balanced_grid_state_v1",
    combined_file_md5 = unname(tools::md5sum(combined_file)),
    combined_payload_checksum = combined_payload_checksum_value,
    saved_source_signature = combined_source_signature,
    base_runtime_id = base_fit$runtime_id,
    model_implementation_signature = sbt:::.sbt_model_signature("sbt_model"),
    draws = n_draws
  )

  cached <- read_rds_cache(cache_file)
  required_columns <- c(
    "Draw", "Cell", "h", "psi", "Year",
    "relative_spawning_biomass", "recruitment", "B10_plus",
    "relative_B10_plus"
  )
  if (!is.list(cached) || !identical(cached$identity, cache_identity) ||
      !is.data.frame(cached$draws) ||
      length(setdiff(required_columns, names(cached$draws))) ||
      nrow(cached$draws) != n_draws * (base_fit$data$n_year + 1L) ||
      any(!is.finite(as.matrix(cached$draws[c(
        "relative_spawning_biomass", "recruitment", "B10_plus",
        "relative_B10_plus"
      )])))) {
    stop(
      "The compatible combined-grid state cache is unavailable. Run ",
      "`Rscript scripts/run-esc31-status-msy.R .` after installing the ",
      "matching sbt package.",
      call. = FALSE
    )
  }
  cached$draws
}

combined_biomass_2000 <- combined_grid_state(combined_mcmc_draws)

combined_biomass_2000_summary <- combined_biomass_2000 |>
  group_by(Year) |>
  summarise(
    lower = quantile(relative_spawning_biomass, 0.025, na.rm = TRUE),
    median = median(relative_spawning_biomass, na.rm = TRUE),
    upper = quantile(relative_spawning_biomass, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

combined_recruitment_2000_summary <- combined_biomass_2000 |>
  group_by(Year) |>
  summarise(
    lower = quantile(recruitment / 1e6, 0.025, na.rm = TRUE),
    median = median(recruitment / 1e6, na.rm = TRUE),
    upper = quantile(recruitment / 1e6, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

Show code

set.seed(42)
combined_biomass_2000 |>
  filter(Draw %in% sample(unique(Draw), min(200L, n_distinct(Draw)))) |>
  mutate(Cell = factor(Cell)) |>
  ggplot(aes(x = Year, y = relative_spawning_biomass, group = Draw, color = Cell)) +
  geom_line(alpha = 0.25, linewidth = 0.25) +
  labs(x = "Year", y = "Relative TRO", color = "Grid cell") +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero()

Figure 4: Sample of relative total reproductive output trajectories from the balanced 2,000-draw projection posterior, colored by source grid cell.

Show code

ggplot(combined_biomass_2000_summary, aes(x = Year, y = median)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = fit_expected_color,
    alpha = 0.18,
    linewidth = 0
  ) +
  geom_line(color = fit_expected_color, linewidth = 0.75) +
  labs(x = "Year", y = "Relative TRO") +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero()

Figure 5: Relative total reproductive output trajectory for the combined 2,000-draw grid posterior. The line shows the posterior median and the shaded ribbon shows the 95% interval across draws from all nine grid cells.

Show code

ggplot(combined_recruitment_2000_summary, aes(x = Year, y = median)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = fit_expected_color,
    alpha = 0.18,
    linewidth = 0
  ) +
  geom_line(color = fit_expected_color, linewidth = 0.75) +
  labs(x = "Year", y = "Recruitment (millions)") +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero(labels = label_number(accuracy = 1))

Figure 6: Recruitment trajectory for the balanced combined 2,000-draw MCMC grid posterior. The line shows the posterior median and the shaded ribbon shows the equal-tailed 95% credible interval across draws from all nine grid cells.

Combined-grid Kobe status

The Kobe plot combines the posterior ratios of total reproductive output to the MSY level and fishing mortality to the MSY level. It therefore preserves the paired uncertainty in stock and fishing status, rather than comparing separate marginal summaries. The 2023 and 2025 estimates use the same 2,000 balanced draws from the nine-cell MCMC grid.

Show code

kobe_cache_file <- file.path(
  run_dir, "msy", "grid_posterior_msy_2023_2025_2000.rds"
)
if (!file.exists(kobe_cache_file)) {
  stop(
    "Could not find the accepted combined-grid MSY cache: ",
    kobe_cache_file,
    call. = FALSE
  )
}

kobe_cache <- readRDS(kobe_cache_file)
if (!is.list(kobe_cache) ||
    !identical(
      kobe_cache$identity$contract,
      "esc31_balanced_grid_msy_2023_2025_v1"
    ) ||
    !is.data.frame(kobe_cache$result$summary)) {
  stop("The combined-grid MSY cache does not satisfy its identity contract.", call. = FALSE)
}

kobe_status <- kobe_cache$result$summary |>
  filter(valid, year %in% c(2023L, 2025L)) |>
  transmute(
    Draw = as.integer(draw),
    Year = factor(year, levels = c(2023L, 2025L)),
    TRO_TROMSY = B_Bmsy,
    F_FMSY = F_Fmsy
  )

if (nrow(kobe_status) != 2L * combined_mcmc_draws ||
    any(!is.finite(kobe_status$TRO_TROMSY)) ||
    any(!is.finite(kobe_status$F_FMSY)) ||
    any(table(kobe_status$Draw) != 2L)) {
  stop("The accepted Kobe cache does not contain two valid states per draw.", call. = FALSE)
}

kobe_summary <- kobe_status |>
  group_by(Year) |>
  summarise(
    x_lower = quantile(TRO_TROMSY, 0.025),
    x_median = median(TRO_TROMSY),
    x_upper = quantile(TRO_TROMSY, 0.975),
    y_lower = quantile(F_FMSY, 0.025),
    y_median = median(F_FMSY),
    y_upper = quantile(F_FMSY, 0.975),
    .groups = "drop"
  )

kobe_pairs <- kobe_status |>
  mutate(Year = as.character(Year)) |>
  pivot_wider(
    names_from = Year,
    values_from = c(TRO_TROMSY, F_FMSY),
    names_sep = "_"
  )

set.seed(73015)
kobe_pairs_sample <- kobe_pairs |>
  slice_sample(n = min(350L, nrow(kobe_pairs)))

kobe_quadrants_2025 <- kobe_status |>
  filter(Year == "2025") |>
  summarise(
    yellow = mean(TRO_TROMSY < 1 & F_FMSY <= 1),
    green = mean(TRO_TROMSY >= 1 & F_FMSY <= 1),
    red = mean(TRO_TROMSY < 1 & F_FMSY > 1),
    orange = mean(TRO_TROMSY >= 1 & F_FMSY > 1)
  )

kobe_hdr_contours <- function(data, probabilities = c(0.50, 0.80, 0.95), n = 180L) {
  density <- MASS::kde2d(
    data$TRO_TROMSY,
    data$F_FMSY,
    n = n,
    lims = c(0, 2, 0, 1.2)
  )
  dx <- diff(density$x[1:2])
  dy <- diff(density$y[1:2])
  ordered_density <- sort(as.vector(density$z), decreasing = TRUE)
  cumulative_probability <- cumsum(ordered_density * dx * dy)
  levels <- vapply(
    probabilities,
    function(probability) {
      ordered_density[which(cumulative_probability >= probability)[1L]]
    },
    numeric(1)
  )

  purrr::map2_dfr(levels, probabilities, function(level, probability) {
    lines <- contourLines(density$x, density$y, density$z, levels = level)
    purrr::imap_dfr(lines, function(line, piece) {
      tibble(
        x = line$x,
        y = line$y,
        coverage = factor(
          percent(probability, accuracy = 1),
          levels = c("50%", "80%", "95%")
        ),
        piece = paste0(percent(probability, accuracy = 1), "-", piece)
      )
    })
  })
}

kobe_contours_2025 <- kobe_hdr_contours(
  filter(kobe_status, Year == "2025")
)

Show code

kobe_year_colours <- c("2023" = "#4D4D4D", "2025" = "#005B96")
kobe_year_shapes <- c("2023" = 21, "2025" = 24)

ggplot() +
  annotate("rect", xmin = 0, xmax = 1, ymin = 0, ymax = 1,
           fill = "#F4D35E", alpha = 0.32) +
  annotate("rect", xmin = 1, xmax = 2, ymin = 0, ymax = 1,
           fill = "#5CB85C", alpha = 0.28) +
  annotate("rect", xmin = 0, xmax = 1, ymin = 1, ymax = 1.2,
           fill = "#D9534F", alpha = 0.34) +
  annotate("rect", xmin = 1, xmax = 2, ymin = 1, ymax = 1.2,
           fill = "#F0AD4E", alpha = 0.32) +
  geom_segment(
    data = kobe_pairs_sample,
    aes(
      x = TRO_TROMSY_2023,
      y = F_FMSY_2023,
      xend = TRO_TROMSY_2025,
      yend = F_FMSY_2025
    ),
    colour = "#3A506B",
    linewidth = 0.25,
    alpha = 0.10
  ) +
  geom_point(
    data = kobe_status,
    aes(
      x = TRO_TROMSY,
      y = F_FMSY,
      colour = Year,
      fill = Year,
      shape = Year
    ),
    size = 1.05,
    stroke = 0.25,
    alpha = 0.14
  ) +
  geom_path(
    data = kobe_contours_2025,
    aes(x = x, y = y, group = piece, linetype = coverage),
    colour = kobe_year_colours[["2025"]],
    linewidth = 0.75,
    alpha = 0.92
  ) +
  geom_segment(
    data = kobe_summary,
    aes(x = x_lower, xend = x_upper, y = y_median, yend = y_median, colour = Year),
    linewidth = 1.05,
    show.legend = FALSE
  ) +
  geom_segment(
    data = kobe_summary,
    aes(x = x_median, xend = x_median, y = y_lower, yend = y_upper, colour = Year),
    linewidth = 1.05,
    show.legend = FALSE
  ) +
  geom_segment(
    data = tibble(
      x = kobe_summary$x_median[[1L]],
      y = kobe_summary$y_median[[1L]],
      xend = kobe_summary$x_median[[2L]],
      yend = kobe_summary$y_median[[2L]]
    ),
    aes(x = x, y = y, xend = xend, yend = yend),
    colour = "#172A3A",
    linewidth = 1.15,
    arrow = arrow(length = unit(0.18, "cm"), type = "closed")
  ) +
  geom_point(
    data = kobe_summary,
    aes(x = x_median, y = y_median, colour = Year, fill = Year, shape = Year),
    size = 4.3,
    stroke = 1.1
  ) +
  geom_text(
    data = kobe_summary,
    aes(x = x_median, y = y_median, label = Year, colour = Year),
    nudge_x = c(-0.055, 0.055),
    nudge_y = c(0.040, -0.040),
    fontface = "bold",
    size = 3.6,
    show.legend = FALSE
  ) +
  annotate("text", x = 0.08, y = 1.15,
           label = "TRO < TRO[MSY]~~'|'~~F > F[MSY]", parse = TRUE,
           hjust = 0, colour = "#6B1E1E", size = 3.4) +
  annotate("text", x = 1.92, y = 1.15,
           label = "TRO >= TRO[MSY]~~'|'~~F > F[MSY]", parse = TRUE,
           hjust = 1, colour = "#7A4700", size = 3.4) +
  annotate("text", x = 0.08, y = 0.08,
           label = paste0(percent(kobe_quadrants_2025$yellow, accuracy = 0.1),
                          " of 2025 draws"),
           hjust = 0, colour = "#5F5200", fontface = "bold", size = 3.5) +
  annotate("text", x = 1.92, y = 0.08,
           label = paste0(percent(kobe_quadrants_2025$green, accuracy = 0.1),
                          " of 2025 draws"),
           hjust = 1, colour = "#175D2A", fontface = "bold", size = 3.5) +
  geom_vline(xintercept = 1, colour = "white", linewidth = 1.05) +
  geom_hline(yintercept = 1, colour = "white", linewidth = 1.05) +
  geom_vline(xintercept = 1, colour = "#272727", linewidth = 0.55) +
  geom_hline(yintercept = 1, colour = "#272727", linewidth = 0.55) +
  scale_colour_manual(values = kobe_year_colours, name = "Posterior year") +
  scale_fill_manual(values = kobe_year_colours, name = "Posterior year") +
  scale_shape_manual(values = kobe_year_shapes, name = "Posterior year") +
  scale_linetype_manual(
    values = c("50%" = "solid", "80%" = "22", "95%" = "42"),
    name = "2025 joint region"
  ) +
  scale_x_continuous(
    limits = c(0, 2),
    breaks = seq(0, 2, by = 0.25),
    expand = expansion(mult = 0)
  ) +
  scale_y_continuous(
    limits = c(0, 1.2),
    breaks = seq(0, 1.2, by = 0.2),
    expand = expansion(mult = 0)
  ) +
  labs(
    x = expression(TRO / TRO[MSY]),
    y = expression(F / F[MSY])
  ) +
  coord_cartesian(clip = "off") +
  guides(
    colour = guide_legend(order = 1, override.aes = list(alpha = 1, size = 3)),
    fill = guide_legend(order = 1, override.aes = list(alpha = 1, size = 3)),
    shape = guide_legend(order = 1, override.aes = list(alpha = 1, size = 3)),
    linetype = guide_legend(order = 2)
  ) +
  theme_bw(base_size = 12) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(colour = "white", linewidth = 0.35, alpha = 0.75),
    legend.position = "bottom",
    legend.box = "horizontal",
    aspect.ratio = 0.63
  )

Figure 7: Kobe status plot for the balanced combined MCMC grid. Small points are all 2,000 paired posterior draws in each year; faint paths show a reproducible sample of paired movement from 2023 to 2025. Large symbols are posterior medians, horizontal and vertical bars are marginal equal-tailed 95% credible intervals, and the arrow connects the two medians. Contours enclose the 50%, 80%, and 95% highest-density regions of the 2025 joint posterior.

The 2025 posterior is split between the lower-left and lower-right Kobe quadrants: 51.9% of draws have \(TRO < TRO_{MSY}\), while 48.1% have \(TRO \geq TRO_{MSY}\). All 2,000 draws remain below \(F_{MSY}\).

FAO annual Kobe trajectory

The joint-posterior Kobe plot above is retained as the primary assessment diagnostic. For continuity with the annual plot supplied for FAO reporting, the figure below follows the full fitted trajectory from 1952 to 2025. The black line joins annual posterior medians, and the grey crosshairs show the 25th to 75th percentiles for each annual stock- and fishing-status ratio.

Show code

fao_kobe_file <- file.path(
  esc_dir, "report_data", "grid_kobe_annual_2026.csv"
)
fao_kobe_provenance_file <- file.path(
  esc_dir, "report_data", "grid_kobe_annual_2026_provenance.csv"
)
if (!file.exists(fao_kobe_file) || !file.exists(fao_kobe_provenance_file)) {
  stop(
    "The committed annual Kobe data or its provenance record is missing. ",
    "Run scripts/run-esc31-grid-kobe-msy.R from the repository root.",
    call. = FALSE
  )
}

fao_kobe <- read_csv(fao_kobe_file, show_col_types = FALSE)
fao_kobe_provenance <- read_csv(
  fao_kobe_provenance_file,
  show_col_types = FALSE
)
required_fao_kobe_fields <- c(
  "Year", "TRO_lower", "TRO_median", "TRO_upper",
  "F_lower", "F_median", "F_upper", "Draws"
)
if (!identical(names(fao_kobe), required_fao_kobe_fields) ||
    nrow(fao_kobe) != 74L ||
    !identical(as.integer(fao_kobe$Year), 1952:2025) ||
    any(fao_kobe$Draws != combined_mcmc_draws) ||
    any(!is.finite(as.matrix(fao_kobe[setdiff(
      names(fao_kobe), c("Year", "Draws")
    )]))) ||
    nrow(fao_kobe_provenance) != 1L ||
    !identical(
      fao_kobe_provenance$contract,
      "esc31_balanced_grid_kobe_plot_data_v1"
    ) ||
    fao_kobe_provenance$annual_csv_md5 != unname(tools::md5sum(fao_kobe_file)) ||
    fao_kobe_provenance$valid_draw_years !=
      combined_mcmc_draws * nrow(fao_kobe) ||
    fao_kobe_provenance$invalid_draw_years != 0L ||
    fao_kobe_provenance$phase1_failures != 0L ||
    fao_kobe_provenance$final_failures != 0L) {
  stop("The annual Kobe plotting data failed its identity or validity gate.",
       call. = FALSE)
}

fao_kobe_x_max <- max(
  1.25,
  ceiling(max(fao_kobe$TRO_upper) * 4) / 4
)
fao_kobe_y_max <- max(
  1.25,
  ceiling(max(fao_kobe$F_upper) * 4) / 4
)
fao_kobe_labels <- fao_kobe |>
  filter(Year %in% c(1952L, seq(1960L, 2020L, by = 10L), 2025L)) |>
  left_join(
    tribble(
      ~Year, ~label_x, ~label_y,
      1952L, 3.50, 0.12,
      1960L, 2.28, 1.16,
      1970L, 2.31, 0.84,
      1980L, 1.78, 1.56,
      1990L, 0.63, 1.44,
      2000L, 0.53, 1.31,
      2010L, 0.27, 0.65,
      2020L, 0.96, 0.45,
      2025L, 0.84, 0.66
    ),
    by = "Year"
  ) |>
  mutate(
    label_face = if_else(Year %in% c(1952L, 2025L), "bold", "plain")
  )
if (anyNA(fao_kobe_labels[c("label_x", "label_y")])) {
  stop("The annual Kobe label-position registry is incomplete.", call. = FALSE)
}

fao_kobe_plot <- ggplot(fao_kobe, aes(x = TRO_median, y = F_median)) +
  annotate(
    "rect", xmin = 0, xmax = 1, ymin = 0, ymax = 1,
    fill = "#F4D35E", alpha = 0.42
  ) +
  annotate(
    "rect", xmin = 1, xmax = fao_kobe_x_max, ymin = 0, ymax = 1,
    fill = "#5CB85C", alpha = 0.36
  ) +
  annotate(
    "rect", xmin = 0, xmax = 1, ymin = 1, ymax = fao_kobe_y_max,
    fill = "#D9534F", alpha = 0.40
  ) +
  annotate(
    "rect", xmin = 1, xmax = fao_kobe_x_max,
    ymin = 1, ymax = fao_kobe_y_max,
    fill = "#F0AD4E", alpha = 0.38
  ) +
  geom_vline(xintercept = 1, colour = "white", linewidth = 1.05) +
  geom_hline(yintercept = 1, colour = "white", linewidth = 1.05) +
  geom_vline(xintercept = 1, colour = "black", linewidth = 0.55) +
  geom_hline(yintercept = 1, colour = "black", linewidth = 0.55) +
  geom_segment(
    aes(x = TRO_lower, xend = TRO_upper, yend = F_median),
    colour = "grey52",
    linewidth = 0.42,
    alpha = 0.72
  ) +
  geom_segment(
    aes(xend = TRO_median, y = F_lower, yend = F_upper),
    colour = "grey52",
    linewidth = 0.42,
    alpha = 0.72
  ) +
  geom_path(
    colour = "black",
    linewidth = 0.72,
    lineend = "round",
    linejoin = "round"
  ) +
  geom_point(colour = "black", size = 1.35) +
  geom_point(
    data = filter(fao_kobe, Year == max(Year)),
    shape = 21,
    fill = "#005B96",
    colour = "white",
    stroke = 1.05,
    size = 4.3
  ) +
  geom_segment(
    data = fao_kobe_labels,
    aes(
      x = TRO_median, y = F_median,
      xend = label_x, yend = label_y
    ),
    inherit.aes = FALSE,
    colour = "grey25",
    linewidth = 0.3
  ) +
  geom_label(
    data = fao_kobe_labels,
    aes(x = label_x, y = label_y, label = Year, fontface = label_face),
    inherit.aes = FALSE,
    colour = "black",
    fill = alpha("white", 0.82),
    linewidth = 0,
    label.padding = unit(0.08, "lines"),
    size = 3
  ) +
  scale_x_continuous(
    limits = c(0, fao_kobe_x_max),
    breaks = pretty_breaks(n = 7),
    expand = expansion(mult = 0)
  ) +
  scale_y_continuous(
    limits = c(0, fao_kobe_y_max),
    breaks = pretty_breaks(n = 7),
    expand = expansion(mult = 0)
  ) +
  labs(
    x = expression(TRO / TRO[MSY]),
    y = expression(F / F[MSY]~"(ages 2–15)")
  ) +
  coord_cartesian(
    xlim = c(0, fao_kobe_x_max),
    ylim = c(0, fao_kobe_y_max),
    expand = FALSE,
    clip = "off"
  ) +
  theme_bw(base_size = 12) +
  theme(
    panel.grid = element_blank(),
    aspect.ratio = 1,
    plot.margin = margin(8, 12, 8, 8)
  )

ggsave(
  filename = file.path(esc_dir, "report_data", "Kobe_2026.png"),
  plot = fao_kobe_plot,
  width = 8.4,
  height = 8.0,
  units = "in",
  dpi = 320,
  bg = "white"
)

Show code

fao_kobe_plot

Figure 8: Annual Kobe trajectory for the balanced combined MCMC grid, in the established FAO reporting style. The black line connects annual posterior medians from 1952 to 2025; grey horizontal and vertical bars show the interquartile ranges across the 2,000 balanced posterior draws. The terminal 2025 estimate is highlighted in blue. Fishing mortality is biomass-weighted over ages 2–15.

Sensitivity and MCMC-grid comparison

The accepted sensitivity runs and all nine accepted MCMC grid cells are shown together on one set of axes. Each trajectory is a posterior median; color distinguishes the two sources without assigning a separate color to every model. The terminal-year table retains each model’s 95% posterior interval.

Show code

sensitivity_grid_comparison_files <- c(
  base = model_file,
  no_uam = file.path(
    run_dir, "sens", "esc31_sens_no_uam.sbt.rds"
  ),
  drop_5yrs = file.path(
    run_dir, "sens", "esc31_sens_drop_5yrs.sbt.rds"
  ),
  cpue_omega_075 = file.path(
    run_dir, "sens", "esc31_sens_cpue_omega_075.sbt.rds"
  ),
  q2008 = file.path(
    run_dir, "sens", "esc31_sens_q2008.sbt.rds"
  ),
  ll1_terminal_3yr = file.path(
    run_dir, "sens", "esc31_sens_ll1_terminal_3yr.sbt.rds"
  ),
  no_pop_hsp = file.path(
    run_dir, "sens", "esc31_sens_no_pop_hsp.sbt.rds"
  ),
  no_hsp = file.path(
    run_dir, "sens", "esc31_sens_no_hsp.sbt.rds"
  ),
  troll = file.path(
    run_dir, "sens", "esc31_sens_troll.sbt.rds"
  ),
  cpue_ll1_sel = file.path(
    run_dir, "sens", "esc31_sens_cpue_ll1_sel.sbt.rds"
  ),
  constant_cpue_cv = file.path(
    run_dir, "sens", "esc31_sens_constant_cpue_cv.sbt.rds"
  ),
  length_m = file.path(
    run_dir, "sens", "esc31_sens_length_m.sbt.rds"
  ),
  indo_sel = file.path(
    run_dir, "sens", "esc31_sens_indo_sel.sbt.rds"
  ),
  estimate_m_slope = file.path(
    run_dir, "sens", "esc31_sens_estimate_m_slope.sbt.rds"
  )
)
sensitivity_grid_comparison_labels <- c(
  base = "Base",
  no_uam = "No UAM",
  drop_5yrs = "Drop CPUE 2021-2025",
  cpue_omega_075 = "CPUE power 0.75",
  q2008 = "CPUE q split from 2008",
  ll1_terminal_3yr = "Annual LL1 selectivity 2023-2025",
  no_pop_hsp = "No POP or HSP",
  no_hsp = "No HSP",
  troll = "GTI with base aerial tau 0.59",
  cpue_ll1_sel = "CPUE index with LL1 selectivity",
  constant_cpue_cv = "Constant CPUE CV 0.2",
  length_m = "Direct four-anchor M",
  indo_sel = "More flexible Indonesian selectivity",
  estimate_m_slope = "Estimated M-at-age slope"
)
if (!identical(
      names(sensitivity_grid_comparison_files),
      names(sensitivity_grid_comparison_labels)
    ) ||
    any(!file.exists(sensitivity_grid_comparison_files))) {
  stop(
    "The complete accepted sensitivity posterior registry is unavailable.",
    call. = FALSE
  )
}

sensitivity_grid_comparison_fits <- lapply(
  sensitivity_grid_comparison_files,
  sbt_fit_read,
  strict = TRUE,
  rebuild = FALSE
)
sensitivity_posterior_accepted <- function(fit) {
  diagnostics <- fit$fit$diagnostics$mcmc
  if (is.data.frame(diagnostics) && nrow(diagnostics) == 1L &&
      isTRUE(diagnostics$passes)) {
    return(TRUE)
  }
  decision <- fit$provenance$metadata$mcmc_acceptance
  is.list(decision) &&
    identical(
      decision$status,
      "accepted_with_explicit_divergence_exception"
    ) &&
    identical(
      as.integer(diagnostics$divergences),
      as.integer(decision$accepted_divergences)
    ) &&
    identical(
      diagnostics$posterior_payload_checksum,
      decision$posterior_payload_checksum
    ) &&
    identical(diagnostics$state_checksum, decision$state_checksum) &&
    isTRUE(diagnostics$state_passes)
}
if (!all(vapply(
      sensitivity_grid_comparison_fits,
      sensitivity_posterior_accepted,
      logical(1)
    ))) {
  stop(
    "At least one sensitivity posterior lacks a current acceptance decision.",
    call. = FALSE
  )
}

grid_comparison_fits <- lapply(
  unname(grid_files),
  sbt_fit_read,
  strict = TRUE,
  rebuild = FALSE
)
if (length(grid_comparison_fits) != 9L ||
    !all(vapply(
      grid_comparison_fits,
      function(fit) isTRUE(fit$fit$diagnostics$mcmc$passes),
      logical(1)
    ))) {
  stop("All nine accepted MCMC grid fits are required.", call. = FALSE)
}

posterior_tro_comparison_frame <- function(fit, model, source) {
  summary <- get_posterior_summary(
    fit,
    reports = c(
      "lf_pred", "cpue_lf_pred", "af_pred", "cpue_pred",
      "troll_pred", "aerial_pred", "catch_pred_ysf", "tag_pred",
      "M_a", "B0", "spawning_biomass_y", "recruitment_y"
    ),
    probs = c(0.025, 0.975),
    transform = "plots"
  )
  relative_tro <- summary$derived$relative_tro
  required_probabilities <- c("0.025", "0.500", "0.975")
  if (!is.matrix(relative_tro) ||
      any(!required_probabilities %in% rownames(relative_tro))) {
    stop(
      "A posterior TRO summary is missing its 2.5%, median, or 97.5% row.",
      call. = FALSE
    )
  }
  tibble(
    Source = source,
    Model = model,
    Year = fit$data$first_yr + seq_len(ncol(relative_tro)) - 1L,
    lower = as.numeric(relative_tro["0.025", ]),
    median = as.numeric(relative_tro["0.500", ]),
    upper = as.numeric(relative_tro["0.975", ])
  )
}

sensitivity_grid_tro_comparison <- bind_rows(
  imap_dfr(
    sensitivity_grid_comparison_fits,
    ~ posterior_tro_comparison_frame(
      .x,
      sensitivity_grid_comparison_labels[[.y]],
      "Sensitivity runs"
    )
  ),
  imap_dfr(
    grid_comparison_fits,
    ~ posterior_tro_comparison_frame(
      .x,
      paste0(
        "Grid cell ", .y, " (h=", format_decimal(grid_values$h[[.y]], 2),
        ", psi=", format_decimal(grid_values$psi[[.y]], 2), ")"
      ),
      "MCMC grid cells"
    )
  )
)

sensitivity_grid_terminal_rows <- sensitivity_grid_tro_comparison |>
  group_by(Source, Model) |>
  slice_max(Year, n = 1L, with_ties = FALSE) |>
  ungroup()
sensitivity_grid_terminal_spread <- sensitivity_grid_terminal_rows |>
  group_by(Source) |>
  summarise(
    Year = unique(Year),
    minimum = min(median),
    maximum = max(median),
    spread = maximum - minimum,
    minimum_model = Model[which.min(median)],
    maximum_model = Model[which.max(median)],
    .groups = "drop"
  )
sensitivity_terminal_spread <- sensitivity_grid_terminal_spread |>
  filter(Source == "Sensitivity runs")
grid_terminal_spread <- sensitivity_grid_terminal_spread |>
  filter(Source == "MCMC grid cells")
sensitivity_grid_terminal_spread_ratio <-
  sensitivity_terminal_spread$spread / grid_terminal_spread$spread

In 2026, the sensitivity posterior medians span 0.239 to 0.339 relative TRO, versus 0.244 to 0.318 across the nine MCMC grid cells. The sensitivity span is therefore 37.0% wider. CPUE q split from 2008 is below the lowest grid median, and Direct four-anchor M is above the highest grid median. The \(h\)/\(\psi\) grid is consequently substantial, but it does not envelop the full structural sensitivity range.

Show code

sensitivity_grid_tro_comparison |>
  ggplot(aes(
    x = Year,
    y = median,
    colour = Source,
    group = interaction(Source, Model)
  )) +
  geom_line(linewidth = 0.65, alpha = 0.72) +
  scale_colour_manual(values = c(
    "MCMC grid cells" = fit_expected_color,
    "Sensitivity runs" = "#D55E00"
  )) +
  labs(x = "Year", y = "Relative TRO", colour = NULL) +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero() +
  theme(legend.position = "bottom")

Figure 9: Posterior median relative total reproductive output for the accepted sensitivity runs and all nine accepted h/psi MCMC grid cells. The two colors distinguish sensitivity runs from MCMC-grid cells; every line is one fitted model.

Show code

sensitivity_grid_terminal_rows |>
  transmute(
    Source,
    Model,
    Year = as.character(as.integer(Year)),
    `2.5%` = format_decimal(lower, digits = 3),
    Median = format_decimal(median, digits = 3),
    `97.5%` = format_decimal(upper, digits = 3)
  ) |>
  kable()

Table 6: Terminal-year relative TRO posterior summaries for the accepted sensitivity fits and all nine MCMC grid cells.

Source	Model	Year	2.5%	Median	97.5%
MCMC grid cells	Grid cell 1 (h=0.60, psi=1.50)	2026	0.192	0.244	0.311
MCMC grid cells	Grid cell 2 (h=0.70, psi=1.50)	2026	0.222	0.280	0.351
MCMC grid cells	Grid cell 3 (h=0.80, psi=1.50)	2026	0.256	0.317	0.388
MCMC grid cells	Grid cell 4 (h=0.60, psi=1.75)	2026	0.194	0.247	0.314
MCMC grid cells	Grid cell 5 (h=0.70, psi=1.75)	2026	0.222	0.282	0.353
MCMC grid cells	Grid cell 6 (h=0.80, psi=1.75)	2026	0.257	0.317	0.390
MCMC grid cells	Grid cell 7 (h=0.60, psi=2.00)	2026	0.194	0.248	0.313
MCMC grid cells	Grid cell 8 (h=0.70, psi=2.00)	2026	0.224	0.280	0.351
MCMC grid cells	Grid cell 9 (h=0.80, psi=2.00)	2026	0.258	0.318	0.392
Sensitivity runs	Annual LL1 selectivity 2023-2025	2026	0.225	0.283	0.356
Sensitivity runs	Base	2026	0.222	0.282	0.353
Sensitivity runs	CPUE index with LL1 selectivity	2026	0.233	0.294	0.363
Sensitivity runs	CPUE power 0.75	2026	0.233	0.292	0.366
Sensitivity runs	CPUE q split from 2008	2026	0.182	0.239	0.305
Sensitivity runs	Constant CPUE CV 0.2	2026	0.224	0.283	0.352
Sensitivity runs	Direct four-anchor M	2026	0.266	0.339	0.432
Sensitivity runs	Drop CPUE 2021-2025	2026	0.214	0.270	0.338
Sensitivity runs	Estimated M-at-age slope	2026	0.228	0.287	0.362
Sensitivity runs	GTI with base aerial tau 0.59	2026	0.229	0.287	0.357
Sensitivity runs	More flexible Indonesian selectivity	2026	0.224	0.279	0.345
Sensitivity runs	No HSP	2026	0.226	0.287	0.363
Sensitivity runs	No POP or HSP	2026	0.201	0.266	0.349
Sensitivity runs	No UAM	2026	0.233	0.290	0.362

MLE grid

The previous OMMP meeting report described a deterministic MLE grid using natural-mortality values informed by MCMC posterior distributions. For the ESC31 workflow, the completed stacked length-based MCMC grid supplies both natural-mortality coordinates. Fixed M0 uses the 1%, 50%, and 95% posterior quantiles. M10 retains the existing three evenly spaced values spanning its 1% to 95% posterior range. The historical steepness levels are taken verbatim from sbt/data-raw/base22/base22.dat: 0.55, 0.63, 0.72, and 0.80.

This dependent grid deliberately changes mortality parameterization. Every cell uses M_switch = 1 and the package get_M() function, fixes M0 and M10, and estimates M4 and M30. Combining the four 2023 ADMB steepness values, three psi values, three M0 values, and three posterior-informed M10 values gives exactly 108 distinct fitted models. No inactive coordinate is duplicated or reused.

This is the executable implementation of the OMMP16 Og updated-MLE-grid test. Its natural-mortality values and results are generated only after the nine-cell MCMC grid passes the production gates; no earlier grid artifact is accepted as a final input. Missing or stale derived grid caches are never rebuilt during a normal render. ESC31_RUN_GRID_POSTPROCESSING=true enables accepted-grid postprocessing, and ESC31_RUN_MLE_GRID=true separately enables the 108-cell MLE run and resample. Convergence code zero alone is not accepted: each cell must also pass the maximum-gradient, estimability, Hessian, and biological-state gates. A cell that misses one of those gates is continued from its fitted state using a Hessian-diagonal-scaled bounded-gradient pass, a safeguarded bounded Newton correction when needed, and a final bounded-gradient certification pass. At most two staged recovery attempts are allowed. Already accepted cells remain cached, and the acceptance thresholds are unchanged. The bounded-Newton Armijo comparison allows an absolute 1e-8 objective tolerance solely for floating-point noise; that tolerance is about 1e-12 relative to the fitted grid objectives and does not replace the requirement for a decreasing gradient or the final strict certification gates.

Within each \(h\)/\(\psi\) group, resampling weights are calculated from differences in the full penalized sbt objective, not a likelihood component alone. Those conditional weights are then combined with equal prior weights over \(h\) and the specified 0.25/0.50/0.25 prior weights over \(\psi\). The resulting labels below therefore use full-objective-weighted rather than likelihood-only terminology.

Show code

# MLE-grid parameter construction, summarisation, state reconstruction, and
# validation are package APIs so the assessment page contains only the
# ESC31-specific grid design and acceptance policy.

Show code

dir.create(mle_grid_root, recursive = TRUE, showWarnings = FALSE)

m10_grid_range <- as.numeric(quantile(
  m10_posterior,
  probs = m10_grid_range_probs,
  na.rm = TRUE
))
m10_grid_values <- round(seq(m10_grid_range[1], m10_grid_range[2], length.out = m10_grid_n), 3)
m0_grid_values <- as.numeric(quantile(
  m0_posterior,
  probs = m0_grid_probs,
  na.rm = TRUE
))
if (length(m0_grid_values) != m0_grid_n ||
    any(!is.finite(m0_grid_values)) ||
    anyDuplicated(m0_grid_values)) {
  stop("The stacked posterior did not produce three distinct finite M0 quantiles.",
       call. = FALSE)
}
natural_mortality_grid_values <- bind_rows(
  tibble(
    parameter = "M0", position = seq_along(m0_grid_values),
    n_values = length(m0_grid_values), value = m0_grid_values,
    source = "1%, 50%, and 95% combined MCMC posterior quantiles"
  ),
  tibble(
    parameter = "M10", position = seq_along(m10_grid_values),
    n_values = length(m10_grid_values), value = m10_grid_values,
    source = "1%-95% combined MCMC posterior range"
  )
)

m10_mle_data <- base_fit$data
m10_mle_data$M_switch <-
  esc31_direct_mle_grid_specification$M_switch
m10_mle_parameters <- get_parameters(m10_mle_data)
shared_mle_parameter_names <- intersect(
  names(m10_mle_parameters),
  names(base_fit$parameters)
)
for (parameter_name in shared_mle_parameter_names) {
  if (!identical(
        length(m10_mle_parameters[[parameter_name]]),
        length(base_fit$parameters[[parameter_name]])
      ) ||
      !identical(
        dim(m10_mle_parameters[[parameter_name]]),
        dim(base_fit$parameters[[parameter_name]])
      )) {
    stop(
      "The direct-M grid and base parameter layouts differ for `",
      parameter_name, "`.",
      call. = FALSE
    )
  }
  m10_mle_parameters[[parameter_name]] <-
    base_fit$parameters[[parameter_name]]
}
base_mortality_report <- sbt_fit_report(base_fit)
m10_mle_parameters$par_log_m0 <-
  log(as.numeric(base_mortality_report$par_m0))
m10_mle_parameters$par_log_m4 <-
  log(as.numeric(base_mortality_report$par_m4))
m10_mle_parameters$par_log_m10 <-
  log(as.numeric(base_mortality_report$par_m10))
m10_mle_parameters$par_log_m30 <-
  log(as.numeric(base_mortality_report$par_m30))
m10_mle_data$priors <- get_priors(m10_mle_parameters)

m10_mle_grid <- expand.grid(
  h = mle_grid_h_values,
  psi = mle_grid_psi_values,
  m0 = m0_grid_values,
  m10 = m10_grid_values
) |>
  arrange(psi, h, m0, m10) |>
  mutate(Cell = row_number(), .before = 1)
if (nrow(m10_mle_grid) !=
      esc31_direct_mle_grid_specification$cells ||
    nrow(distinct(m10_mle_grid, h, psi, m0, m10)) !=
      esc31_direct_mle_grid_specification$cells) {
  stop("The direct-M dependent grid must contain 108 distinct cells.",
       call. = FALSE)
}

m10_mle_map <- base_fit$map
m10_mle_map <- m10_mle_map[
  intersect(names(m10_mle_map), names(m10_mle_parameters))
]
m10_mle_map$par_log_h <- factor(NA)
m10_mle_map$par_log_psi <- factor(NA)
m10_mle_map$par_log_m0 <- factor(NA)
m10_mle_map$par_log_m10 <- factor(NA)
direct_mle_fixed_parameters <-
  esc31_direct_mle_grid_specification$fixed_parameters
direct_mle_estimated_parameters <-
  esc31_direct_mle_grid_specification$estimated_mortality_parameters
direct_mle_parameter_is_fixed <- function(name) {
  name %in% names(m10_mle_map) &&
    length(m10_mle_map[[name]]) > 0L &&
    all(is.na(m10_mle_map[[name]]))
}
if (!all(c(
      direct_mle_fixed_parameters,
      direct_mle_estimated_parameters
    ) %in% names(m10_mle_parameters)) ||
    !all(vapply(
      direct_mle_fixed_parameters,
      direct_mle_parameter_is_fixed,
      logical(1)
    )) ||
    any(vapply(
      direct_mle_estimated_parameters,
      direct_mle_parameter_is_fixed,
      logical(1)
    ))) {
  stop(
    "The direct-M grid must fix h, psi, M0, and M10 while estimating M4 and M30.",
    call. = FALSE
  )
}
m10_mle_cells <- as.integer(m10_mle_grid$Cell)

m10_mle_signature <- esc31_object_md5(list(
  grid_specification = esc31_grid_specification,
  direct_mle_grid_specification =
    esc31_direct_mle_grid_specification,
  prerequisite_acceptance = grid_prerequisite_validation[
    setdiff(
      names(grid_prerequisite_validation),
      c(
        "manifest_signature",
        "validator_identity_signature",
        "base_posterior_source_checksum"
      )
    )
  ],
  grid_files = grid_signature(grid_diagnostics$file),
  combined_source_signature = combined_source_signature,
  base_fit_md5 = base_fit_md5,
  base_run_signature = base_fit$provenance$metadata$run_identity$signature,
  mortality_model = list(
    M_switch = m10_mle_data$M_switch,
    function_name =
      esc31_direct_mle_grid_specification$mortality_function,
    fixed = esc31_direct_mle_grid_specification$fixed_parameters,
    estimated =
      esc31_direct_mle_grid_specification$estimated_mortality_parameters
  ),
  m0_grid_values = m0_grid_values,
  m10_grid_values = m10_grid_values,
  h_grid_values = mle_grid_h_values,
  psi_grid_values = mle_grid_psi_values,
  # Retain this cache-schema field name, but store all 108 distinct cells.
  # No representative-cell expansion or fit reuse is performed.
  representative_cells = m10_mle_cells,
  check_estimability = check_m10_mle_estimability,
  max_gradient_limit = m10_mle_gradient_limit,
  biological_state_contract = biological_state_contract(),
  harvest_wall_contract = grid_harvest_wall_identity$executable,
  harvest_wall_decision = grid_harvest_wall_identity$decision,
  implementation =
    "direct_get_M_posterior_m0_m10_estimated_m4_m30_mle_grid_108_v14"
))

m10_mle_cell_summary_passes <- function(summary) {
  contract <- biological_state_contract()
  is.data.frame(summary) &&
    nrow(summary) == 1L &&
    is.finite(summary$objective[[1L]]) &&
    identical(as.integer(summary$convergence[[1L]]), 0L) &&
    is.finite(summary$max_gradient[[1L]]) &&
    summary$max_gradient[[1L]] <= m10_mle_gradient_limit &&
    isTRUE(summary$estimable[[1L]]) &&
    isTRUE(summary$state_passes[[1L]]) &&
    identical(as.integer(summary$state_invalid_cells[[1L]]), 0L) &&
    is.finite(summary$state_max_raw_harvest[[1L]]) &&
    summary$state_max_raw_harvest[[1L]] <=
      contract$hrate_limit + contract$hrate_tolerance &&
    is.finite(summary$state_min_number[[1L]]) &&
    summary$state_min_number[[1L]] > contract$number_tolerance &&
    is.finite(summary$state_min_spawning_biomass[[1L]]) &&
    summary$state_min_spawning_biomass[[1L]] > 0 &&
    is.finite(summary$state_min_recruitment[[1L]]) &&
    summary$state_min_recruitment[[1L]] > 0 &&
    is.finite(summary$state_max_harvest_penalty[[1L]]) &&
    abs(summary$state_max_harvest_penalty[[1L]]) <=
      contract$penalty_tolerance &&
    is.finite(summary$state_lp_penalty[[1L]]) &&
    abs(summary$state_lp_penalty[[1L]]) <=
      contract$penalty_tolerance &&
    is.finite(summary$state_max_catch_relative_error[[1L]]) &&
    summary$state_max_catch_relative_error[[1L]] <=
      contract$catch_relative_tolerance &&
    is.character(summary$state_checksum) &&
    grepl("^[0-9a-f]{32}$", summary$state_checksum[[1L]])
}

m10_mle_optimizer_recovery <- list(
  implementation =
    "scaled_bounded_gradient_newton_certification_v3",
  attempts = 2L,
  predecessor_methods = c(
    "hessian_diagonal_scaled_bounded_gradient_nlminb_v1",
    "scaled_bounded_gradient_newton_certification_v2"
  ),
  scaled = list(
    max_passes = 1L,
    control = list(
      eval.max = 15000L,
      iter.max = 5000L,
      rel.tol = 1e-10,
      x.tol = 1e-8
    ),
    hessian_diagonal_floor = 1e-8,
    scale_min = 1e-4,
    scale_max = 10
  ),
  newton = list(
    minimum_reciprocal_condition = 1e-12,
    initial_alpha = 1,
    fraction_to_boundary = 0.995,
    backtrack_factor = 0.5,
    maximum_backtracks = 12L,
    armijo_c1 = 1e-4,
    objective_numerical_tolerance = 1e-8,
    require_gradient_improvement = TRUE
  ),
  certification = list(
    max_passes = 1L,
    control = list(eval.max = 2000L, iter.max = 1000L)
  )
)

m10_mle_results <- if (file.exists(m10_mle_grid_file) && !rebuild_m10_mle_grid) {
  read_rds_cache(m10_mle_grid_file)
} else {
  NULL
}

m10_mle_cache_signature_valid <-
  is.list(m10_mle_results) &&
    is.character(m10_mle_results$signature) &&
    length(m10_mle_results$signature) == 1L &&
    !is.na(m10_mle_results$signature) &&
    grepl("^[0-9a-f]{32}$", m10_mle_results$signature)
m10_mle_cache_current <- m10_mle_cache_signature_valid &&
    all(c(
      "signature", "m0_quantile_probs", "m10_range_probs",
      "posterior_ranges",
      "natural_mortality_grid_values", "grid", "summary", "biomass",
      "state_records", "optimizer_recovery"
    ) %in% names(m10_mle_results)) &&
    identical(
      m10_mle_results$optimizer_recovery,
      m10_mle_optimizer_recovery
    ) &&
    isTRUE(all.equal(
      as.numeric(m10_mle_results$m0_quantile_probs),
      as.numeric(m0_grid_probs),
      tolerance = 0,
      check.attributes = FALSE
    )) &&
    isTRUE(all.equal(
      as.numeric(m10_mle_results$m10_range_probs),
      as.numeric(m10_grid_range_probs),
      tolerance = 0,
      check.attributes = FALSE
    )) &&
    isTRUE(all.equal(
      as.data.frame(m10_mle_results$natural_mortality_grid_values),
      as.data.frame(natural_mortality_grid_values),
      tolerance = 0,
      check.attributes = FALSE
    )) &&
    isTRUE(all.equal(
      as.data.frame(m10_mle_results$posterior_ranges),
      as.data.frame(bind_rows(
        tibble(
          parameter = "M0",
          probability = m0_grid_probs,
          value = m0_grid_values
        ),
        tibble(
          parameter = "M10",
          probability = m10_grid_range_probs,
          value = m10_grid_range
        )
      )),
      tolerance = 0,
      check.attributes = FALSE
    )) &&
    isTRUE(all.equal(
      as.data.frame(m10_mle_results$grid),
      as.data.frame(m10_mle_grid),
      tolerance = 0,
      check.attributes = FALSE
    )) &&
    is.data.frame(m10_mle_results$summary) &&
    nrow(m10_mle_results$summary) == length(m10_mle_cells) &&
    all(vapply(
      seq_len(nrow(m10_mle_results$summary)),
      function(i) {
        m10_mle_cell_summary_passes(
          m10_mle_results$summary[i, , drop = FALSE]
        )
      },
      logical(1)
    )) &&
    is.data.frame(m10_mle_results$biomass)
if (isTRUE(m10_mle_cache_current)) {
  m10_mle_cache_current <- sbt::validate_mle_grid_state_records(
    records = m10_mle_results$state_records,
    saved_summary = m10_mle_results$summary,
    saved_biomass = m10_mle_results$biomass,
    grid = m10_mle_grid,
    data = m10_mle_data,
    parameters = m10_mle_parameters,
    map = m10_mle_map,
    cores = grid_state_diagnostic_cores
  )
}
m10_mle_signature_audit <- if (
  isTRUE(m10_mle_cache_current) &&
    !identical(m10_mle_results$signature, m10_mle_signature)
) {
  list(
    source_signature = m10_mle_results$signature,
    current_scientific_signature = m10_mle_signature,
    compatibility = paste(
      "complete 108-cell input, optimum, gradient, estimability, Hessian,",
      "biological-state, biomass, and payload revalidation under the",
      "unchanged direct-M scientific contract"
    )
  )
} else {
  NULL
}

if (!isTRUE(m10_mle_cache_current)) {
  if (!run_m10_mle_grid) {
    stop(
      "The compatible 108-cell MLE grid is unavailable. Set ",
      "ESC31_RUN_MLE_GRID=true only after the approved MCMC grid passes every production gate.",
      call. = FALSE
    )
  }
  base_start <- base_fit$fit$opt$par
  m10_mle_cell_dir <- file.path(
    mle_grid_root,
    "direct_get_M_posterior_m0_m10_mle_grid_108_cells_v2"
  )
  dir.create(m10_mle_cell_dir, recursive = TRUE, showWarnings = FALSE)
  m10_mle_cell_file <- function(cell) {
    file.path(
      m10_mle_cell_dir,
      sprintf("cell%03d.rds", as.integer(cell))
    )
  }
  m10_mle_cell_record_usable <- function(record, cell) {
    is.list(record) &&
      identical(record$signature, m10_mle_signature) &&
      identical(as.integer(record$cell), as.integer(cell)) &&
      is.list(record$result) &&
      is.data.frame(record$result$summary) &&
      nrow(record$result$summary) == 1L &&
      identical(
        as.integer(record$result$summary$Cell),
        as.integer(cell)
      ) &&
      is.data.frame(record$result$biomass) &&
      is.list(record$result$state_records) &&
      length(record$result$state_records) == 1L
  }
  m10_mle_cell_record_current <- function(record, cell) {
    m10_mle_cell_record_usable(record, cell) &&
      m10_mle_cell_summary_passes(record$result$summary)
  }
  m10_mle_cell_record_start <- function(record, fallback) {
    candidate <- tryCatch(
      record$result$state_records[[1L]]$fitted_par,
      error = function(error) NULL
    )
    if (is.numeric(candidate) && length(candidate) &&
        !is.null(names(candidate)) && all(is.finite(candidate))) {
      candidate
    } else {
      fallback
    }
  }
  m10_mle_align_start <- function(cell_start, saved_start) {
    if (is.null(saved_start)) return(cell_start)
    cell_names <- sbt::expand_parameter_names(names(cell_start))
    saved_names <- sbt::expand_parameter_names(names(saved_start))
    common <- intersect(cell_names, saved_names)
    if (!length(common)) {
      stop("The saved cell state has no active parameters in common.",
           call. = FALSE)
    }
    cell_start[match(common, cell_names)] <-
      saved_start[match(common, saved_names)]
    cell_start
  }
  m10_mle_staged_refinement <- function(
      grid_row, saved_start, prior_history = NULL,
      scaled_already_completed = FALSE) {
    cell_parameters <- sbt::make_mle_grid_parameters(
      grid_row,
      m10_mle_parameters,
      allow_new_parameters = FALSE
    )[[1L]]
    obj <- RTMB::MakeADFun(
      func = sbt::cmb(sbt::sbt_model, m10_mle_data),
      parameters = cell_parameters,
      map = m10_mle_map,
      random = character()
    )
    cell_bounds <- sbt::get_bounds(
      obj = obj,
      parameters = cell_parameters
    )
    start <- m10_mle_align_start(obj$par, saved_start)
    point_diagnostics <- function(par) {
      finite_parameters <- is.numeric(par) &&
        length(par) == length(obj$par) &&
        identical(names(par), names(obj$par)) &&
        all(is.finite(par))
      within_bounds <- finite_parameters &&
        all(par >= cell_bounds$lower) &&
        all(par <= cell_bounds$upper)
      objective <- if (finite_parameters && within_bounds) {
        tryCatch(as.numeric(obj$fn(par)), error = function(error) Inf)
      } else {
        Inf
      }
      finite_objective <- length(objective) == 1L &&
        is.finite(objective)
      gradient <- if (finite_objective) {
        tryCatch(as.numeric(obj$gr(par)), error = function(error) numeric())
      } else {
        numeric()
      }
      finite_gradient <- length(gradient) == length(par) &&
        all(is.finite(gradient))
      list(
        valid = finite_parameters && within_bounds &&
          finite_objective && finite_gradient,
        objective = if (finite_objective) objective else Inf,
        gradient = gradient,
        max_gradient = if (finite_gradient) {
          max(abs(gradient))
        } else {
          Inf
        }
      )
    }
    start_point <- point_diagnostics(start)
    if (!isTRUE(start_point$valid)) {
      stop("The saved grid-cell state is not a valid refinement start.",
           call. = FALSE)
    }

    best <- list(
      par = start,
      objective = start_point$objective,
      gradient = start_point$gradient,
      max_gradient = start_point$max_gradient,
      opt = list(
        par = start,
        objective = start_point$objective,
        convergence = 1L,
        iterations = 0L,
        evaluations = c("function" = 1L, "gradient" = 1L),
        message = "saved state awaiting staged recovery"
      )
    )
    history <- list()
    nlminb_calls <- 0L
    append_history <- function(
        stage, point, result = NULL, accepted = FALSE,
        scale = NULL, nonpositive_hessian_diagonal = NA_integer_,
        hessian_reciprocal_condition = NA_real_, alpha = NA_real_,
        armijo_pass = NA, gradient_improvement = NA) {
      history[[length(history) + 1L]] <<- data.frame(
        pass = length(history) + 1L,
        stage = as.character(stage),
        objective = as.numeric(point$objective),
        max_gradient = as.numeric(point$max_gradient),
        convergence = as.integer(
          if (is.null(result)) NA_integer_ else result$convergence
        ),
        accepted = isTRUE(accepted),
        scale_min = if (is.null(scale)) NA_real_ else min(scale),
        scale_max = if (is.null(scale)) NA_real_ else max(scale),
        nonpositive_hessian_diagonal =
          as.integer(nonpositive_hessian_diagonal),
        hessian_reciprocal_condition =
          as.numeric(hessian_reciprocal_condition),
        alpha = as.numeric(alpha),
        armijo_pass = as.logical(armijo_pass),
        gradient_improvement = as.logical(gradient_improvement),
        stringsAsFactors = FALSE
      )
      invisible(NULL)
    }
    retain_candidate <- function(result, point) {
      accepted <- isTRUE(point$valid) &&
        point$objective <= best$objective
      if (accepted) {
        result$objective <- point$objective
        best <<- list(
          par = result$par,
          objective = point$objective,
          gradient = point$gradient,
          max_gradient = point$max_gradient,
          opt = result
        )
      }
      accepted
    }
    run_bounded_nlminb <- function(stage, control, scale = NULL) {
      result <- if (is.null(scale)) {
        nlminb(
          start = best$par,
          objective = obj$fn,
          gradient = obj$gr,
          lower = cell_bounds$lower,
          upper = cell_bounds$upper,
          control = control
        )
      } else {
        nlminb(
          start = best$par,
          objective = obj$fn,
          gradient = obj$gr,
          scale = scale,
          lower = cell_bounds$lower,
          upper = cell_bounds$upper,
          control = control
        )
      }
      nlminb_calls <<- nlminb_calls + 1L
      point <- point_diagnostics(result$par)
      accepted <- retain_candidate(result, point)
      append_history(
        stage = stage,
        point = point,
        result = result,
        accepted = accepted,
        scale = scale
      )
      list(result = result, point = point, accepted = accepted)
    }
    best_certified <- function() {
      identical(as.integer(best$opt$convergence), 0L) &&
        is.finite(best$max_gradient) &&
        best$max_gradient <= m10_mle_gradient_limit
    }
    run_certification <- function(stage) {
      for (pass in seq_len(
        m10_mle_optimizer_recovery$certification$max_passes
      )) {
        run_bounded_nlminb(
          stage = paste(stage, pass),
          control = m10_mle_optimizer_recovery$certification$control
        )
        if (best_certified()) break
      }
      invisible(best_certified())
    }

    if (best$max_gradient <= m10_mle_gradient_limit) {
      run_certification("saved-state bounded certification")
    }

    if (!best_certified() && !isTRUE(scaled_already_completed)) {
      for (pass in seq_len(
        m10_mle_optimizer_recovery$scaled$max_passes
      )) {
        hessian_point <- best$par
        scaling_hessian <- obj$he(hessian_point)
        if (!is.matrix(scaling_hessian) ||
            !identical(
              dim(scaling_hessian),
              c(length(best$par), length(best$par))
            ) ||
            any(!is.finite(scaling_hessian))) {
          stop("The refinement scaling Hessian is invalid.", call. = FALSE)
        }
        hessian_diagonal <- diag(scaling_hessian)
        parameter_scale <- 1 / sqrt(pmax(
          abs(hessian_diagonal),
          m10_mle_optimizer_recovery$scaled$hessian_diagonal_floor
        ))
        parameter_scale <- pmin(
          pmax(
            parameter_scale,
            m10_mle_optimizer_recovery$scaled$scale_min
          ),
          m10_mle_optimizer_recovery$scaled$scale_max
        )
        stage_result <- run_bounded_nlminb(
          stage = paste("scaled bounded gradient nlminb", pass),
          control = m10_mle_optimizer_recovery$scaled$control,
          scale = parameter_scale
        )
        history[[length(history)]]$
          nonpositive_hessian_diagonal <-
            as.integer(sum(hessian_diagonal <= 0))
        if (best_certified()) break
        if (isTRUE(stage_result$accepted) &&
            best$max_gradient <= m10_mle_gradient_limit) {
          run_certification("post-scaled bounded certification")
          if (best_certified()) break
        }
      }
    }

    if (!best_certified()) {
      newton_contract <- m10_mle_optimizer_recovery$newton
      pre_newton <- best
      newton_hessian <- obj$he(pre_newton$par)
      valid_newton_hessian <- is.matrix(newton_hessian) &&
        identical(
          dim(newton_hessian),
          c(length(pre_newton$par), length(pre_newton$par))
        ) &&
        all(is.finite(newton_hessian))
      if (valid_newton_hessian) {
        newton_hessian <- (newton_hessian + t(newton_hessian)) / 2
        hessian_cholesky <- tryCatch(
          chol(newton_hessian),
          error = function(error) error
        )
        hessian_rcond <- tryCatch(
          as.numeric(rcond(newton_hessian)),
          error = function(error) NA_real_
        )
      } else {
        hessian_cholesky <- NULL
        hessian_rcond <- NA_real_
      }
      hessian_acceptable <- valid_newton_hessian &&
        !inherits(hessian_cholesky, "error") &&
        length(hessian_rcond) == 1L &&
        is.finite(hessian_rcond) &&
        hessian_rcond >=
          newton_contract$minimum_reciprocal_condition
      if (!hessian_acceptable) {
        append_history(
          stage = "bounded Newton correction unavailable",
          point = pre_newton,
          hessian_reciprocal_condition = hessian_rcond
        )
      } else {
        newton_step <- tryCatch(
          as.numeric(backsolve(
            hessian_cholesky,
            forwardsolve(
              t(hessian_cholesky),
              matrix(pre_newton$gradient, ncol = 1L)
            )
          )),
          error = function(error) numeric()
        )
        direction <- -newton_step
        descent_slope <- if (
          length(direction) == length(pre_newton$par) &&
          all(is.finite(direction))
        ) {
          sum(pre_newton$gradient * direction)
        } else {
          NA_real_
        }
        valid_direction <- length(direction) ==
          length(pre_newton$par) &&
          all(is.finite(direction)) &&
          any(direction != 0) &&
          is.finite(descent_slope) &&
          descent_slope < 0
        newton_accepted <- FALSE
        if (!valid_direction) {
          append_history(
            stage = "bounded Newton direction unavailable",
            point = pre_newton,
            hessian_reciprocal_condition = hessian_rcond
          )
        } else {
          upper_limited <- direction > 0 &
            is.finite(cell_bounds$upper)
          lower_limited <- direction < 0 &
            is.finite(cell_bounds$lower)
          feasible_limits <- c(
            (
              cell_bounds$upper[upper_limited] -
                pre_newton$par[upper_limited]
            ) / direction[upper_limited],
            (
              cell_bounds$lower[lower_limited] -
                pre_newton$par[lower_limited]
            ) / direction[lower_limited]
          )
          feasible_limits <- feasible_limits[
            is.finite(feasible_limits) & feasible_limits >= 0
          ]
          maximum_feasible_alpha <- if (length(feasible_limits)) {
            min(feasible_limits)
          } else {
            Inf
          }
          initial_alpha <- newton_contract$initial_alpha
          if (is.finite(maximum_feasible_alpha)) {
            initial_alpha <- min(
              initial_alpha,
              newton_contract$fraction_to_boundary *
                maximum_feasible_alpha
            )
          }
          if (is.finite(initial_alpha) && initial_alpha > 0) {
            for (backtrack in 0:newton_contract$maximum_backtracks) {
              alpha <- initial_alpha *
                newton_contract$backtrack_factor^backtrack
              candidate_par <- pre_newton$par + alpha * direction
              candidate_point <- point_diagnostics(candidate_par)
              armijo_pass <- isTRUE(candidate_point$valid) &&
                candidate_point$objective <=
                  pre_newton$objective +
                  newton_contract$armijo_c1 *
                    alpha * descent_slope +
                  newton_contract$objective_numerical_tolerance
              gradient_improvement <-
                is.finite(candidate_point$max_gradient) &&
                candidate_point$max_gradient <
                  pre_newton$max_gradient
              gradient_requirement_pass <-
                !isTRUE(
                  newton_contract$require_gradient_improvement
                ) ||
                gradient_improvement
              accepted_newton <- armijo_pass &&
                gradient_requirement_pass
              append_history(
                stage = "bounded Newton line search",
                point = candidate_point,
                accepted = accepted_newton,
                hessian_reciprocal_condition = hessian_rcond,
                alpha = alpha,
                armijo_pass = armijo_pass,
                gradient_improvement = gradient_improvement
              )
              if (accepted_newton) {
                best <- list(
                  par = candidate_par,
                  objective = candidate_point$objective,
                  gradient = candidate_point$gradient,
                  max_gradient = candidate_point$max_gradient,
                  opt = list(
                    par = candidate_par,
                    objective = candidate_point$objective,
                    convergence = 1L,
                    iterations = 0L,
                    evaluations = c(
                      "function" = 1L,
                      "gradient" = 1L
                    ),
                    message =
                      "accepted Newton correction awaiting certification"
                  )
                )
                newton_accepted <- TRUE
                break
              }
            }
          }
        }
        if (isTRUE(newton_accepted)) {
          run_certification("post-Newton bounded certification")
        }
      }
    }

    start_objective <- start_point$objective
    prior_history <- if (is.data.frame(prior_history) &&
        nrow(prior_history)) {
      prior_history$recovery_source <- "predecessor"
      prior_history
    } else {
      NULL
    }
    current_history <- bind_rows(history)
    if (nrow(current_history)) {
      current_history$recovery_source <-
        m10_mle_optimizer_recovery$implementation
    }

    obj$par <- best$par
    obj$env$last.par.best <- best$par
    obj$fn(best$par)
    best$opt$par <- best$par
    best$opt$objective <- best$objective
    best$opt$convergence <- as.integer(best$opt$convergence)
    obj$opt <- best$opt
    obj$grid_start_fn <- start_objective
    obj$grid_b0_start_multiplier <- 1
    obj$grid_nlminb_passes <- nlminb_calls
    obj$grid_refinement_history <- bind_rows(
      prior_history,
      current_history
    )
    list(refined_cell = obj)
  }
  fit_m10_mle_cell <- function(cell) {
    path <- m10_mle_cell_file(cell)
    record <- if (file.exists(path) && !rebuild_m10_mle_grid) {
      tryCatch(read_rds_cache(path), error = function(error) NULL)
    } else {
      NULL
    }
    if (m10_mle_cell_record_current(record, cell)) {
      return(record)
    }

    grid_row <- m10_mle_grid[
      m10_mle_grid$Cell == cell,
      ,
      drop = FALSE
    ]
    if (!m10_mle_cell_record_usable(record, cell)) {
      started <- Sys.time()
      record <- tryCatch({
        fits <- run_grid(
          data = m10_mle_data,
          grid_parameters = sbt::make_mle_grid_parameters(
            grid_row,
            m10_mle_parameters,
            allow_new_parameters = FALSE
          ),
          bounds = NULL,
          map = m10_mle_map,
          control = base_fit$control,
          n_passes = 3L,
          start = base_start,
          b0_start_step = 1.25,
          b0_start_max = 10
        )
        result <- sbt::summarise_mle_grid(
          grid_fits = fits,
          grid = grid_row,
          data = m10_mle_data,
          check_estimability_cells = check_m10_mle_estimability
        )
        list(
          signature = m10_mle_signature,
          cell = as.integer(cell),
          primary_attempt = 1L,
          recovery_method = NULL,
          recovery_attempt = 0L,
          start = "base_mle",
          elapsed_seconds = as.numeric(
            difftime(Sys.time(), started, units = "secs")
          ),
          refinement_history = data.frame(),
          result = result,
          error = NULL
        )
      }, error = function(error) {
        list(
          signature = m10_mle_signature,
          cell = as.integer(cell),
          primary_attempt = 1L,
          recovery_method = NULL,
          recovery_attempt = 0L,
          start = "base_mle",
          elapsed_seconds = as.numeric(
            difftime(Sys.time(), started, units = "secs")
          ),
          refinement_history = data.frame(),
          result = NULL,
          error = conditionMessage(error)
        )
      })
      atomic_save_rds(record, path)
      if (m10_mle_cell_record_current(record, cell) ||
          !m10_mle_cell_record_usable(record, cell)) {
        return(record)
      }
    }

    previous_recovery_attempt <- if (
      identical(
        record$recovery_method,
        m10_mle_optimizer_recovery$implementation
      ) &&
      is.numeric(record$recovery_attempt) &&
      length(record$recovery_attempt) == 1L &&
      is.finite(record$recovery_attempt)
    ) {
      as.integer(record$recovery_attempt)
    } else {
      0L
    }
    if (previous_recovery_attempt >=
        m10_mle_optimizer_recovery$attempts) {
      return(record)
    }

    for (recovery_attempt in seq.int(
      previous_recovery_attempt + 1L,
      m10_mle_optimizer_recovery$attempts
    )) {
      started <- Sys.time()
      prior_record <- record
      continuation_start <- m10_mle_cell_record_start(
        prior_record,
        base_start
      )
      scaled_already_completed <- is.list(prior_record) &&
        is.character(prior_record$recovery_method) &&
        length(prior_record$recovery_method) == 1L &&
        prior_record$recovery_method %in% c(
          m10_mle_optimizer_recovery$implementation,
          m10_mle_optimizer_recovery$predecessor_methods
        ) &&
        is.data.frame(prior_record$refinement_history) &&
        nrow(prior_record$refinement_history) > 0L
      record <- tryCatch({
        fits <- m10_mle_staged_refinement(
          grid_row = grid_row,
          saved_start = continuation_start,
          prior_history = prior_record$refinement_history,
          scaled_already_completed = scaled_already_completed
        )
        result <- sbt::summarise_mle_grid(
          grid_fits = fits,
          grid = grid_row,
          data = m10_mle_data,
          check_estimability_cells = check_m10_mle_estimability
        )
        list(
          signature = m10_mle_signature,
          cell = as.integer(cell),
          primary_attempt = 1L,
          recovery_method =
            m10_mle_optimizer_recovery$implementation,
          recovery_attempt = as.integer(recovery_attempt),
          start = "previous_cell_state",
          elapsed_seconds = as.numeric(
            difftime(Sys.time(), started, units = "secs")
          ),
          refinement_history =
            fits[[1L]]$grid_refinement_history,
          result = result,
          error = NULL
        )
      }, error = function(error) {
        list(
          signature = m10_mle_signature,
          cell = as.integer(cell),
          primary_attempt = 1L,
          recovery_method =
            m10_mle_optimizer_recovery$implementation,
          recovery_attempt = as.integer(recovery_attempt),
          start = "previous_cell_state",
          elapsed_seconds = as.numeric(
            difftime(Sys.time(), started, units = "secs")
          ),
          refinement_history = data.frame(),
          result = prior_record$result,
          error = conditionMessage(error)
        )
      })
      atomic_save_rds(record, path)
      if (m10_mle_cell_record_current(record, cell)) break
    }
    record
  }

  m10_mle_cell_records <- if (
    .Platform$OS.type != "windows" && m10_mle_grid_cores > 1L
  ) {
    parallel::mclapply(
      m10_mle_cells,
      fit_m10_mle_cell,
      mc.cores = min(
        m10_mle_grid_cores,
        length(m10_mle_cells)
      ),
      mc.preschedule = FALSE,
      mc.set.seed = FALSE
    )
  } else {
    lapply(m10_mle_cells, fit_m10_mle_cell)
  }
  failed_m10_mle_cells <- m10_mle_cells[!vapply(
    seq_along(m10_mle_cell_records),
    function(i) {
      m10_mle_cell_record_current(
        m10_mle_cell_records[[i]],
        m10_mle_cells[[i]]
      )
    },
    logical(1)
  )]
  if (length(failed_m10_mle_cells)) {
    failure_record_index <- match(
      failed_m10_mle_cells,
      m10_mle_cells
    )
    failure_messages <- vapply(
      failure_record_index,
      function(i) {
        m10_mle_cell_records[[i]]$error %||%
          "invalid cell cache"
      },
      character(1)
    )
    stop(
      "Production MLE-grid cell failures: ",
      paste0(
        failed_m10_mle_cells, " (", failure_messages, ")",
        collapse = "; "
      ),
      call. = FALSE
    )
  }
  names(m10_mle_cell_records) <- as.character(m10_mle_cells)
  m10_mle_grid_summary <- list(
    summary = bind_rows(lapply(
      m10_mle_cell_records,
      function(record) record$result$summary
    )),
    biomass = bind_rows(lapply(
      m10_mle_cell_records,
      function(record) record$result$biomass
    )),
    state_records = lapply(
      m10_mle_cell_records,
      function(record) record$result$state_records[[1L]]
    )
  )
  if (!sbt::validate_mle_grid_state_records(
        records = m10_mle_grid_summary$state_records,
        saved_summary = m10_mle_grid_summary$summary,
        saved_biomass = m10_mle_grid_summary$biomass,
        grid = m10_mle_grid,
        data = m10_mle_data,
        parameters = m10_mle_parameters,
        map = m10_mle_map,
        cores = grid_state_diagnostic_cores
      )) {
    stop(
      "The direct-M production MLE-grid records failed independent ",
      "108-cell reconstruction.",
      call. = FALSE
    )
  }

  m10_mle_results <- list(
    signature = m10_mle_signature,
    m0_quantile_probs = m0_grid_probs,
    m10_range_probs = m10_grid_range_probs,
    posterior_ranges = bind_rows(
      tibble(
        parameter = "M0",
        probability = m0_grid_probs,
        value = m0_grid_values
      ),
      tibble(
        parameter = "M10",
        probability = m10_grid_range_probs,
        value = m10_grid_range
      )
    ),
    natural_mortality_grid_values = natural_mortality_grid_values,
    grid = m10_mle_grid,
    summary = m10_mle_grid_summary$summary,
    biomass = m10_mle_grid_summary$biomass,
    state_records = m10_mle_grid_summary$state_records,
    optimizer_recovery = m10_mle_optimizer_recovery
  )
  atomic_save_rds(m10_mle_results, m10_mle_grid_file)
}

m10_mle_summary <- m10_mle_results$summary |>
  mutate(delta_nll = objective - min(objective, na.rm = TRUE))
required_mle_state_fields <- c(
  "state_passes", "state_invalid_cells", "state_max_raw_harvest",
  "state_min_number", "state_min_spawning_biomass",
  "state_min_recruitment", "state_max_harvest_penalty",
  "state_lp_penalty", "state_max_catch_relative_error", "state_checksum"
)
if (!all(required_mle_state_fields %in% names(m10_mle_summary)) ||
    nrow(m10_mle_summary) != 108L ||
    !identical(sort(as.integer(m10_mle_summary$Cell)), seq_len(108L)) ||
    any(!is.finite(m10_mle_summary$objective)) ||
    any(!is.finite(m10_mle_summary$max_gradient)) ||
    any(m10_mle_summary$convergence != 0L) ||
    any(m10_mle_summary$max_gradient > m10_mle_gradient_limit) ||
    any(is.na(m10_mle_summary$estimable) | !m10_mle_summary$estimable) ||
    any(is.na(m10_mle_summary$state_passes) | !m10_mle_summary$state_passes) ||
    any(!is.finite(m10_mle_summary$state_invalid_cells)) ||
    any(m10_mle_summary$state_invalid_cells != 0L) ||
    any(!is.finite(m10_mle_summary$state_max_raw_harvest)) ||
    any(m10_mle_summary$state_max_raw_harvest >
      biological_state_contract()$hrate_limit +
      biological_state_contract()$hrate_tolerance) ||
    any(!is.finite(m10_mle_summary$state_min_number)) ||
    any(m10_mle_summary$state_min_number <=
      biological_state_contract()$number_tolerance) ||
    any(!is.finite(m10_mle_summary$state_min_spawning_biomass) |
      m10_mle_summary$state_min_spawning_biomass <= 0) ||
    any(!is.finite(m10_mle_summary$state_min_recruitment) |
      m10_mle_summary$state_min_recruitment <= 0) ||
    any(!is.finite(m10_mle_summary$state_max_harvest_penalty)) ||
    any(abs(m10_mle_summary$state_max_harvest_penalty) >
      biological_state_contract()$penalty_tolerance) ||
    any(!is.finite(m10_mle_summary$state_lp_penalty)) ||
    any(abs(m10_mle_summary$state_lp_penalty) >
      biological_state_contract()$penalty_tolerance) ||
    any(!is.finite(m10_mle_summary$state_max_catch_relative_error)) ||
    any(m10_mle_summary$state_max_catch_relative_error >
      biological_state_contract()$catch_relative_tolerance) ||
    anyNA(m10_mle_summary$state_checksum) ||
    any(!grepl("^[0-9a-f]{32}$", m10_mle_summary$state_checksum))) {
  stop(
    "The 108-cell MLE grid contains missing, failed, non-estimable, ",
    "high-gradient, or biologically invalid cells and cannot be resampled.",
    call. = FALSE
  )
}
m10_mle_biomass <- m10_mle_results$biomass
m10_mle_payload_checksum <- esc31_object_md5(list(
  summary = m10_mle_results$summary,
  biomass = m10_mle_results$biomass,
  state_records = m10_mle_results$state_records
))
mle_grid_projection_map <- m10_mle_map
mle_grid_projection_map$par_log_h <- NULL
mle_grid_projection_map$par_log_psi <- NULL
mle_grid_projection_map$par_log_m0 <- NULL
mle_grid_projection_map$par_log_m10 <- NULL
obj_mle_grid_projection <- MakeADFun(
  func = cmb(sbt_model, m10_mle_data),
  parameters = m10_mle_parameters,
  map = mle_grid_projection_map
)
mle_grid_target_par_names <-
  expand_parameter_names(names(obj_mle_grid_projection$par))
mle_grid_target_sample_names <- c(mle_grid_target_par_names, "lp__")
mle_grid_conversion_identity <- esc31_sbt_function_identity(c(
  "grid_to_tmbfit",
  "make_mle_grid_parameters",
  "expand_parameter_names"
))
mle_grid_tmbfit_signature <- esc31_object_md5(list(
  # Bind the resample to the fully revalidated source artifact. The current
  # scientific signature is audited separately, so package-version changes do
  # not rewrite accepted source files or invalidate downstream projections.
  m10_mle_signature = m10_mle_results$signature,
  m10_mle_payload_checksum = m10_mle_payload_checksum,
  n_samples = mle_grid_resample_n,
  seed = mle_grid_resample_seed,
  target_parameter_names = mle_grid_target_par_names,
  target_map = mle_grid_projection_map,
  optimizer_recovery = m10_mle_optimizer_recovery,
  harvest_wall_contract = grid_harvest_wall_identity$executable,
  harvest_wall_decision = grid_harvest_wall_identity$decision,
  implementation =
    "direct_get_M_mle_grid_full_state_to_tmbfit_v14_posterior_m0"
))
legacy_mle_grid_tmbfit_signature <- esc31_object_md5(list(
  m10_mle_signature = m10_mle_results$signature,
  m10_mle_payload_checksum = m10_mle_payload_checksum,
  n_samples = mle_grid_resample_n,
  seed = mle_grid_resample_seed,
  target_parameter_names = mle_grid_target_par_names,
  target_map = mle_grid_projection_map,
  conversion_identity = mle_grid_conversion_identity,
  harvest_wall_contract = grid_harvest_wall_identity$executable,
  harvest_wall_decision = grid_harvest_wall_identity$decision,
  implementation =
    "direct_get_M_mle_grid_full_state_to_tmbfit_v12_canonical_checksum_parallel_projection_map_harvest_wall_bound"
))
mle_grid_tmbfit_audit <- list(
  source_mle_grid_signature = m10_mle_results$signature,
  conversion_identity = mle_grid_conversion_identity
)
expected_mle_grid_sample <- sample_grid(
  grid = m10_mle_summary,
  n_samples = mle_grid_resample_n,
  seed = mle_grid_resample_seed
)
mle_grid_resample_payload <- function(sample, fit) {
  if (!is.list(sample) ||
      !identical(names(sample), names(expected_mle_grid_sample)) ||
      !inherits(fit, "tmbfit") || !is.list(fit) ||
      !is.array(fit$samples) || anyNA(fit$samples) ||
      any(!is.finite(fit$samples)) ||
      !identical(fit$par_names, mle_grid_target_par_names) ||
      !identical(fit$sample_names, mle_grid_target_sample_names) ||
      !identical(
        dimnames(fit$samples)[[3L]],
        mle_grid_target_sample_names
      ) ||
      !identical(as.integer(fit$grid_cells),
                 as.integer(sample$grid_cells)) ||
      !is.data.frame(fit$draw_metadata) ||
      nrow(fit$draw_metadata) != length(sample$grid_cells)) {
    return(NULL)
  }
  list(
    sample = sample,
    fit = list(
      samples = fit$samples,
      sampler_params = fit$sampler_params,
      model = fit$model,
      metric = fit$metric,
      par_names = fit$par_names,
      sample_names = fit$sample_names,
      grid_cells = fit$grid_cells,
      draw_metadata = fit$draw_metadata,
      max_treedepth = fit$max_treedepth,
      warmup = fit$warmup,
      iter = fit$iter,
      thin = fit$thin,
      algorithm = fit$algorithm,
      class = class(fit)
    )
  )
}
mle_grid_resample_checksum <- function(payload) {
  if (is.null(payload)) return(NA_character_)
  canonical_payload <- unserialize(serialize(
    payload,
    connection = NULL,
    version = 3
  ))
  esc31_object_md5(canonical_payload)
}
build_mle_grid_tmbfit <- function(sample = expected_mle_grid_sample) {
  sbt::grid_to_tmbfit(
    data = m10_mle_data,
    parameters = m10_mle_parameters,
    grid = m10_mle_grid,
    grid_parameters = sbt::make_mle_grid_parameters(
      m10_mle_grid,
      m10_mle_parameters,
      allow_new_parameters = FALSE
    ),
    grid_cells = sample,
    fitted_parameters = m10_mle_results$state_records,
    source_map = m10_mle_map,
    target_map = mle_grid_projection_map,
    cores = grid_state_diagnostic_cores
  )
}
mle_grid_tmbfit_cache <- if (file.exists(m10_mle_tmbfit_file)) {
  read_rds_cache(m10_mle_tmbfit_file)
} else {
  NULL
}
cached_mle_grid_resample_payload <- tryCatch(
  mle_grid_resample_payload(
    mle_grid_tmbfit_cache$sample,
    mle_grid_tmbfit_cache$fit
  ),
  error = function(error) NULL
)
mle_grid_tmbfit_cache_record_valid <- is.list(mle_grid_tmbfit_cache) &&
  all(c(
    "signature", "sample", "fit", "payload_checksum"
  ) %in% names(mle_grid_tmbfit_cache))
# The accepted cached draw sequence is part of the signed payload. Recreating
# `sample()` is not a stable identity check across R sampling implementations;
# validate the complete cached payload and its source-sensitive signature.
mle_grid_tmbfit_cache_content_current <-
  mle_grid_tmbfit_cache_record_valid &&
  !is.null(cached_mle_grid_resample_payload) &&
  identical(
    mle_grid_tmbfit_cache$payload_checksum,
    mle_grid_resample_checksum(cached_mle_grid_resample_payload)
  )
mle_grid_tmbfit_signature_current <-
  mle_grid_tmbfit_cache_record_valid &&
  is.character(mle_grid_tmbfit_cache$signature) &&
  length(mle_grid_tmbfit_cache$signature) == 1L &&
  !is.na(mle_grid_tmbfit_cache$signature) &&
  mle_grid_tmbfit_cache$signature %in% c(
    mle_grid_tmbfit_signature,
    legacy_mle_grid_tmbfit_signature
  )
mle_grid_tmbfit_cache_current <-
  mle_grid_tmbfit_cache_content_current &&
  mle_grid_tmbfit_signature_current
mle_grid_tmbfit_migration_reason <- NULL
if (mle_grid_tmbfit_cache_content_current &&
    !mle_grid_tmbfit_signature_current &&
    is.character(mle_grid_tmbfit_cache$signature) &&
    length(mle_grid_tmbfit_cache$signature) == 1L &&
    !is.na(mle_grid_tmbfit_cache$signature) &&
    grepl("^[0-9a-f]{32}$", mle_grid_tmbfit_cache$signature)) {
  rebuilt_mle_grid_tmbfit <- build_mle_grid_tmbfit()
  rebuilt_mle_grid_payload_checksum <- mle_grid_resample_checksum(
    mle_grid_resample_payload(
      expected_mle_grid_sample,
      rebuilt_mle_grid_tmbfit
    )
  )
  if (identical(
        mle_grid_tmbfit_cache$payload_checksum,
        rebuilt_mle_grid_payload_checksum
      )) {
    mle_grid_tmbfit_cache_current <- TRUE
    mle_grid_tmbfit_migration_reason <-
      "verified_source_sensitive_signature_migration"
  }
}
if (isTRUE(mle_grid_tmbfit_cache_current)) {
  mle_grid_sample <- mle_grid_tmbfit_cache$sample
  mle_grid_tmbfit <- mle_grid_tmbfit_cache$fit
} else {
  if (!run_m10_mle_grid) {
    stop(
      "The compatible 108-cell MLE-grid resample is unavailable. Set ",
      "ESC31_RUN_MLE_GRID=true only after the approved MCMC grid passes every ",
      "production gate.",
      call. = FALSE
    )
  }
  mle_grid_sample <- expected_mle_grid_sample
  mle_grid_tmbfit <- build_mle_grid_tmbfit(mle_grid_sample)
  mle_grid_resample_checksum_value <- mle_grid_resample_checksum(
    mle_grid_resample_payload(mle_grid_sample, mle_grid_tmbfit)
  )
  atomic_save_rds(
    list(
      signature = mle_grid_tmbfit_signature,
      sample = mle_grid_sample,
      fit = mle_grid_tmbfit,
      payload_checksum = mle_grid_resample_checksum_value,
      audit = mle_grid_tmbfit_audit
    ),
    m10_mle_tmbfit_file
  )
}
mle_grid_freq <- mle_grid_sample$grid_freq |>
  mutate(delta_nll = nll - min(nll, na.rm = TRUE))
m_posterior_df <- bind_rows(
  tibble(parameter = "M0", value = m0_posterior),
  tibble(parameter = "M10", value = m10_posterior)
)

All 108 direct-M cells are distinct and pass. Objectives range from 7,610.100 to 7,617.852. The largest maximum gradient is 0.009225 (cell 51), maximum raw harvest is 0.668, and minimum populated abundance is 1,168.756. Estimated M4 ranges from 0.137 to 0.152, and estimated M30 ranges from 0.427 to 0.465. The validated projection-ready object contains 2,000 draws and samples 104 of the 108 cells with nonzero frequency.

Show code

natural_mortality_grid_values |>
  mutate(
    Parameter = parameter,
    Position = paste0(position, " of ", n_values),
    value = format_decimal(value, digits = 3)
  ) |>
  select(Parameter, Position, Value = value) |>
  kable(caption = "Natural-mortality values used in the direct get_M grid. M0 uses the 1%, 50%, and 95% stacked-MCMC posterior quantiles; M10 retains the existing three evenly spaced values across the 1%-95% combined-MCMC posterior range.")

Table 7: Natural-mortality values used in the direct get_M grid. M0 uses the 1%, 50%, and 95% stacked-MCMC posterior quantiles; M10 retains the existing three evenly spaced values across the 1%-95% combined-MCMC posterior range.

Parameter	Position	Value
M0	1 of 3	0.339
M0	2 of 3	0.368
M0	3 of 3	0.389
M10	1 of 3	0.099
M10	2 of 3	0.107
M10	3 of 3	0.114

Show code

m10_mle_summary |>
  transmute(
    Cell = format_decimal(Cell, digits = 0),
    h,
    psi,
    `Fixed M0` = fixed_m0,
    `Estimated M4` = estimated_m4,
    `Fixed M10` = fixed_m10,
    `Estimated M30` = estimated_m30,
    Objective = objective,
    `MLE convergence code` = if_else(is.na(convergence), NA_character_, as.character(as.integer(convergence))),
    `Max gradient` = max_gradient,
    `Estimability check` = case_when(
      !is.na(estimable) & estimable ~ "Pass",
      !is.na(estimable) & !estimable ~ "Fail",
      TRUE ~ "Not available"
    ),
    `Biological state` = if_else(state_passes, "Pass", "Fail"),
    `Max raw harvest` = state_max_raw_harvest,
    `Min abundance` = state_min_number,
    `Max continuation penalty` = state_max_harvest_penalty,
    `Max catch relative error` = state_max_catch_relative_error
  ) |>
  mutate(
    h = format_decimal(h, digits = 2),
    psi = format_decimal(psi, digits = 2)
  ) |>
  mutate(across(where(is.numeric), ~ format_decimal(.x, digits = 3))) |>
  replace_missing() |>
  kable(caption = "Production direct-get_M grid: 108 distinct fitted models from four 2023 ADMB h values, three psi values, three fixed posterior-quantile M0 values, and three fixed posterior-informed M10 values. M4 and M30 are estimated in every cell.")

Table 8: Production direct-get_M grid: 108 distinct fitted models from four 2023 ADMB h values, three psi values, three fixed posterior-quantile M0 values, and three fixed posterior-informed M10 values. M4 and M30 are estimated in every cell.

Cell	h	psi	Fixed M0	Estimated M4	Fixed M10	Estimated M30	Objective	Max gradient	Estimability check	Biological state	Max raw harvest	Min abundance
1	0.55	1.50	0.339	0.151	0.099	0.460	7,617.193	0.000	Pass	Pass	0.584	1,794.900
2	0.55	1.50	0.339	0.148	0.107	0.450	7,615.409	0.001	Pass	Pass	0.584	1,528.017
3	0.55	1.50	0.339	0.146	0.114	0.441	7,614.109	0.000	Pass	Pass	0.584	1,319.920
4	0.55	1.50	0.368	0.148	0.099	0.460	7,616.352	0.000	Pass	Pass	0.587	1,810.496
5	0.55	1.50	0.368	0.145	0.107	0.450	7,614.537	0.000	Pass	Pass	0.587	1,541.741
6	0.55	1.50	0.368	0.143	0.114	0.441	7,613.213	0.000	Pass	Pass	0.587	1,332.185
7	0.55	1.50	0.389	0.146	0.099	0.461	7,616.065	0.000	Pass	Pass	0.589	1,821.462
8	0.55	1.50	0.389	0.143	0.107	0.450	7,614.226	0.000	Pass	Pass	0.589	1,551.455
9	0.55	1.50	0.389	0.141	0.114	0.441	7,612.881	0.000	Pass	Pass	0.589	1,340.907
10	0.63	1.50	0.339	0.151	0.099	0.455	7,615.740	0.000	Pass	Pass	0.586	1,673.646
11	0.63	1.50	0.339	0.149	0.107	0.445	7,614.239	0.000	Pass	Pass	0.586	1,424.117
12	0.63	1.50	0.339	0.146	0.114	0.436	7,613.175	0.000	Pass	Pass	0.586	1,230.105
13	0.63	1.50	0.368	0.148	0.099	0.456	7,614.852	0.000	Pass	Pass	0.589	1,688.017
14	0.63	1.50	0.368	0.146	0.107	0.445	7,613.325	0.000	Pass	Pass	0.589	1,436.799
15	0.63	1.50	0.368	0.144	0.114	0.436	7,612.238	0.000	Pass	Pass	0.589	1,241.450
16	0.63	1.50	0.389	0.147	0.099	0.456	7,614.532	0.000	Pass	Pass	0.591	1,698.175
17	0.63	1.50	0.389	0.144	0.107	0.446	7,612.984	0.000	Pass	Pass	0.591	1,445.814
18	0.63	1.50	0.389	0.142	0.114	0.436	7,611.880	0.000	Pass	Pass	0.591	1,249.555
19	0.72	1.50	0.339	0.152	0.099	0.451	7,616.286	0.000	Pass	Pass	0.614	1,601.723
20	0.72	1.50	0.339	0.149	0.107	0.440	7,614.973	0.000	Pass	Pass	0.615	1,364.327
21	0.72	1.50	0.339	0.147	0.114	0.431	7,614.059	0.000	Pass	Pass	0.615	1,179.872
22	0.72	1.50	0.368	0.149	0.099	0.451	7,615.364	0.000	Pass	Pass	0.611	1,615.181
23	0.72	1.50	0.368	0.146	0.107	0.441	7,614.029	0.000	Pass	Pass	0.611	1,376.261
24	0.72	1.50	0.368	0.144	0.114	0.431	7,613.096	0.000	Pass	Pass	0.612	1,190.591
25	0.72	1.50	0.389	0.147	0.099	0.451	7,615.022	0.000	Pass	Pass	0.608	1,624.743
26	0.72	1.50	0.389	0.144	0.107	0.441	7,613.668	0.000	Pass	Pass	0.609	1,384.788
27	0.72	1.50	0.389	0.142	0.114	0.432	7,612.720	0.000	Pass	Pass	0.609	1,198.278
28	0.80	1.50	0.339	0.152	0.099	0.447	7,617.852	0.000	Pass	Pass	0.668	1,583.207
29	0.80	1.50	0.339	0.149	0.107	0.437	7,616.606	0.000	Pass	Pass	0.668	1,350.067
30	0.80	1.50	0.339	0.147	0.114	0.427	7,615.740	0.000	Pass	Pass	0.668	1,168.756
31	0.80	1.50	0.368	0.149	0.099	0.448	7,616.916	0.001	Pass	Pass	0.664	1,596.291
32	0.80	1.50	0.368	0.146	0.107	0.437	7,615.650	0.000	Pass	Pass	0.664	1,361.734
33	0.80	1.50	0.368	0.144	0.114	0.428	7,614.767	0.000	Pass	Pass	0.665	1,179.278
34	0.80	1.50	0.389	0.147	0.099	0.448	7,616.563	0.000	Pass	Pass	0.661	1,605.620
35	0.80	1.50	0.389	0.144	0.107	0.438	7,615.281	0.000	Pass	Pass	0.661	1,370.091
36	0.80	1.50	0.389	0.142	0.114	0.428	7,614.384	0.000	Pass	Pass	0.662	1,186.840
37	0.55	1.75	0.339	0.149	0.099	0.462	7,615.699	0.000	Pass	Pass	0.585	1,892.485
38	0.55	1.75	0.339	0.146	0.107	0.452	7,613.773	0.000	Pass	Pass	0.585	1,618.504
39	0.55	1.75	0.339	0.144	0.114	0.443	7,612.347	0.000	Pass	Pass	0.585	1,403.837
40	0.55	1.75	0.368	0.146	0.099	0.462	7,614.849	0.004	Pass	Pass	0.588	1,907.939
41	0.55	1.75	0.368	0.143	0.107	0.452	7,612.892	0.000	Pass	Pass	0.588	1,632.199
42	0.55	1.75	0.368	0.141	0.114	0.443	7,611.440	0.000	Pass	Pass	0.588	1,416.160
43	0.55	1.75	0.389	0.144	0.099	0.463	7,614.554	0.001	Pass	Pass	0.590	1,918.821
44	0.55	1.75	0.389	0.141	0.107	0.453	7,612.573	0.001	Pass	Pass	0.590	1,641.899
45	0.55	1.75	0.389	0.139	0.114	0.444	7,611.101	0.000	Pass	Pass	0.590	1,424.925
46	0.63	1.75	0.339	0.149	0.099	0.458	7,614.243	0.000	Pass	Pass	0.587	1,769.223
47	0.63	1.75	0.339	0.147	0.107	0.448	7,612.602	0.000	Pass	Pass	0.587	1,512.374
48	0.63	1.75	0.339	0.144	0.114	0.438	7,611.411	0.000	Pass	Pass	0.586	1,311.591
49	0.63	1.75	0.368	0.147	0.099	0.458	7,613.346	0.000	Pass	Pass	0.589	1,783.497
50	0.63	1.75	0.368	0.144	0.107	0.448	7,611.678	0.000	Pass	Pass	0.589	1,525.065
51	0.63	1.75	0.368	0.142	0.114	0.439	7,610.465	0.009	Pass	Pass	0.589	1,323.037
52	0.63	1.75	0.389	0.145	0.099	0.458	7,613.019	0.003	Pass	Pass	0.591	1,793.596
53	0.63	1.75	0.389	0.142	0.107	0.448	7,611.331	0.000	Pass	Pass	0.591	1,534.095
54	0.63	1.75	0.389	0.140	0.114	0.439	7,610.100	0.003	Pass	Pass	0.591	1,331.210
55	0.72	1.75	0.339	0.150	0.099	0.453	7,614.773	0.000	Pass	Pass	0.605	1,696.165
56	0.72	1.75	0.339	0.147	0.107	0.443	7,613.322	0.001	Pass	Pass	0.605	1,451.340
57	0.72	1.75	0.339	0.145	0.114	0.434	7,612.283	0.002	Pass	Pass	0.605	1,260.018
58	0.72	1.75	0.368	0.147	0.099	0.454	7,613.842	0.000	Pass	Pass	0.601	1,709.557
59	0.72	1.75	0.368	0.144	0.107	0.444	7,612.368	0.000	Pass	Pass	0.601	1,463.315
60	0.72	1.75	0.368	0.142	0.114	0.434	7,611.310	0.000	Pass	Pass	0.602	1,270.863
61	0.72	1.75	0.389	0.145	0.099	0.454	7,613.493	0.000	Pass	Pass	0.599	1,719.080
62	0.72	1.75	0.389	0.142	0.107	0.444	7,612.001	0.000	Pass	Pass	0.599	1,471.872
63	0.72	1.75	0.389	0.140	0.114	0.435	7,610.927	0.000	Pass	Pass	0.599	1,278.636
64	0.80	1.75	0.339	0.150	0.099	0.450	7,616.323	0.000	Pass	Pass	0.657	1,677.902
65	0.80	1.75	0.339	0.147	0.107	0.440	7,614.942	0.000	Pass	Pass	0.657	1,437.317
66	0.80	1.75	0.339	0.145	0.114	0.430	7,613.954	0.000	Pass	Pass	0.657	1,249.110
67	0.80	1.75	0.368	0.147	0.099	0.450	7,615.377	0.000	Pass	Pass	0.653	1,690.933
68	0.80	1.75	0.368	0.144	0.107	0.440	7,613.976	0.000	Pass	Pass	0.653	1,449.036
69	0.80	1.75	0.368	0.142	0.114	0.431	7,612.971	0.000	Pass	Pass	0.653	1,259.768
70	0.80	1.75	0.389	0.145	0.099	0.451	7,615.018	0.000	Pass	Pass	0.651	1,700.227
71	0.80	1.75	0.389	0.143	0.107	0.441	7,613.600	0.000	Pass	Pass	0.650	1,457.430
72	0.80	1.75	0.389	0.140	0.114	0.431	7,612.581	0.000	Pass	Pass	0.651	1,267.423
73	0.55	2.00	0.339	0.147	0.099	0.464	7,616.146	0.000	Pass	Pass	0.586	1,986.271
74	0.55	2.00	0.339	0.144	0.107	0.454	7,614.084	0.001	Pass	Pass	0.586	1,706.017
75	0.55	2.00	0.339	0.142	0.114	0.445	7,612.535	0.000	Pass	Pass	0.586	1,485.475
76	0.55	2.00	0.368	0.144	0.099	0.464	7,615.286	0.000	Pass	Pass	0.588	2,001.595
77	0.55	2.00	0.368	0.141	0.107	0.455	7,613.193	0.000	Pass	Pass	0.588	1,719.681
78	0.55	2.00	0.368	0.139	0.114	0.446	7,611.618	0.000	Pass	Pass	0.588	1,497.846
79	0.55	2.00	0.389	0.142	0.099	0.465	7,614.984	0.000	Pass	Pass	0.590	2,012.396
80	0.55	2.00	0.389	0.140	0.107	0.455	7,612.866	0.000	Pass	Pass	0.590	1,729.365
81	0.55	2.00	0.389	0.137	0.114	0.446	7,611.272	0.000	Pass	Pass	0.590	1,506.646
82	0.63	2.00	0.339	0.148	0.099	0.460	7,614.689	0.001	Pass	Pass	0.587	1,861.355
83	0.63	2.00	0.339	0.145	0.107	0.450	7,612.913	0.000	Pass	Pass	0.587	1,598.024
84	0.63	2.00	0.339	0.142	0.114	0.441	7,611.599	0.000	Pass	Pass	0.587	1,391.180
85	0.63	2.00	0.368	0.145	0.099	0.460	7,613.783	0.000	Pass	Pass	0.590	1,875.536
86	0.63	2.00	0.368	0.142	0.107	0.450	7,611.979	0.000	Pass	Pass	0.590	1,610.721
87	0.63	2.00	0.368	0.140	0.114	0.441	7,610.644	0.000	Pass	Pass	0.590	1,402.708
88	0.63	2.00	0.389	0.143	0.099	0.461	7,613.449	0.003	Pass	Pass	0.592	1,885.579
89	0.63	2.00	0.389	0.140	0.107	0.451	7,611.624	0.000	Pass	Pass	0.592	1,619.758
90	0.63	2.00	0.389	0.138	0.114	0.442	7,610.271	0.000	Pass	Pass	0.592	1,410.939
91	0.72	2.00	0.339	0.148	0.099	0.456	7,615.200	0.000	Pass	Pass	0.596	1,787.393
92	0.72	2.00	0.339	0.145	0.107	0.446	7,613.616	0.000	Pass	Pass	0.595	1,535.979
93	0.72	2.00	0.339	0.143	0.114	0.437	7,612.457	0.000	Pass	Pass	0.595	1,338.497
94	0.72	2.00	0.368	0.145	0.099	0.456	7,614.261	0.000	Pass	Pass	0.592	1,800.720
95	0.72	2.00	0.368	0.142	0.107	0.446	7,612.654	0.000	Pass	Pass	0.592	1,547.985
96	0.72	2.00	0.368	0.140	0.114	0.437	7,611.475	0.000	Pass	Pass	0.592	1,349.449
97	0.72	2.00	0.389	0.143	0.099	0.457	7,613.905	0.000	Pass	Pass	0.592	1,810.205
98	0.72	2.00	0.389	0.140	0.107	0.447	7,612.279	0.000	Pass	Pass	0.592	1,556.566
99	0.72	2.00	0.389	0.138	0.114	0.437	7,611.084	0.000	Pass	Pass	0.592	1,357.297
100	0.80	2.00	0.339	0.148	0.099	0.453	7,616.732	0.000	Pass	Pass	0.647	1,769.394
101	0.80	2.00	0.339	0.145	0.107	0.443	7,615.220	0.000	Pass	Pass	0.646	1,522.210
102	0.80	2.00	0.339	0.143	0.114	0.434	7,614.114	0.000	Pass	Pass	0.646	1,327.819
103	0.80	2.00	0.368	0.145	0.099	0.453	7,615.777	0.000	Pass	Pass	0.643	1,782.371
104	0.80	2.00	0.368	0.143	0.107	0.443	7,614.245	0.000	Pass	Pass	0.642	1,533.969
105	0.80	2.00	0.368	0.140	0.114	0.434	7,613.122	0.000	Pass	Pass	0.642	1,338.594
106	0.80	2.00	0.389	0.143	0.099	0.454	7,615.410	0.000	Pass	Pass	0.640	1,791.633
107	0.80	2.00	0.389	0.141	0.107	0.443	7,613.862	0.000	Pass	Pass	0.640	1,542.392
108	0.80	2.00	0.389	0.138	0.114	0.434	7,612.724	0.000	Pass	Pass	0.640	1,346.329

Show code

ggplot(m_posterior_df, aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = fit_expected_color, alpha = 0.32, color = "white") +
  geom_vline(
    data = natural_mortality_grid_values,
    aes(xintercept = value),
    color = "#D55E00",
    linewidth = 0.7,
    linetype = "dashed"
  ) +
  facet_wrap(vars(parameter), scales = "free_x") +
  labs(x = "Natural mortality", y = "Posterior density") +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.03))) +
  scale_y_zero()

Figure 10: Combined length-based MCMC posterior for derived M0 and estimated M10, with the direct-get_M grid values overlaid. The direct grid uses the 1%, 50%, and 95% M0 quantiles and retains the existing posterior-informed M10 values; it fixes M0 and M10 while estimating M4 and M30.

Show code

mle_grid_freq |>
  filter(Freq > 0) |>
  transmute(
    Cell = format_decimal(Cell, digits = 0),
    h,
    psi,
    `Fixed M0` = fixed_m0,
    `Fixed M10` = fixed_m10,
    `Delta NLL` = delta_nll,
    Probability = prob,
    Frequency = format_decimal(Freq, digits = 0)
  ) |>
  mutate(
    h = format_decimal(h, digits = 2),
    psi = format_decimal(psi, digits = 2)
  ) |>
  mutate(across(where(is.numeric), ~ format_decimal(.x, digits = 3))) |>
  replace_missing() |>
  kable(caption = paste0("Full-objective-weighted resampling frequencies for ", mle_grid_resample_n, " draws from the 108 distinct direct-get_M fits."))

Table 9: Full-objective-weighted resampling frequencies for 2000 draws from the 108 distinct direct-get_M fits.

Cell	h	psi	Fixed M0	Fixed M10	Delta NLL	Probability	Frequency
1	0.55	1.50	0.339	0.099	7.094	0.000	3
2	0.55	1.50	0.339	0.107	5.309	0.002	3
3	0.55	1.50	0.339	0.114	4.010	0.007	19
5	0.55	1.50	0.368	0.107	4.438	0.005	9
6	0.55	1.50	0.368	0.114	3.113	0.017	37
7	0.55	1.50	0.389	0.099	5.965	0.001	1
8	0.55	1.50	0.389	0.107	4.126	0.006	10
9	0.55	1.50	0.389	0.114	2.782	0.024	45
10	0.63	1.50	0.339	0.099	5.640	0.000	1
11	0.63	1.50	0.339	0.107	4.140	0.002	6
12	0.63	1.50	0.339	0.114	3.075	0.006	11
13	0.63	1.50	0.368	0.099	4.752	0.001	1
14	0.63	1.50	0.368	0.107	3.225	0.005	12
15	0.63	1.50	0.368	0.114	2.138	0.016	22
16	0.63	1.50	0.389	0.099	4.432	0.002	2
17	0.63	1.50	0.389	0.107	2.884	0.007	14
18	0.63	1.50	0.389	0.114	1.780	0.023	39
19	0.72	1.50	0.339	0.099	6.187	0.001	3
20	0.72	1.50	0.339	0.107	4.874	0.002	4
21	0.72	1.50	0.339	0.114	3.959	0.006	11
22	0.72	1.50	0.368	0.099	5.265	0.002	6
23	0.72	1.50	0.368	0.107	3.929	0.006	5
24	0.72	1.50	0.368	0.114	2.996	0.015	22
25	0.72	1.50	0.389	0.099	4.922	0.002	2
26	0.72	1.50	0.389	0.107	3.569	0.008	18
27	0.72	1.50	0.389	0.114	2.620	0.021	40
28	0.80	1.50	0.339	0.099	7.753	0.001	1
29	0.80	1.50	0.339	0.107	6.506	0.002	5
30	0.80	1.50	0.339	0.114	5.640	0.005	10
31	0.80	1.50	0.368	0.099	6.816	0.002	1
32	0.80	1.50	0.368	0.107	5.550	0.006	7
33	0.80	1.50	0.368	0.114	4.667	0.014	22
34	0.80	1.50	0.389	0.099	6.464	0.002	7
35	0.80	1.50	0.389	0.107	5.181	0.009	19
36	0.80	1.50	0.389	0.114	4.284	0.021	50
38	0.55	1.75	0.339	0.107	3.674	0.003	8
39	0.55	1.75	0.339	0.114	2.247	0.014	32
40	0.55	1.75	0.368	0.099	4.749	0.001	2
41	0.55	1.75	0.368	0.107	2.792	0.008	10
42	0.55	1.75	0.368	0.114	1.340	0.035	70
43	0.55	1.75	0.389	0.099	4.454	0.002	5
44	0.55	1.75	0.389	0.107	2.473	0.011	30
45	0.55	1.75	0.389	0.114	1.002	0.049	106
46	0.63	1.75	0.339	0.099	4.144	0.001	2
47	0.63	1.75	0.339	0.107	2.503	0.004	6
48	0.63	1.75	0.339	0.114	1.311	0.013	23
49	0.63	1.75	0.368	0.099	3.246	0.002	3
50	0.63	1.75	0.368	0.107	1.579	0.010	20
51	0.63	1.75	0.368	0.114	0.365	0.033	57
52	0.63	1.75	0.389	0.099	2.920	0.003	4
53	0.63	1.75	0.389	0.107	1.231	0.014	24
54	0.63	1.75	0.389	0.114	0.000	0.047	106
55	0.72	1.75	0.339	0.099	4.673	0.001	2
56	0.72	1.75	0.339	0.107	3.222	0.004	9
57	0.72	1.75	0.339	0.114	2.183	0.012	18
58	0.72	1.75	0.368	0.099	3.742	0.002	8
59	0.72	1.75	0.368	0.107	2.269	0.011	20
60	0.72	1.75	0.368	0.114	1.211	0.031	63
61	0.72	1.75	0.389	0.099	3.393	0.003	7
62	0.72	1.75	0.389	0.107	1.901	0.015	29
63	0.72	1.75	0.389	0.114	0.827	0.045	109
64	0.80	1.75	0.339	0.099	6.223	0.001	1
65	0.80	1.75	0.339	0.107	4.842	0.004	14
66	0.80	1.75	0.339	0.114	3.855	0.011	24
67	0.80	1.75	0.368	0.099	5.278	0.003	6
68	0.80	1.75	0.368	0.107	3.877	0.011	24
69	0.80	1.75	0.368	0.114	2.872	0.030	52
70	0.80	1.75	0.389	0.099	4.918	0.004	9
71	0.80	1.75	0.389	0.107	3.500	0.016	36
72	0.80	1.75	0.389	0.114	2.481	0.045	91
74	0.55	2.00	0.339	0.107	3.984	0.002	2
75	0.55	2.00	0.339	0.114	2.435	0.007	12
76	0.55	2.00	0.368	0.099	5.186	0.000	1
77	0.55	2.00	0.368	0.107	3.093	0.004	10
78	0.55	2.00	0.368	0.114	1.518	0.018	29
79	0.55	2.00	0.389	0.099	4.884	0.001	1
80	0.55	2.00	0.389	0.107	2.767	0.005	7
81	0.55	2.00	0.389	0.114	1.172	0.026	55
82	0.63	2.00	0.339	0.099	4.590	0.000	1
83	0.63	2.00	0.339	0.107	2.813	0.002	4
84	0.63	2.00	0.339	0.114	1.500	0.006	14
85	0.63	2.00	0.368	0.099	3.683	0.001	1
86	0.63	2.00	0.368	0.107	1.880	0.004	13
87	0.63	2.00	0.368	0.114	0.544	0.017	32
88	0.63	2.00	0.389	0.099	3.349	0.001	1
89	0.63	2.00	0.389	0.107	1.525	0.006	16
90	0.63	2.00	0.389	0.114	0.171	0.025	47
92	0.72	2.00	0.339	0.107	3.517	0.002	2
93	0.72	2.00	0.339	0.114	2.357	0.006	8
94	0.72	2.00	0.368	0.099	4.161	0.001	2
95	0.72	2.00	0.368	0.107	2.554	0.005	9
96	0.72	2.00	0.368	0.114	1.375	0.016	40
97	0.72	2.00	0.389	0.099	3.805	0.001	1
98	0.72	2.00	0.389	0.107	2.179	0.007	17
99	0.72	2.00	0.389	0.114	0.984	0.024	57
100	0.80	2.00	0.339	0.099	6.632	0.000	1
101	0.80	2.00	0.339	0.107	5.121	0.002	2
102	0.80	2.00	0.339	0.114	4.015	0.006	12
103	0.80	2.00	0.368	0.099	5.677	0.001	3
104	0.80	2.00	0.368	0.107	4.145	0.005	19
105	0.80	2.00	0.368	0.114	3.022	0.016	25
106	0.80	2.00	0.389	0.099	5.311	0.002	2
107	0.80	2.00	0.389	0.107	3.762	0.007	17
108	0.80	2.00	0.389	0.114	2.624	0.023	39

Show code

mle_grid_sampled_cells <- tibble(
  Draw = seq_along(mle_grid_sample$grid_cells),
  Cell = mle_grid_sample$grid_cells
) |>
  left_join(
    m10_mle_summary |>
      select(Cell, h, psi, fixed_m0, fixed_m10),
    by = "Cell"
  )

if (length(mle_grid_tmbfit$warmup) != 1L ||
    !is.finite(mle_grid_tmbfit$warmup) ||
    mle_grid_tmbfit$warmup < 0L ||
    mle_grid_tmbfit$warmup != as.integer(mle_grid_tmbfit$warmup) ||
    dim(mle_grid_tmbfit$samples)[1L] !=
      mle_grid_tmbfit$warmup + length(mle_grid_sample$grid_cells)) {
  stop("The canonical MLE-grid tmbfit does not align with sampled cells.",
       call. = FALSE)
}
mle_grid_post_iterations <- seq.int(
  as.integer(mle_grid_tmbfit$warmup) + 1L,
  dim(mle_grid_tmbfit$samples)[1L]
)
mle_grid_post <- mle_grid_tmbfit$samples[
  mle_grid_post_iterations, 1L, , drop = FALSE
][, 1L, ]
lev_df <- tibble(
  h = exp(mle_grid_post[, "par_log_h"]),
  M0 = exp(mle_grid_post[, "par_log_m0"]),
  M10 = exp(mle_grid_post[, "par_log_m10"]),
  psi = exp(mle_grid_post[, "par_log_psi"])
)
expected_lev_df <- mle_grid_sampled_cells |>
  transmute(h, M0 = fixed_m0, M10 = fixed_m10, psi)
if (!isTRUE(all.equal(
      unclass(lev_df),
      unclass(expected_lev_df),
      check.attributes = FALSE,
      tolerance = 1e-12
    ))) {
  stop("The MLE-grid tmbfit parameters do not match the sampled cells.",
       call. = FALSE)
}

plot_levs(
  lev_df,
  lev_vars,
  point_alpha = 0.18,
  point_size = 0.25,
  discrete_jitter = 0.42,
  discrete_levels = list(
    h = mle_grid_h_values,
    M0 = m0_grid_values,
    M10 = m10_grid_values,
    psi = mle_grid_psi_values
  )
)

Figure 11: Full-objective-weighted resampling from the 108 distinct direct-get_M fits. The projection-ready tmbfit retains every fitted parameter and all four fixed grid coordinates. Diagonal panels show marginal counts or histograms; off-diagonal panels show jittered sampled grid-cell pairs from 2,000 resampled rows.

Show code

mle_grid_sampled_biomass <- mle_grid_sampled_cells |>
  select(Draw, Cell) |>
  left_join(m10_mle_biomass, by = "Cell", relationship = "many-to-many")

mle_grid_biomass_summary <- mle_grid_sampled_biomass |>
  group_by(Year) |>
  summarise(
    lower = quantile(relative_spawning_biomass, 0.025, na.rm = TRUE),
    median = median(relative_spawning_biomass, na.rm = TRUE),
    upper = quantile(relative_spawning_biomass, 0.975, na.rm = TRUE),
    .groups = "drop"
  )
mle_grid_comparison_terminal_year <- max(intersect(
  combined_biomass_2000_summary$Year,
  mle_grid_biomass_summary$Year
))
mcmc_grid_terminal_biomass <- combined_biomass_2000_summary |>
  filter(Year == mle_grid_comparison_terminal_year) |>
  slice(1)
mle_grid_terminal_biomass <- mle_grid_biomass_summary |>
  filter(Year == mle_grid_comparison_terminal_year) |>
  slice(1)

Show code

mortality_ages <- seq.int(
  base_fit$data$min_age,
  base_fit$data$max_age
)
required_mcmc_mortality_parameters <- c(
  "par_log_m10", "par_log_m30"
)
required_mle_mortality_parameters <- c(
  "par_log_m0", "par_log_m4", "par_log_m10", "par_log_m30"
)
if (!all(required_mcmc_mortality_parameters %in%
      names(combined_grid_post)) ||
    !all(required_mle_mortality_parameters %in%
      colnames(mle_grid_post)) ||
    nrow(combined_grid_post) != combined_mcmc_draws ||
    nrow(mle_grid_post) != mle_grid_resample_n ||
    !identical(as.integer(base_fit$data$M_switch), 2L) ||
    !identical(as.integer(m10_mle_data$M_switch), 1L)) {
  stop(
    "The MCMC- and MLE-grid mortality draws do not match their approved ",
    "parameterizations.",
    call. = FALSE
  )
}

mcmc_mc <- if ("par_mc" %in% names(combined_grid_post)) {
  as.numeric(combined_grid_post$par_mc)
} else {
  rep(
    as.numeric(base_fit$parameters$par_mc),
    nrow(combined_grid_post)
  )
}
if (length(mcmc_mc) != nrow(combined_grid_post) ||
    any(!is.finite(mcmc_mc))) {
  stop("The MCMC-grid mortality exponent is invalid.", call. = FALSE)
}

mcmc_grid_mortality <- vapply(
  seq_len(nrow(combined_grid_post)),
  function(draw) {
    as.numeric(sbt::get_M_length(
      min_age = base_fit$data$min_age,
      max_age = base_fit$data$max_age,
      age_increase_M = base_fit$data$age_increase_M,
      m10 = exp(combined_grid_post$par_log_m10[[draw]]),
      m30 = exp(combined_grid_post$par_log_m30[[draw]]),
      mc = mcmc_mc[[draw]],
      length_mu_ysa = base_fit$data$length_mu_ysa
    ))
  },
  numeric(length(mortality_ages))
)
mle_grid_mortality <- vapply(
  seq_len(nrow(mle_grid_post)),
  function(draw) {
    as.numeric(sbt::get_M(
      min_age = m10_mle_data$min_age,
      max_age = m10_mle_data$max_age,
      age_increase_M = m10_mle_data$age_increase_M,
      m0 = exp(mle_grid_post[draw, "par_log_m0"]),
      m4 = exp(mle_grid_post[draw, "par_log_m4"]),
      m10 = exp(mle_grid_post[draw, "par_log_m10"]),
      m30 = exp(mle_grid_post[draw, "par_log_m30"])
    ))
  },
  numeric(length(mortality_ages))
)
if (any(!is.finite(mcmc_grid_mortality)) ||
    any(mcmc_grid_mortality <= 0) ||
    any(!is.finite(mle_grid_mortality)) ||
    any(mle_grid_mortality <= 0)) {
  stop("The grid mortality comparison contains invalid values.",
       call. = FALSE)
}

mortality_draw_frame <- function(values, source) {
  tibble(
    source = source,
    Draw = rep(seq_len(ncol(values)), each = nrow(values)),
    Age = rep(mortality_ages, times = ncol(values)),
    M = as.vector(values)
  )
}
mortality_comparison_summary <- bind_rows(
  mortality_draw_frame(
    mcmc_grid_mortality,
    "MCMC grid (length-based)"
  ),
  mortality_draw_frame(
    mle_grid_mortality,
    "MLE grid (direct)"
  )
) |>
  group_by(source, Age) |>
  summarise(
    lower = quantile(M, 0.025),
    median = median(M),
    upper = quantile(M, 0.975),
    .groups = "drop"
  )

Show code

mortality_comparison_summary |>
  ggplot(aes(x = Age, y = median, color = source, fill = source)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.18, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_color_manual(values = c(
    "MCMC grid (length-based)" = fit_expected_color,
    "MLE grid (direct)" = "#D55E00"
  )) +
  scale_fill_manual(values = c(
    "MCMC grid (length-based)" = fit_expected_color,
    "MLE grid (direct)" = "#D55E00"
  )) +
  labs(x = "Age", y = "Natural mortality (M)", color = NULL, fill = NULL) +
  scale_x_continuous(
    breaks = pretty_breaks(n = 8),
    limits = range(mortality_ages),
    expand = expansion(mult = c(0, 0.02))
  ) +
  scale_y_zero()

Figure 12: Natural mortality at age comparison between the combined length-based MCMC grid posterior and the full-objective-weighted 108-cell direct-get_M MLE grid. Lines show medians and shaded ribbons show 95% intervals across the respective 2,000-draw samples.

Show code

bind_rows(
  combined_biomass_2000_summary |> mutate(source = "MCMC grid"),
  mle_grid_biomass_summary |> mutate(source = "MLE grid")
) |>
  ggplot(aes(x = Year, y = median, color = source, fill = source)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.18, color = NA) +
  geom_line(linewidth = 0.75) +
  scale_color_manual(values = c("MCMC grid" = fit_expected_color, "MLE grid" = "#D55E00")) +
  scale_fill_manual(values = c("MCMC grid" = fit_expected_color, "MLE grid" = "#D55E00")) +
  labs(x = "Year", y = "Relative TRO", color = NULL, fill = NULL) +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero()

Figure 13: Relative total reproductive output comparison between the combined MCMC grid posterior and the full-objective-weighted 108-cell MLE grid. Lines show medians and shaded ribbons show 95% intervals across 2,000 sampled draws or grid cells.

In 2026, the MCMC-grid relative-TRO median is 0.281 (95% interval 0.209–0.373), whereas the full-objective-weighted MLE-grid median is 0.280 (95% interval 0.228–0.329). The deterministic MLE-grid resample is therefore slightly lower and narrower at the terminal endpoint; it is a robustness comparison, not a replacement for the MCMC posterior.

Previous and current CCSBT base-assessment grids

The previous CCSBT base stock assessment used a 108-cell deterministic grid with four steepness values, three psi values, three M0 values, and three M10 values. The comparison below reads all 108 archived ADMB cell reports directly from the package source and shows them against all 108 cells of the current MLE grid. Neither grid is likelihood weighted in this comparison.

Show code

sbt_source_root <- normalizePath(
  file.path(dirname(dirname(esc_dir)), "sbt"),
  mustWork = TRUE
)
previous_mle_grid_root <- file.path(
  sbt_source_root,
  "data-raw", "base22", "base22sqrt"
)
previous_mle_grid_files <- sort(list.files(
  previous_mle_grid_root,
  pattern = "_lab\\.rep$",
  full.names = TRUE
))
if (length(previous_mle_grid_files) != 108L) {
  stop(
    "The previous CCSBT base-assessment grid must contain all 108 ADMB ",
    "cell reports.",
    call. = FALSE
  )
}

read_admb_numeric_field <- function(lines, field, path) {
  marker <- paste0("$", field)
  start <- which(lines == marker)
  if (length(start) != 1L) {
    stop(
      "Expected one `", marker, "` field in ", basename(path), ".",
      call. = FALSE
    )
  }
  next_marker <- which(
    seq_along(lines) > start &
      grepl("^\\$", lines)
  )
  end <- if (length(next_marker)) {
    next_marker[[1L]] - 1L
  } else {
    length(lines)
  }
  values <- scan(
    text = paste(lines[seq.int(start + 1L, end)], collapse = " "),
    quiet = TRUE
  )
  if (!length(values) || any(!is.finite(values))) {
    stop(
      "Invalid `", marker, "` field in ", basename(path), ".",
      call. = FALSE
    )
  }
  values
}

read_previous_mle_grid_cell <- function(path, cell) {
  lines <- trimws(readLines(path, warn = FALSE))
  assessment_years <- read_admb_numeric_field(lines, "years", path)
  B0 <- read_admb_numeric_field(lines, "B0", path)
  spawning_biomass <- read_admb_numeric_field(lines, "Sbio", path)
  mortality <- read_admb_numeric_field(lines, "M", path)
  h <- read_admb_numeric_field(lines, "steep", path)
  psi <- read_admb_numeric_field(lines, "psi", path)
  if (length(assessment_years) != 2L ||
      length(B0) != 1L || B0 <= 0 ||
      length(mortality) != 31L ||
      length(h) != 1L || length(psi) != 1L ||
      length(spawning_biomass) !=
        diff(as.integer(assessment_years)) + 2L ||
      any(spawning_biomass <= 0)) {
    stop(
      "The archived grid cell has an invalid population-state layout: ",
      basename(path), ".",
      call. = FALSE
    )
  }
  tibble(
    Cell = as.integer(cell),
    h = h,
    psi = psi,
    M0 = mortality[[1L]],
    M10 = mortality[[11L]],
    Year = seq.int(
      as.integer(assessment_years[[1L]]),
      length.out = length(spawning_biomass)
    ),
    relative_spawning_biomass = spawning_biomass / B0
  )
}

previous_mle_grid_biomass <- imap_dfr(
  previous_mle_grid_files,
  ~ read_previous_mle_grid_cell(.x, .y)
)
previous_mle_grid_coordinates <- previous_mle_grid_biomass |>
  distinct(Cell, h, psi, M0, M10)
expected_previous_mle_grid_coordinates <- list(
  h = c(0.55, 0.63, 0.72, 0.80),
  psi = c(1.50, 1.75, 2.00),
  M0 = c(0.40, 0.45, 0.50),
  M10 = c(0.065, 0.085, 0.105)
)
if (nrow(previous_mle_grid_coordinates) != 108L ||
    any(vapply(
      names(expected_previous_mle_grid_coordinates),
      function(parameter) {
        !isTRUE(all.equal(
          sort(unique(previous_mle_grid_coordinates[[parameter]])),
          expected_previous_mle_grid_coordinates[[parameter]],
          tolerance = 1e-12,
          check.attributes = FALSE
        ))
      },
      logical(1)
    ))) {
  stop(
    "The archived ADMB reports do not reconstruct the approved previous ",
    "108-cell CCSBT grid.",
    call. = FALSE
  )
}

previous_current_mle_grid_biomass <- bind_rows(
  previous_mle_grid_biomass |>
    transmute(
      Grid = "Previous CCSBT base assessment",
      Cell,
      Year,
      relative_spawning_biomass
    ),
  m10_mle_biomass |>
    transmute(
      Grid = "Current ESC31 assessment",
      Cell,
      Year,
      relative_spawning_biomass
    )
)
if (n_distinct(
      previous_current_mle_grid_biomass$Cell[
        previous_current_mle_grid_biomass$Grid ==
          "Previous CCSBT base assessment"
      ]
    ) != 108L ||
    n_distinct(
      previous_current_mle_grid_biomass$Cell[
        previous_current_mle_grid_biomass$Grid ==
          "Current ESC31 assessment"
      ]
    ) != 108L ||
    any(!is.finite(
      previous_current_mle_grid_biomass$relative_spawning_biomass
    ))) {
  stop("The all-cell MLE-grid comparison is incomplete.", call. = FALSE)
}
previous_current_mle_grid_summary <-
  previous_current_mle_grid_biomass |>
  group_by(Grid, Year) |>
  summarise(
    median = median(relative_spawning_biomass),
    .groups = "drop"
  )

format_grid_values <- function(x, digits) {
  paste(
    formatC(
      sort(unique(as.numeric(x))),
      format = "f",
      digits = digits
    ),
    collapse = ", "
  )
}
previous_current_mle_grid_coordinates <- bind_rows(
  previous_mle_grid_coordinates |>
    summarise(
      Grid = "Previous CCSBT base assessment",
      Cells = n(),
      h = format_grid_values(h, 2),
      psi = format_grid_values(psi, 2),
      M0 = format_grid_values(M0, 3),
      M10 = format_grid_values(M10, 3),
      `Final population state` = max(previous_mle_grid_biomass$Year)
    ),
  m10_mle_grid |>
    summarise(
      Grid = "Current ESC31 assessment",
      Cells = n(),
      h = format_grid_values(h, 2),
      psi = format_grid_values(psi, 2),
      M0 = format_grid_values(m0, 3),
      M10 = format_grid_values(m10, 3),
      `Final population state` = max(m10_mle_biomass$Year)
    )
)

Show code

previous_current_mle_grid_coordinates |>
  mutate(
    Cells = format_decimal(Cells, digits = 0),
    `Final population state` =
      format_decimal(`Final population state`, digits = 0)
  ) |>
  kable()

Table 10: All-cell design comparison for the previous CCSBT base-assessment grid and the current ESC31 MLE grid.

Grid	Cells	h	psi	M0	M10	Final population state
Previous CCSBT base assessment	108	0.55, 0.63, 0.72, 0.80	1.50, 1.75, 2.00	0.400, 0.450, 0.500	0.065, 0.085, 0.105	2,023
Current ESC31 assessment	108	0.55, 0.63, 0.72, 0.80	1.50, 1.75, 2.00	0.339, 0.368, 0.389	0.099, 0.107, 0.114	2,026

Show code

previous_current_mle_grid_biomass |>
  ggplot(aes(
    x = Year,
    y = relative_spawning_biomass,
    color = Grid,
    group = interaction(Grid, Cell)
  )) +
  geom_line(alpha = 0.14, linewidth = 0.28) +
  geom_line(
    data = previous_current_mle_grid_summary,
    aes(x = Year, y = median, color = Grid, group = Grid),
    inherit.aes = FALSE,
    linewidth = 1.05
  ) +
  scale_color_manual(values = c(
    "Previous CCSBT base assessment" = "#7A7A7A",
    "Current ESC31 assessment" = fit_expected_color
  )) +
  labs(x = "Year", y = "Relative TRO", color = NULL) +
  scale_x_continuous(breaks = pretty_breaks(n = 8)) +
  scale_y_zero() +
  theme(legend.position = "bottom")

Figure 14: Relative total reproductive output for all 108 cells in the previous CCSBT base-assessment grid and all 108 cells in the current ESC31 MLE grid. Thin lines are individual unweighted MLE cells and thick lines are the unweighted median within each grid. The previous grid reports population states through 2023 and the current grid through 2026.

Appendix: MCMC Diagnostics

Sampler Parameters

Show code

plot_sampler_cell <- function(cell) {
  row <- grid_diagnostics |>
    filter(Cell == cell) |>
    slice(1)

  if (!nrow(row)) {
    stop("Grid cell ", cell, " is not available in the completed grid diagnostics.", call. = FALSE)
  }

  cell_output <- load_grid_output(row$file, fallback_cell = cell)
  plot_sampler_params(fit = cell_output$mcmc, plot = TRUE)
}

Show code

plot_sampler_cell(1)

Figure 15: Sampler-parameter diagnostics for grid cell 1.

Show code

plot_sampler_cell(2)

Figure 16: Sampler-parameter diagnostics for grid cell 2.

Show code

plot_sampler_cell(3)

Figure 17: Sampler-parameter diagnostics for grid cell 3.

Show code

plot_sampler_cell(4)

Figure 18: Sampler-parameter diagnostics for grid cell 4.

Show code

plot_sampler_cell(5)

Figure 19: Sampler-parameter diagnostics for grid cell 5.

Show code

plot_sampler_cell(6)

Figure 20: Sampler-parameter diagnostics for grid cell 6.

Show code

plot_sampler_cell(7)

Figure 21: Sampler-parameter diagnostics for grid cell 7.

Show code

plot_sampler_cell(8)

Figure 22: Sampler-parameter diagnostics for grid cell 8.

Show code

plot_sampler_cell(9)

Figure 23: Sampler-parameter diagnostics for grid cell 9.

MCMC Pair Plots

Show code

plot_pairs_cell <- function(cell, order = "slow", pars = 1:5) {
  row <- grid_diagnostics |>
    filter(Cell == cell) |>
    slice(1)

  if (!nrow(row)) {
    stop("Grid cell ", cell, " is not available in the completed grid diagnostics.", call. = FALSE)
  }

  cell_output <- load_grid_output(row$file, fallback_cell = cell)
  pairs_rtmb(
    fit = cell_output$mcmc,
    order = order,
    pars = pars,
    diag = "trace"
  )
}

Show code

plot_pairs_cell(1)

Figure 24: pairs_rtmb diagnostic plot for grid cell 1, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(2)

Figure 25: pairs_rtmb diagnostic plot for grid cell 2, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(3)

Figure 26: pairs_rtmb diagnostic plot for grid cell 3, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(4)

Figure 27: pairs_rtmb diagnostic plot for grid cell 4, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(5)

Figure 28: pairs_rtmb diagnostic plot for grid cell 5, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(6)

Figure 29: pairs_rtmb diagnostic plot for grid cell 6, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(7)

Figure 30: pairs_rtmb diagnostic plot for grid cell 7, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(8)

Figure 31: pairs_rtmb diagnostic plot for grid cell 8, showing the five slowest-mixing parameters.

Show code

plot_pairs_cell(9)

Figure 32: pairs_rtmb diagnostic plot for grid cell 9, showing the five slowest-mixing parameters.

References

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.

Monnahan, Cole C. 2024. “Toward Good Practices for Bayesian Data-Rich Fisheries Stock Assessments Using a Modern Statistical Workflow.” Fisheries Research 275: 107024. https://doi.org/10.1016/j.fishres.2024.107024.

Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved R-Hat for Assessing Convergence of MCMC.” Bayesian Analysis 16 (2): 667–718. https://doi.org/10.1214/20-BA1221.