---
title: "ESC31 Grid"
format:
  html:
    toc: true
    page-layout: full
    theme: default
    code-fold: true
    code-summary: "Show code"
    lightbox: true
    mainfont: system-ui
    css: esc31.css
    favicon: favicon.ico
    include-in-header:
      text: |
        <link rel="icon" href="favicon.ico" sizes="any">
    include-before-body: nav.html
    embed-resources: true
execute:
  warning: false
  message: false
bibliography: references.bib
link-citations: true
editor_options:
  chunk_output_type: console
---

```{r}
#| label: setup
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))
  )
}
```

::: {.callout-important}
## 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
[@HillaryEtAl2023StockAssessment]. 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

```{r}
#| label: paths
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

```{r}
#| label: tbl-approved-grid-cells
#| tbl-cap: "Approved nine-cell MCMC grid. Cell 5 is the accepted base posterior; the other eight cells require a new MLE and MCMC."
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)
```

```{r}
#| label: tbl-approved-grid-sampling
#| tbl-cap: "Approved mixed-source grid sampling and downstream design."
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
  )
```

## Load Base Fit

```{r}
#| label: load-base-fit
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

```{r}
#| label: helper-functions
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.

```{r}
#| label: run-grid-mcmc
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

```{r}
#| label: read-grid-summary
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)
  )
```

```{r}
#| label: tbl-grid-status
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),
    "."
  ))
```

For production inference, the acceptance limits used in
@tbl-grid-diagnostics 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 [@Monnahan2024]. These limits follow the rank-normalized
$\widehat{R}$ and bulk/tail ESS guidance of @VehtariEtAl2021, 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.

```{r}
#| label: tbl-grid-diagnostics
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."
  )
```

<!--
The archived failure-only MLE-grid diagnostic remains executable below for a
failed-grid checkpoint, but it is not part of the accepted-results page
structure and is never used downstream.
-->

```{r}
#| label: load-diagnostic-mle-grid
#| echo: false
failed_grid_cells <- grid_diagnostic_checks |>
  filter(!diagnostics_pass) |>
  pull(Cell)

if (length(failed_grid_cells)) {
diagnostic_mle_grid_file <- file.path(
  esc_dir,
  "archive",
  "exploratory",
  "diagnostic_mle_grid_failed_mcmc_20260723",
  "diagnostic_mle_grid_108.rds"
)
diagnostic_mle_grid <- read_rds_cache(diagnostic_mle_grid_file)
diagnostic_mle_grid_fields <- c(
  "status", "reason", "source_grid_diagnostics", "failed_source_cells",
  "grid", "completed_cells", "error_cells", "failed_gate_cells", "summary",
  "biomass", "state_records", "errors"
)
diagnostic_mle_grid_diagnostic_fields <- c(
  "Cell", "h", "psi", "convergence", "objective", "estimable", "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_mle_grid_biomass_fields <- c(
  "Cell", "h", "psi", "fixed_m0", "fixed_m10", "Year",
  "relative_spawning_biomass"
)

diagnostic_mle_grid_current <-
  is.list(diagnostic_mle_grid) &&
  all(diagnostic_mle_grid_fields %in% names(diagnostic_mle_grid)) &&
  identical(diagnostic_mle_grid$status, "diagnostic_non_production") &&
  identical(
    diagnostic_mle_grid$reason,
    "source_grid_failed_production_rhat_gate"
  ) &&
  identical(
    as.integer(diagnostic_mle_grid$failed_source_cells),
    as.integer(failed_grid_cells)
  ) &&
  is.data.frame(diagnostic_mle_grid$source_grid_diagnostics) &&
  all(diagnostic_mle_grid_diagnostic_fields %in%
        names(diagnostic_mle_grid$source_grid_diagnostics)) &&
  isTRUE(all.equal(
    diagnostic_mle_grid$source_grid_diagnostics[
      diagnostic_mle_grid_diagnostic_fields
    ],
    as.data.frame(grid_diagnostics)[
      diagnostic_mle_grid_diagnostic_fields
    ],
    check.attributes = FALSE,
    tolerance = 1e-10
  )) &&
  identical(as.integer(diagnostic_mle_grid$completed_cells), 1:108) &&
  !length(diagnostic_mle_grid$error_cells) &&
  !length(diagnostic_mle_grid$failed_gate_cells) &&
  !length(diagnostic_mle_grid$errors) &&
  is.data.frame(diagnostic_mle_grid$summary) &&
  nrow(diagnostic_mle_grid$summary) == 108L &&
  all(diagnostic_mle_grid$summary$diagnostic_passes) &&
  is.data.frame(diagnostic_mle_grid$biomass) &&
  identical(
    names(diagnostic_mle_grid$biomass),
    diagnostic_mle_grid_biomass_fields
  ) &&
  nrow(diagnostic_mle_grid$biomass) ==
    108L * (as.integer(base_fit$data$n_year) + 1L) &&
  all(is.finite(as.matrix(diagnostic_mle_grid$biomass)))

if (!isTRUE(diagnostic_mle_grid_current)) {
  stop(
    "The diagnostic 108-cell MLE grid is missing or does not match the ",
    "selected failed MCMC grid.",
    call. = FALSE
  )
}

diagnostic_mle_grid_tro <- diagnostic_mle_grid$biomass |>
  mutate(
    h = factor(
      formatC(.data$h, format = "f", digits = 2L),
      levels = formatC(
        sort(unique(diagnostic_mle_grid$biomass$h)),
        format = "f",
        digits = 2L
      )
    ),
    psi = factor(
      paste0(
        "\u03c8 = ",
        formatC(.data$psi, format = "f", digits = 2L)
      ),
      levels = paste0(
        "\u03c8 = ",
        formatC(
          sort(unique(diagnostic_mle_grid$biomass$psi)),
          format = "f",
          digits = 2L
        )
      )
    ),
    M10 = factor(
      formatC(.data$fixed_m10, format = "f", digits = 3L),
      levels = formatC(
        sort(unique(.data$fixed_m10)),
        format = "f",
        digits = 3L
      )
    )
  ) |>
  distinct(.data$h, .data$psi, .data$M10, .data$Year, .keep_all = TRUE)
} else {
  diagnostic_mle_grid_tro <- tibble()
}
```

```{r}
#| label: fig-diagnostic-mle-grid-tro
#| echo: false
#| fig-cap: "Archived diagnostic relative total reproductive output for the 108-row MLE-grid layout. This figure is displayed only when a source MCMC grid cell fails; it is not a posterior uncertainty result or a substitute for MCMC-grid acceptance."
#| fig-width: 12
#| fig-height: 6.5
if (length(failed_grid_cells)) {
  print(
    ggplot(
      diagnostic_mle_grid_tro,
      aes(
        x = .data$Year,
        y = .data$relative_spawning_biomass,
        colour = .data$h,
        linetype = .data$M10,
        group = interaction(.data$h, .data$M10)
      )
    ) +
      geom_line(linewidth = 0.7, alpha = 0.9) +
      facet_wrap(vars(.data$psi), nrow = 1L) +
      scale_colour_brewer(palette = "Dark2") +
      scale_x_continuous(breaks = pretty_breaks(n = 6L)) +
      scale_y_zero() +
      labs(
        x = "Year",
        y = "Relative TRO",
        colour = "h",
        linetype = "M10"
      ) +
      theme(legend.position = "bottom")
  )
}
```

```{r}
#| label: failed-grid-diagnostic-checkpoint
#| echo: false
#| results: asis
if (length(failed_grid_cells)) {
  if (!render_failed_grid_diagnostics) {
    stop(
      "Grid cells ", paste(failed_grid_cells, collapse = ", "),
      " fail the production thresholds. Set ",
      "ESC31_RENDER_FAILED_GRID_DIAGNOSTICS=true to render a diagnostic ",
      "checkpoint without combining cells or running dependent analyses.",
      call. = FALSE
    )
  }
  cat(
    "::: {.callout-warning}\n",
    "## Diagnostic checkpoint only\n\n",
    "Grid cells **", paste(failed_grid_cells, collapse = ", "),
    "** fail at least one production threshold. This page therefore ends ",
    "after the embedded payload-bound diagnostics. No failed MCMC cell has ",
    "been accepted or combined, and no balanced posterior or projection ",
    "result is shown.\n",
    ":::\n",
    sep = ""
  )
  knitr::knit_exit()
}
```

### MCMC-grid TRO trajectories

Figures @fig-grid-relative-spawning-biomass and
@fig-grid-relative-spawning-biomass-overlay 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.

```{r}
#| label: grid-biomass-all-draw-summary
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)
  )
```

```{r}
#| label: fig-grid-relative-spawning-biomass
#| fig-cap: "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."
#| fig-width: 12
#| fig-height: 9
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()
```

```{r}
#| label: fig-grid-relative-spawning-biomass-overlay
#| fig-cap: "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."
#| fig-width: 11
#| fig-height: 6
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")
```

### 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.

```{r}
#| label: combine-grid-mcmc
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)
}
}
```

```{r}
#| label: tbl-grid-draws
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."
    )
  )
```

```{r}
#| label: combined-grid-mortality-posterior
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
}
```

```{r}
#| label: fig-mcmc-grid-lev-resample
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 10
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
  )
)
```

```{r}
#| label: combined-grid-state
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"
  )
```

```{r}
#| label: fig-grid-combined-relative-spawning-biomass-by-cell
#| fig-cap: "Sample of relative total reproductive output trajectories from the balanced 2,000-draw projection posterior, colored by source grid cell."
#| fig-width: 11
#| fig-height: 6
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()
```

```{r}
#| label: fig-grid-combined-relative-spawning-biomass
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 5.5
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()
```

```{r}
#| label: fig-grid-combined-recruitment
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 5.5
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))
```

### 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.

```{r}
#| label: prepare-grid-combined-kobe
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")
)
```

```{r}
#| label: fig-grid-combined-kobe
#| fig-cap: "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."
#| fig-width: 11
#| fig-height: 7.5
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
  )
```

The 2025 posterior is split between the lower-left and lower-right Kobe
quadrants: `r percent(kobe_quadrants_2025$yellow, accuracy = 0.1)` of draws
have $TRO < TRO_{MSY}$, while `r percent(kobe_quadrants_2025$green, accuracy = 0.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.

```{r}
#| label: prepare-grid-fao-kobe
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 = 0.88,
    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 = 7.6,
  units = "in",
  dpi = 320,
  bg = "white"
)
```

```{r}
#| label: fig-grid-fao-kobe
#| fig-cap: "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."
#| fig-width: 8.4
#| fig-height: 7.6
fao_kobe_plot
```

### 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.

```{r}
#| label: prepare-sensitivity-grid-comparison
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 `r sensitivity_terminal_spread$Year`, the sensitivity posterior medians span
`r format_decimal(sensitivity_terminal_spread$minimum, 3)` to
`r format_decimal(sensitivity_terminal_spread$maximum, 3)` relative TRO, versus
`r format_decimal(grid_terminal_spread$minimum, 3)` to
`r format_decimal(grid_terminal_spread$maximum, 3)` across the nine MCMC grid
cells. The sensitivity span is therefore
`r format_decimal(100 * (sensitivity_grid_terminal_spread_ratio - 1), 1)`%
wider. `r sensitivity_terminal_spread$minimum_model` is below the lowest grid
median, and `r sensitivity_terminal_spread$maximum_model` is above the highest
grid median. The $h$/$\psi$ grid is consequently substantial, but it does not
envelop the full structural sensitivity range.

```{r}
#| label: fig-sensitivity-grid-mcmc-tro-comparison
#| fig-cap: "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."
#| fig-width: 11
#| fig-height: 6
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")
```

```{r}
#| label: tbl-sensitivity-grid-terminal-tro
#| tbl-cap: "Terminal-year relative TRO posterior summaries for the accepted sensitivity fits and all nine MCMC grid cells."
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()
```

## 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.

```{r}
#| label: m10-mle-grid-functions
# 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.

```

```{r}
#| label: run-m10-mle-grid
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 `r nrow(m10_mle_summary)` direct-M cells are distinct and pass. Objectives
range from `r format_decimal(min(m10_mle_summary$objective), 3)` to
`r format_decimal(max(m10_mle_summary$objective), 3)`. The largest maximum
gradient is `r format_decimal(max(m10_mle_summary$max_gradient), 6)` (cell
`r m10_mle_summary$Cell[which.max(m10_mle_summary$max_gradient)]`), maximum
raw harvest is
`r format_decimal(max(m10_mle_summary$state_max_raw_harvest), 3)`, and minimum
populated abundance is
`r format_decimal(min(m10_mle_summary$state_min_number), 3)`. Estimated M4
ranges from `r format_decimal(min(m10_mle_summary$estimated_m4), 3)` to
`r format_decimal(max(m10_mle_summary$estimated_m4), 3)`, and estimated M30
ranges from `r format_decimal(min(m10_mle_summary$estimated_m30), 3)` to
`r format_decimal(max(m10_mle_summary$estimated_m30), 3)`. The validated
projection-ready object contains
`r format_decimal(length(mle_grid_sample$grid_cells), 0)` draws and samples
`r sum(mle_grid_sample$grid_freq$Freq > 0L)` of the 108 cells with nonzero
frequency.

```{r}
#| label: tbl-m10-posterior-grid
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.")
```

```{r}
#| label: tbl-m10-mle-grid
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.")
```

```{r}
#| label: fig-m10-posterior-grid
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 5
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()
```

```{r}
#| label: tbl-mle-grid-resampling
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."))
```

```{r}
#| label: fig-mle-grid-lev-resample
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 10
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
  )
)
```

```{r}
#| label: mle-grid-biomass-summary
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)
```

```{r}
#| label: mle-grid-mortality-summary
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"
  )
```

```{r}
#| label: fig-grid-mle-mcmc-mortality-at-age
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 5.5
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()
```

```{r}
#| label: fig-grid-mle-mcmc-relative-spawning-biomass
#| fig-cap: "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."
#| fig-width: 10
#| fig-height: 5.5
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()
```

In `r mle_grid_comparison_terminal_year`, the MCMC-grid relative-TRO median is
`r format_decimal(mcmc_grid_terminal_biomass$median, 3)` (95% interval
`r format_decimal(mcmc_grid_terminal_biomass$lower, 3)`--`r format_decimal(mcmc_grid_terminal_biomass$upper, 3)`),
whereas the full-objective-weighted MLE-grid median is
`r format_decimal(mle_grid_terminal_biomass$median, 3)` (95% interval
`r format_decimal(mle_grid_terminal_biomass$lower, 3)`--`r format_decimal(mle_grid_terminal_biomass$upper, 3)`).
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.

```{r}
#| label: prepare-previous-current-mle-grid-comparison
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)
    )
)
```

```{r}
#| label: tbl-previous-current-mle-grid-coordinates
#| tbl-cap: "All-cell design comparison for the previous CCSBT base-assessment grid and the current ESC31 MLE grid."
previous_current_mle_grid_coordinates |>
  mutate(
    Cells = format_decimal(Cells, digits = 0),
    `Final population state` =
      format_decimal(`Final population state`, digits = 0)
  ) |>
  kable()
```

```{r}
#| label: fig-previous-current-mle-grid-all-cells
#| fig-cap: "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."
#| fig-width: 11
#| fig-height: 6
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")
```

## Appendix: MCMC Diagnostics

### Sampler Parameters

```{r}
#| label: sampler-params-helper
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)
}
```

::: {.panel-tabset}

#### Cell 1

```{r}
#| label: fig-grid-sampler-cell-1
#| fig-cap: "Sampler-parameter diagnostics for grid cell 1."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(1)
```

#### Cell 2

```{r}
#| label: fig-grid-sampler-cell-2
#| fig-cap: "Sampler-parameter diagnostics for grid cell 2."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(2)
```

#### Cell 3

```{r}
#| label: fig-grid-sampler-cell-3
#| fig-cap: "Sampler-parameter diagnostics for grid cell 3."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(3)
```

#### Cell 4

```{r}
#| label: fig-grid-sampler-cell-4
#| fig-cap: "Sampler-parameter diagnostics for grid cell 4."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(4)
```

#### Cell 5

```{r}
#| label: fig-grid-sampler-cell-5
#| fig-cap: "Sampler-parameter diagnostics for grid cell 5."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(5)
```

#### Cell 6

```{r}
#| label: fig-grid-sampler-cell-6
#| fig-cap: "Sampler-parameter diagnostics for grid cell 6."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(6)
```

#### Cell 7

```{r}
#| label: fig-grid-sampler-cell-7
#| fig-cap: "Sampler-parameter diagnostics for grid cell 7."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(7)
```

#### Cell 8

```{r}
#| label: fig-grid-sampler-cell-8
#| fig-cap: "Sampler-parameter diagnostics for grid cell 8."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(8)
```

#### Cell 9

```{r}
#| label: fig-grid-sampler-cell-9
#| fig-cap: "Sampler-parameter diagnostics for grid cell 9."
#| fig-width: 10
#| fig-height: 8
plot_sampler_cell(9)
```

:::

### MCMC Pair Plots

```{r}
#| label: pairs-helper
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"
  )
}
```

::: {.panel-tabset}

#### Cell 1

```{r}
#| label: fig-grid-pairs-cell-1
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 1, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(1)
```

#### Cell 2

```{r}
#| label: fig-grid-pairs-cell-2
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 2, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(2)
```

#### Cell 3

```{r}
#| label: fig-grid-pairs-cell-3
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 3, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(3)
```

#### Cell 4

```{r}
#| label: fig-grid-pairs-cell-4
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 4, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(4)
```

#### Cell 5

```{r}
#| label: fig-grid-pairs-cell-5
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 5, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(5)
```

#### Cell 6

```{r}
#| label: fig-grid-pairs-cell-6
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 6, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(6)
```

#### Cell 7

```{r}
#| label: fig-grid-pairs-cell-7
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 7, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(7)
```

#### Cell 8

```{r}
#| label: fig-grid-pairs-cell-8
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 8, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(8)
```

#### Cell 9

```{r}
#| label: fig-grid-pairs-cell-9
#| fig-cap: "pairs_rtmb diagnostic plot for grid cell 9, showing the five slowest-mixing parameters."
#| fig-width: 12
#| fig-height: 10
plot_pairs_cell(9)
```

:::