args <- commandArgs(trailingOnly = FALSE)
file_arg <- grep("^--file=", args, value = TRUE)
analysis_dir <- if (length(file_arg)) dirname(normalizePath(sub("^--file=", "", file_arg[[1]]), mustWork = FALSE)) else file.path(getwd(), "analysis")
source(file.path(analysis_dir, "00_config.R"))
source(file.path(analysis_dir, "lib_timss.R"))

replace_file <- function(staged, destination) {
  dir.create(dirname(destination), recursive = TRUE, showWarnings = FALSE)
  backup <- paste0(destination, ".backup-", Sys.getpid())
  if (file.exists(backup)) unlink(backup)
  had_destination <- file.exists(destination)
  if (had_destination && !file.rename(destination, backup)) {
    stop("Could not stage existing file for replacement: ", destination)
  }
  installed <- file.rename(staged, destination)
  if (!installed) {
    if (had_destination) file.rename(backup, destination)
    stop("Could not install validated file: ", destination)
  }
  if (had_destination) unlink(backup)
  invisible(destination)
}

specification_systems <- function(specification, systems) {
  set <- standardization_specifications[[specification]]$systems
  if (identical(set, "all")) return(systems)
  if (identical(set, "without_korea")) return(setdiff(systems, "Korea, Rep. of"))
  stop("Unknown standardization system set: ", set)
}

make_design <- function(math_t0, science_t0, threshold, surface = "quadratic_interaction") {
  math_deficit <- (threshold - math_t0) / 100
  science_deficit <- (threshold - science_t0) / 100
  if (surface == "quadratic_interaction") {
    return(cbind(
      intercept = 1,
      math_deficit = math_deficit,
      science_deficit = science_deficit,
      math_deficit_sq = math_deficit^2,
      science_deficit_sq = science_deficit^2,
      deficit_interaction = math_deficit * science_deficit
    ))
  }
  if (surface == "additive_linear") {
    return(cbind(
      intercept = 1,
      math_deficit = math_deficit,
      science_deficit = science_deficit
    ))
  }
  if (surface == "linear_interaction") {
    return(cbind(
      intercept = 1,
      math_deficit = math_deficit,
      science_deficit = science_deficit,
      deficit_interaction = math_deficit * science_deficit
    ))
  }
  stop("Unknown recovery surface: ", surface)
}

fit_recovery_model <- function(design, outcome, weight, collect_diagnostics = FALSE) {
  valid <- is.finite(outcome) & is.finite(weight) & weight > 0 & rowSums(!is.finite(design)) == 0
  if (sum(valid) <= ncol(design)) stop("Too few positive-weight records for standardization model")
  model_weight <- weight[valid]
  # Constant rescaling leaves the weighted score equations unchanged and avoids
  # binomial initialization collapsing to machine-level probabilities.
  model_weight <- model_weight / mean(model_weight)
  warning_messages <- character(0)
  fit <- withCallingHandlers(
    stats::glm.fit(
      x = design[valid, , drop = FALSE],
      y = outcome[valid],
      weights = model_weight,
      family = stats::quasibinomial(),
      intercept = FALSE,
      control = stats::glm.control(maxit = 200, epsilon = 1e-9)
    ),
    warning = function(warning) {
      warning_messages <<- c(warning_messages, conditionMessage(warning))
      invokeRestart("muffleWarning")
    }
  )
  coefficients <- fit$coefficients
  if (fit$rank != ncol(design)) stop("Rank-deficient standardization model")
  if (!isTRUE(fit$converged)) stop("Standardization model did not converge")
  if (any(!is.finite(coefficients))) stop("Non-finite standardization coefficient")

  result <- list(
    coefficients = coefficients,
    converged = isTRUE(fit$converged),
    iterations = fit$iter,
    rank = fit$rank,
    warning_count = length(warning_messages)
  )
  if (collect_diagnostics) {
    weighted_design <- design[valid, , drop = FALSE] * sqrt(model_weight)
    result$valid <- valid
    result$condition_number <- kappa(weighted_design, exact = TRUE)
    result$positive_weight_n <- sum(valid)
    result$event_n <- sum(outcome[valid] == 1)
  }
  result
}

predict_recovery <- function(model, design) {
  linear_predictor <- as.numeric(design %*% model$coefficients)
  clipped <- pmax(-config$standardization_eta_limit, pmin(config$standardization_eta_limit, linear_predictor))
  list(
    probability = stats::plogis(clipped),
    linear_predictor = linear_predictor,
    clipped = linear_predictor != clipped
  )
}

system_bounds <- function(math_t0, science_t0, weight, probabilities) {
  c(
    math_lower = weighted_quantile(math_t0, weight, probabilities[1]),
    math_upper = weighted_quantile(math_t0, weight, probabilities[2]),
    science_lower = weighted_quantile(science_t0, weight, probabilities[1]),
    science_upper = weighted_quantile(science_t0, weight, probabilities[2])
  )
}

common_bounds <- function(bounds_matrix) {
  bounds <- c(
    math_lower = max(bounds_matrix[, "math_lower"]),
    math_upper = min(bounds_matrix[, "math_upper"]),
    science_lower = max(bounds_matrix[, "science_lower"]),
    science_upper = min(bounds_matrix[, "science_upper"])
  )
  if (any(!is.finite(bounds)) || bounds["math_lower"] >= bounds["math_upper"] ||
      bounds["science_lower"] >= bounds["science_upper"]) {
    stop("No finite rectangular common support")
  }
  bounds
}

target_from_eligibility <- function(pooled, analysis_weight, eligible, systems, target) {
  reference <- numeric(nrow(pooled))
  if (target == "equal_system") {
    for (system in systems) {
      included <- eligible & pooled$system == system & analysis_weight > 0
      included_total <- sum(analysis_weight[included])
      if (!is.finite(included_total) || included_total <= 0) stop("Empty common support for ", system)
      reference[included] <- analysis_weight[included] / included_total / length(systems)
    }
  } else if (target == "pooled_population") {
    included <- eligible & pooled$system %in% systems & analysis_weight > 0
    included_total <- sum(analysis_weight[included])
    if (!is.finite(included_total) || included_total <= 0) stop("Empty pooled-population target")
    reference[included] <- analysis_weight[included] / included_total
  } else {
    stop("Unknown standardization target: ", target)
  }
  if (abs(sum(reference) - 1) > 1e-10) stop("Reference weights do not sum to one")
  reference
}

rectangle_support <- function(pooled, analysis_weight, systems, probabilities, target) {
  bounds_matrix <- do.call(rbind, lapply(systems, function(system) {
    index <- pooled$system == system
    system_bounds(pooled$math_t0[index], pooled$science_t0[index], analysis_weight[index], probabilities)
  }))
  rownames(bounds_matrix) <- systems
  bounds <- common_bounds(bounds_matrix)
  eligible <- pooled$system %in% systems &
    pooled$math_t0 >= bounds["math_lower"] & pooled$math_t0 <= bounds["math_upper"] &
    pooled$science_t0 >= bounds["science_lower"] & pooled$science_t0 <= bounds["science_upper"]
  reference <- target_from_eligibility(pooled, analysis_weight, eligible, systems, target)
  list(
    systems = systems,
    probabilities = probabilities,
    bounds_matrix = bounds_matrix,
    bounds = bounds,
    eligible = eligible,
    target_weight = reference,
    target = target,
    support = "rectangle_01_99",
    target_mass_retained = 1,
    grid_cells_retained = NA_integer_,
    grid_cells_total = NA_integer_,
    minimum_common_cell_n = NA_integer_,
    minimum_common_cell_ess = NA_real_,
    math_breaks = rep(NA_real_, config$strict_cell_grid + 1L),
    science_breaks = rep(NA_real_, config$strict_cell_grid + 1L)
  )
}

strict_cell_support <- function(pooled, analysis_weight, primary_support) {
  grid <- config$strict_cell_grid
  positive_target <- primary_support$target_weight > 0
  probabilities <- seq(0, 1, length.out = grid + 1L)
  math_breaks <- weighted_quantile(
    pooled$math_t0[positive_target], primary_support$target_weight[positive_target], probabilities
  )
  science_breaks <- weighted_quantile(
    pooled$science_t0[positive_target], primary_support$target_weight[positive_target], probabilities
  )
  if (any(!is.finite(c(math_breaks, science_breaks))) ||
      anyDuplicated(math_breaks) || anyDuplicated(science_breaks)) {
    stop("Degenerate strict common-cell grid")
  }
  math_breaks[c(1L, length(math_breaks))] <-
    math_breaks[c(1L, length(math_breaks))] + c(-1e-8, 1e-8)
  science_breaks[c(1L, length(science_breaks))] <-
    science_breaks[c(1L, length(science_breaks))] + c(-1e-8, 1e-8)
  math_bin <- cut(pooled$math_t0, breaks = math_breaks, labels = FALSE, include.lowest = TRUE)
  science_bin <- cut(pooled$science_t0, breaks = science_breaks, labels = FALSE, include.lowest = TRUE)
  cell <- (math_bin - 1L) * grid + science_bin

  cell_rows <- list()
  row <- 1L
  supported_cells <- integer(0)
  for (candidate in seq_len(grid^2L)) {
    candidate_supported <- TRUE
    for (system in primary_support$systems) {
      index <- primary_support$eligible & pooled$system == system & cell == candidate & analysis_weight > 0
      unweighted_n <- sum(index)
      effective_n <- weighted_effective_n(analysis_weight[index])
      cell_rows[[row]] <- data.frame(
        cell = candidate,
        system = system,
        unweighted_n = unweighted_n,
        effective_n = effective_n,
        stringsAsFactors = FALSE
      )
      row <- row + 1L
      if (unweighted_n < config$strict_cell_min_n || effective_n < config$strict_cell_min_ess) {
        candidate_supported <- FALSE
      }
    }
    if (candidate_supported) supported_cells <- c(supported_cells, candidate)
  }
  if (!length(supported_cells)) stop("Strict common-cell rule retained no cells")
  cell_diagnostics <- do.call(rbind, cell_rows)
  eligible <- primary_support$eligible & cell %in% supported_cells
  reference <- target_from_eligibility(
    pooled, analysis_weight, eligible, primary_support$systems, "equal_system"
  )
  supported_diagnostics <- cell_diagnostics[cell_diagnostics$cell %in% supported_cells, ]
  list(
    systems = primary_support$systems,
    probabilities = primary_support$probabilities,
    bounds_matrix = primary_support$bounds_matrix,
    bounds = primary_support$bounds,
    eligible = eligible,
    target_weight = reference,
    target = "equal_system",
    support = "strict_common_cells",
    target_mass_retained = sum(primary_support$target_weight[eligible]),
    grid_cells_retained = length(supported_cells),
    grid_cells_total = grid^2L,
    minimum_common_cell_n = min(supported_diagnostics$unweighted_n),
    minimum_common_cell_ess = min(supported_diagnostics$effective_n),
    math_breaks = math_breaks,
    science_breaks = science_breaks,
    cell = cell,
    supported_cells = supported_cells,
    cell_diagnostics = cell_diagnostics
  )
}

build_supports <- function(pooled, analysis_weight, systems, probabilities) {
  if (length(probabilities) != 2L || any(!is.finite(probabilities)) ||
      probabilities[1] < 0 || probabilities[2] > 1 || probabilities[1] >= probabilities[2]) {
    stop("Common-support probabilities must be two ordered finite values in [0, 1]")
  }
  systems_without_korea <- setdiff(systems, "Korea, Rep. of")
  primary <- rectangle_support(
    pooled, analysis_weight, systems, probabilities, "equal_system"
  )
  population <- rectangle_support(
    pooled, analysis_weight, systems, probabilities, "pooled_population"
  )
  leave_korea_out <- rectangle_support(
    pooled, analysis_weight, systems_without_korea, probabilities, "equal_system"
  )
  strict <- strict_cell_support(pooled, analysis_weight, primary)
  list(
    primary = primary,
    population_target = population,
    leave_korea_out = leave_korea_out,
    strict_common_cells = strict,
    additive_linear = primary,
    linear_interaction = primary
  )
}

prepare_draw <- function(data, draw, threshold, systems) {
  map <- pv_variable_map("g45", draw)
  math_t0 <- data[[map$math_t0]]
  science_t0 <- data[[map$science_t0]]
  state_t0 <- achievement_state(math_t0, science_t0, threshold)
  state_t1 <- achievement_state(data[[map$math_t1]], data[[map$science_t1]], threshold)
  risk <- state_t0 == "BB"
  pooled <- data.frame(
    system = data$system[risk],
    math_t0 = math_t0[risk],
    science_t0 = science_t0[risk],
    recovered = as.numeric(state_t1[risk] == "AA"),
    weight = data$weight[risk],
    jkzone = data$jkzone[risk],
    jkrep = data$jkrep[risk],
    stringsAsFactors = FALSE
  )
  pooled <- pooled[pooled$system %in% systems, ]
  pooled$system <- factor(pooled$system, levels = systems)
  if (any(!is.finite(pooled$recovered)) || any(!is.finite(pooled$weight)) || any(pooled$weight <= 0)) {
    stop("Invalid baseline-BB analysis records")
  }
  pooled
}

raw_parameter <- function(system) paste0("raw_full::", system)
component_parameter <- function(specification, component, system) {
  paste(specification, component, system, sep = "::")
}

fit_full_surfaces <- function(pooled, designs, systems, draw, common_target_weight) {
  if (length(common_target_weight) != nrow(pooled) || any(!is.finite(common_target_weight)) ||
      any(common_target_weight < 0) || abs(sum(common_target_weight) - 1) > 1e-10) {
    stop("Invalid primary common-target weights for model diagnostics")
  }
  predictions <- list()
  models <- list()
  diagnostics <- list()
  coefficients <- list()
  diagnostic_row <- coefficient_row <- 1L
  for (surface in names(designs)) {
    prediction_matrix <- matrix(
      NA_real_, nrow = nrow(pooled), ncol = length(systems),
      dimnames = list(NULL, systems)
    )
    surface_models <- setNames(vector("list", length(systems)), systems)
    for (system in systems) {
      index <- pooled$system == system
      model <- fit_recovery_model(
        designs[[surface]][index, , drop = FALSE], pooled$recovered[index], pooled$weight[index],
        collect_diagnostics = TRUE
      )
      prediction <- predict_recovery(model, designs[[surface]])
      prediction_matrix[, system] <- prediction$probability
      surface_models[[system]] <- model
      own_prediction <- prediction$probability[index]
      pooled_clipped <- prediction$clipped
      own_clipped <- pooled_clipped[index]
      common_target_clipped <- pooled_clipped & common_target_weight > 0
      diagnostics[[diagnostic_row]] <- data.frame(
        draw = draw,
        surface = surface,
        system = system,
        risk_n = sum(index),
        event_n = sum(pooled$recovered[index] == 1),
        design_columns = ncol(designs[[surface]]),
        model_rank = model$rank,
        condition_number = model$condition_number,
        model_iterations = model$iterations,
        converged = model$converged,
        warning_count = model$warning_count,
        maximum_absolute_coefficient = max(abs(model$coefficients)),
        weighted_observed_rate = weighted_probability(pooled$recovered[index] == 1, pooled$weight[index]),
        weighted_fitted_rate = weighted_mean(own_prediction, pooled$weight[index]),
        weighted_calibration_difference = weighted_mean(own_prediction, pooled$weight[index]) -
          weighted_probability(pooled$recovered[index] == 1, pooled$weight[index]),
        weighted_brier_score = weighted_mean((pooled$recovered[index] - own_prediction)^2, pooled$weight[index]),
        linear_predictor_minimum = min(prediction$linear_predictor[index]),
        linear_predictor_maximum = max(prediction$linear_predictor[index]),
        predicted_probability_minimum = min(own_prediction),
        predicted_probability_maximum = max(own_prediction),
        pooled_clipped_prediction_n = sum(pooled_clipped),
        pooled_clipped_prediction_weighted_share = weighted_probability(pooled_clipped, pooled$weight),
        own_system_clipped_prediction_n = sum(own_clipped),
        own_system_clipped_prediction_weighted_share = weighted_probability(own_clipped, pooled$weight[index]),
        common_target_clipped_prediction_n = sum(common_target_clipped),
        common_target_clipped_prediction_weighted_share = as.numeric(crossprod(common_target_weight, pooled_clipped)),
        stringsAsFactors = FALSE
      )
      diagnostic_row <- diagnostic_row + 1L
      for (term in names(model$coefficients)) {
        coefficients[[coefficient_row]] <- data.frame(
          draw = draw,
          surface = surface,
          system = system,
          term = term,
          estimate = unname(model$coefficients[term]),
          stringsAsFactors = FALSE
        )
        coefficient_row <- coefficient_row + 1L
      }
    }
    predictions[[surface]] <- prediction_matrix
    models[[surface]] <- surface_models
  }
  list(
    predictions = predictions,
    models = models,
    diagnostics = do.call(rbind, diagnostics),
    coefficients = do.call(rbind, coefficients)
  )
}

evaluate_parameters <- function(pooled, analysis_weight, supports, predictions, systems) {
  values <- setNames(numeric(0), character(0))
  for (system in systems) {
    index <- pooled$system == system
    values[raw_parameter(system)] <- weighted_probability(
      pooled$recovered[index] == 1, analysis_weight[index]
    )
  }
  for (specification in names(standardization_specifications)) {
    definition <- standardization_specifications[[specification]]
    support <- supports[[specification]]
    specification_system_set <- specification_systems(specification, systems)
    prediction_matrix <- predictions[[definition$surface]]
    for (system in specification_system_set) {
      index <- support$eligible & pooled$system == system & analysis_weight > 0
      if (!any(index)) stop("Empty supported system distribution: ", specification, " / ", system)
      raw_support <- weighted_probability(pooled$recovered[index] == 1, analysis_weight[index])
      model_own_support <- weighted_mean(prediction_matrix[index, system], analysis_weight[index])
      standardized_common <- as.numeric(crossprod(support$target_weight, prediction_matrix[, system]))
      values[component_parameter(specification, "raw_support", system)] <- raw_support
      values[component_parameter(specification, "model_own_support", system)] <- model_own_support
      values[component_parameter(specification, "standardized_common", system)] <- standardized_common
    }
  }
  values
}

support_diagnostic_rows <- function(pooled, analysis_weight, supports, systems, draw) {
  rows <- list()
  row <- 1L
  for (specification in names(standardization_specifications)) {
    support <- supports[[specification]]
    specification_system_set <- specification_systems(specification, systems)
    for (system in specification_system_set) {
      in_system <- pooled$system == system & analysis_weight > 0
      included <- support$eligible & in_system
      denominator <- sum(analysis_weight[in_system])
      included_total <- sum(analysis_weight[included])
      rows[[row]] <- data.frame(
        draw = draw,
        specification = specification,
        system = system,
        surface = standardization_specifications[[specification]]$surface,
        target = standardization_specifications[[specification]]$target,
        support = standardization_specifications[[specification]]$support,
        risk_n = sum(in_system),
        supported_n = sum(included),
        support_retained_weight = included_total / denominator,
        support_effective_n = weighted_effective_n(analysis_weight[included]),
        target_system_mass = sum(support$target_weight[pooled$system == system]),
        target_weight_sum = sum(support$target_weight),
        target_mass_retained_before_renormalization = support$target_mass_retained,
        math_lower = support$bounds["math_lower"],
        math_upper = support$bounds["math_upper"],
        science_lower = support$bounds["science_lower"],
        science_upper = support$bounds["science_upper"],
        grid_cells_retained = support$grid_cells_retained,
        grid_cells_total = support$grid_cells_total,
        minimum_common_cell_n = support$minimum_common_cell_n,
        minimum_common_cell_ess = support$minimum_common_cell_ess,
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
  do.call(rbind, rows)
}

standardize_draw <- function(data, draw, threshold, systems, probabilities) {
  pooled <- prepare_draw(data, draw, threshold, systems)
  indices <- lapply(systems, function(system) which(pooled$system == system))
  names(indices) <- systems
  surfaces <- unique(vapply(standardization_specifications, `[[`, character(1), "surface"))
  specifications <- names(standardization_specifications)
  designs <- setNames(lapply(surfaces, function(surface) {
    make_design(pooled$math_t0, pooled$science_t0, threshold, surface)
  }), surfaces)

  supports_full <- build_supports(pooled, pooled$weight, systems, probabilities)
  full_surfaces <- fit_full_surfaces(
    pooled, designs, systems, draw, supports_full$primary$target_weight
  )
  point <- evaluate_parameters(
    pooled, pooled$weight, supports_full, full_surfaces$predictions, systems
  )
  support_diagnostics <- support_diagnostic_rows(
    pooled, pooled$weight, supports_full, systems, draw
  )

  zone_count <- vapply(indices, function(index) {
    length(unique(pooled$jkzone[index][is.finite(pooled$jkzone[index])]))
  }, integer(1))
  if (any(zone_count <= 0)) stop("A system has no finite JRR zones in the baseline-BB risk set")
  replicate_total <- 2L * sum(zone_count)
  expected_affected_replicates <- setNames(vapply(specifications, function(specification) {
    2L * sum(zone_count[specification_systems(specification, systems)])
  }, integer(1)), specifications)
  deviations <- matrix(
    NA_real_, nrow = replicate_total, ncol = length(point),
    dimnames = list(NULL, names(point))
  )

  replicate_attempted <- 0L
  replicate_completed <- 0L
  support_rebuild_attempted <- 0L
  support_rebuild_completed <- 0L
  affected_replicates_attempted <- setNames(integer(length(specifications)), specifications)
  affected_replicates_completed <- setNames(integer(length(specifications)), specifications)
  surface_refits_attempted <- setNames(integer(length(specifications)), specifications)
  surface_refits_completed <- setNames(integer(length(specifications)), specifications)
  refits_with_warnings <- setNames(integer(length(specifications)), specifications)
  warning_count_across_refits <- setNames(integer(length(specifications)), specifications)
  refits_with_any_pooled_clipping <- setNames(integer(length(specifications)), specifications)
  pooled_clipped_rows <- setNames(integer(length(specifications)), specifications)
  maximum_pooled_clipped_share <- setNames(numeric(length(specifications)), specifications)
  refits_with_any_own_clipping <- setNames(integer(length(specifications)), specifications)
  own_clipped_rows <- setNames(integer(length(specifications)), specifications)
  maximum_own_clipped_share <- setNames(numeric(length(specifications)), specifications)
  refits_with_any_target_clipping <- setNames(integer(length(specifications)), specifications)
  target_clipped_rows <- setNames(integer(length(specifications)), specifications)
  maximum_target_clipped_share <- setNames(numeric(length(specifications)), specifications)
  replicate_row <- 1L

  for (perturbed_system in systems) {
    index <- indices[[perturbed_system]]
    zones <- sort(unique(pooled$jkzone[index][is.finite(pooled$jkzone[index])]))
    for (zone in zones) {
      for (included_member in c(0, 1)) {
        affected_specifications <- specifications[vapply(
          specifications,
          function(specification) perturbed_system %in% specification_systems(specification, systems),
          logical(1)
        )]
        replicate_attempted <- replicate_attempted + 1L
        affected_replicates_attempted[affected_specifications] <-
          affected_replicates_attempted[affected_specifications] + 1L
        replicate_weights <- pooled$weight
        replicate_weights[index] <- replicate_weight(
          pooled$weight[index], pooled$jkzone[index], pooled$jkrep[index], zone, included_member
        )
        support_rebuild_attempted <- support_rebuild_attempted + 1L
        replicate_supports <- tryCatch(
          build_supports(pooled, replicate_weights, systems, probabilities),
          error = function(error) stop(
            "Support rebuild failed in draw ", draw, ", ", perturbed_system,
            ", zone ", zone, ", member ", included_member, ": ", conditionMessage(error)
          )
        )
        support_rebuild_completed <- support_rebuild_completed + 1L
        replicate_predictions <- full_surfaces$predictions
        for (surface in surfaces) {
          relevant_specifications <- affected_specifications[vapply(
            affected_specifications,
            function(specification) identical(standardization_specifications[[specification]]$surface, surface),
            logical(1)
          )]
          surface_refits_attempted[relevant_specifications] <-
            surface_refits_attempted[relevant_specifications] + 1L
          model <- tryCatch(
            fit_recovery_model(
              designs[[surface]][index, , drop = FALSE],
              pooled$recovered[index], replicate_weights[index]
            ),
            error = function(error) stop(
              "Model refit failed in draw ", draw, ", surface ", surface, ", ",
              perturbed_system, ", zone ", zone, ", member ", included_member,
              ": ", conditionMessage(error)
            )
          )
          surface_refits_completed[relevant_specifications] <-
            surface_refits_completed[relevant_specifications] + 1L
          refits_with_warnings[relevant_specifications] <-
            refits_with_warnings[relevant_specifications] + as.integer(model$warning_count > 0L)
          warning_count_across_refits[relevant_specifications] <-
            warning_count_across_refits[relevant_specifications] + model$warning_count

          prediction <- predict_recovery(model, designs[[surface]])
          replicate_predictions[[surface]][, perturbed_system] <- prediction$probability
          pooled_indicator <- prediction$clipped
          own_indicator <- pooled_indicator[index]
          pooled_share <- weighted_probability(pooled_indicator, replicate_weights)
          own_share <- weighted_probability(own_indicator, replicate_weights[index])
          for (specification in relevant_specifications) {
            target_weight <- replicate_supports[[specification]]$target_weight
            target_indicator <- pooled_indicator & target_weight > 0
            target_share <- as.numeric(crossprod(target_weight, pooled_indicator))
            refits_with_any_pooled_clipping[specification] <-
              refits_with_any_pooled_clipping[specification] + as.integer(any(pooled_indicator))
            pooled_clipped_rows[specification] <-
              pooled_clipped_rows[specification] + sum(pooled_indicator)
            maximum_pooled_clipped_share[specification] <- max(
              maximum_pooled_clipped_share[specification], pooled_share
            )
            refits_with_any_own_clipping[specification] <-
              refits_with_any_own_clipping[specification] + as.integer(any(own_indicator))
            own_clipped_rows[specification] <- own_clipped_rows[specification] + sum(own_indicator)
            maximum_own_clipped_share[specification] <- max(
              maximum_own_clipped_share[specification], own_share
            )
            refits_with_any_target_clipping[specification] <-
              refits_with_any_target_clipping[specification] + as.integer(any(target_indicator))
            target_clipped_rows[specification] <-
              target_clipped_rows[specification] + sum(target_indicator)
            maximum_target_clipped_share[specification] <- max(
              maximum_target_clipped_share[specification], target_share
            )
          }
        }
        theta_replicate <- evaluate_parameters(
          pooled, replicate_weights, replicate_supports, replicate_predictions, systems
        )
        if (!identical(names(theta_replicate), names(point))) stop("Replicate parameter order changed")
        deviations[replicate_row, ] <- theta_replicate - point
        replicate_completed <- replicate_completed + 1L
        affected_replicates_completed[affected_specifications] <-
          affected_replicates_completed[affected_specifications] + 1L
        replicate_row <- replicate_row + 1L
      }
    }
  }
  if (replicate_row - 1L != replicate_total || replicate_attempted != replicate_total ||
      replicate_completed != replicate_total || support_rebuild_attempted != replicate_total ||
      support_rebuild_completed != replicate_total || any(!is.finite(deviations)) ||
      any(affected_replicates_attempted != expected_affected_replicates) ||
      any(affected_replicates_completed != expected_affected_replicates) ||
      any(surface_refits_attempted != expected_affected_replicates) ||
      any(surface_refits_completed != expected_affected_replicates)) {
    stop("Incomplete or nonfinite fail-fast stacked replicate execution")
  }
  variance <- crossprod(deviations) / 2
  dimnames(variance) <- list(names(point), names(point))

  replicate_diagnostics <- do.call(rbind, lapply(specifications, function(specification) {
    data.frame(
      draw = draw,
      specification = specification,
      surface = standardization_specifications[[specification]]$surface,
      stacked_replicates_expected = replicate_total,
      stacked_replicates_attempted = replicate_attempted,
      stacked_replicates_completed = replicate_completed,
      affected_replicates_expected = expected_affected_replicates[specification],
      affected_replicates_attempted = affected_replicates_attempted[specification],
      affected_replicates_completed = affected_replicates_completed[specification],
      support_rebuilds_expected = replicate_total,
      support_rebuilds_attempted = support_rebuild_attempted,
      support_rebuilds_completed = support_rebuild_completed,
      shared_surface_refits_expected = expected_affected_replicates[specification],
      shared_surface_refits_attempted = surface_refits_attempted[specification],
      shared_surface_refits_completed = surface_refits_completed[specification],
      refits_with_warnings = refits_with_warnings[specification],
      warning_count_across_refits = warning_count_across_refits[specification],
      refits_with_any_pooled_clipping = refits_with_any_pooled_clipping[specification],
      total_pooled_clipped_prediction_rows_across_refits = pooled_clipped_rows[specification],
      maximum_pooled_clipped_prediction_weighted_share = maximum_pooled_clipped_share[specification],
      refits_with_any_own_system_clipping = refits_with_any_own_clipping[specification],
      total_own_system_clipped_prediction_rows_across_refits = own_clipped_rows[specification],
      maximum_own_system_clipped_prediction_weighted_share = maximum_own_clipped_share[specification],
      refits_with_any_common_target_clipping = refits_with_any_target_clipping[specification],
      total_common_target_clipped_prediction_rows_across_refits = target_clipped_rows[specification],
      maximum_common_target_clipped_prediction_weighted_share = maximum_target_clipped_share[specification],
      stringsAsFactors = FALSE
    )
  }))

  list(
    estimate = point,
    variance = variance,
    support_diagnostics = support_diagnostics,
    model_diagnostics = full_surfaces$diagnostics,
    coefficients = full_surfaces$coefficients,
    replicate_diagnostics = replicate_diagnostics,
    replicate_expected = replicate_total,
    replicate_attempted = replicate_attempted,
    replicate_completed = replicate_completed
  )
}

extract_primary_combined <- function(combined, systems) {
  source_names <- c(
    vapply(systems, raw_parameter, character(1)),
    vapply(systems, function(system) {
      component_parameter("primary", "standardized_common", system)
    }, character(1))
  )
  target_names <- c(paste0("raw::", systems), paste0("standardized::", systems))
  extract_vector <- function(value) {
    output <- value[source_names]
    names(output) <- target_names
    output
  }
  extract_matrix <- function(value) {
    output <- value[source_names, source_names, drop = FALSE]
    dimnames(output) <- list(target_names, target_names)
    output
  }
  draw_estimates <- combined$draw_estimates[, source_names, drop = FALSE]
  colnames(draw_estimates) <- target_names
  list(
    estimate = extract_vector(combined$estimate),
    within = extract_matrix(combined$within),
    between = extract_matrix(combined$between),
    variance = extract_matrix(combined$variance),
    se = extract_vector(combined$se),
    draw_estimates = draw_estimates
  )
}

linear_summary <- function(combined, coefficients) {
  coefficients <- coefficients[names(coefficients) %in% names(combined$estimate)]
  full <- setNames(numeric(length(combined$estimate)), names(combined$estimate))
  full[names(coefficients)] <- coefficients
  estimate <- sum(full * combined$estimate)
  variance <- as.numeric(t(full) %*% combined$variance %*% full)
  if (!is.finite(variance) || variance < -1e-10) stop("Invalid contrast variance")
  standard_error <- sqrt(max(0, variance))
  c(
    estimate = estimate,
    standard_error = standard_error,
    ci_lower = estimate - 1.96 * standard_error,
    ci_upper = estimate + 1.96 * standard_error
  )
}

message("Running frozen common-severity specifications with stacked cross-specification JRR")
g45 <- readRDS(file.path(paths$derived, "timss_l_g45_analysis.rds"))
systems <- config$g45_systems
if (!setequal(systems, unique(g45$system))) stop("Configured systems do not match data")
previous_path <- file.path(paths$derived, "standardization_results.rds")
previous_standardization <- if (file.exists(previous_path)) readRDS(previous_path) else NULL

draw_results <- lapply(seq_len(config$pv_count), function(draw) {
  message("  plausible value draw ", draw, " of ", config$pv_count)
  standardize_draw(g45, draw, config$primary_benchmark, systems, config$common_support_quantiles)
})
combined_all <- combine_pv_results(draw_results)
combined_primary <- extract_primary_combined(combined_all, systems)
support_diagnostics <- do.call(rbind, lapply(draw_results, `[[`, "support_diagnostics"))
model_diagnostics <- do.call(rbind, lapply(draw_results, `[[`, "model_diagnostics"))
coefficient_diagnostics <- do.call(rbind, lapply(draw_results, `[[`, "coefficients"))
replicate_diagnostics <- do.call(rbind, lapply(draw_results, `[[`, "replicate_diagnostics"))
rownames(support_diagnostics) <- rownames(model_diagnostics) <-
  rownames(coefficient_diagnostics) <- rownames(replicate_diagnostics) <- NULL

specification_rows <- list()
row <- 1L
for (specification in names(standardization_specifications)) {
  for (system in specification_systems(specification, systems)) {
    raw_full_name <- raw_parameter(system)
    raw_support_name <- component_parameter(specification, "raw_support", system)
    model_own_name <- component_parameter(specification, "model_own_support", system)
    standardized_name <- component_parameter(specification, "standardized_common", system)
    support_component <- linear_summary(combined_all, setNames(c(-1, 1), c(raw_full_name, raw_support_name)))
    calibration_component <- linear_summary(combined_all, setNames(c(-1, 1), c(raw_support_name, model_own_name)))
    target_component <- linear_summary(combined_all, setNames(c(-1, 1), c(model_own_name, standardized_name)))
    total_difference <- linear_summary(combined_all, setNames(c(-1, 1), c(raw_full_name, standardized_name)))
    specification_rows[[row]] <- data.frame(
      specification = specification,
      system = system,
      surface = standardization_specifications[[specification]]$surface,
      target = standardization_specifications[[specification]]$target,
      support = standardization_specifications[[specification]]$support,
      raw_full = combined_all$estimate[raw_full_name],
      raw_full_se = combined_all$se[raw_full_name],
      raw_support = combined_all$estimate[raw_support_name],
      raw_support_se = combined_all$se[raw_support_name],
      model_own_support = combined_all$estimate[model_own_name],
      model_own_support_se = combined_all$se[model_own_name],
      standardized_common = combined_all$estimate[standardized_name],
      standardized_common_se = combined_all$se[standardized_name],
      support_restriction_component = support_component["estimate"],
      support_restriction_se = support_component["standard_error"],
      model_calibration_component = calibration_component["estimate"],
      model_calibration_se = calibration_component["standard_error"],
      target_distribution_component = target_component["estimate"],
      target_distribution_se = target_component["standard_error"],
      standardized_minus_raw_full = total_difference["estimate"],
      standardized_minus_raw_full_se = total_difference["standard_error"],
      total_ci_lower = total_difference["ci_lower"],
      total_ci_upper = total_difference["ci_upper"],
      decomposition_residual = total_difference["estimate"] -
        sum(support_component["estimate"], calibration_component["estimate"], target_component["estimate"]),
      stringsAsFactors = FALSE
    )
    row <- row + 1L
  }
}
specification_results <- do.call(rbind, specification_rows)
rownames(specification_results) <- NULL

alternative_rows <- list()
row <- 1L
for (specification in setdiff(names(standardization_specifications), "primary")) {
  for (system in specification_systems(specification, systems)) {
    alternative_name <- component_parameter(specification, "standardized_common", system)
    primary_name <- component_parameter("primary", "standardized_common", system)
    difference <- linear_summary(
      combined_all, setNames(c(-1, 1), c(primary_name, alternative_name))
    )
    alternative_rows[[row]] <- data.frame(
      specification = specification,
      system = system,
      alternative = combined_all$estimate[alternative_name],
      primary = combined_all$estimate[primary_name],
      alternative_minus_primary = difference["estimate"],
      standard_error = difference["standard_error"],
      ci_lower = difference["ci_lower"],
      ci_upper = difference["ci_upper"],
      stringsAsFactors = FALSE
    )
    row <- row + 1L
  }
}
alternative_contrasts <- do.call(rbind, alternative_rows)
rownames(alternative_contrasts) <- NULL

primary_output <- specification_results[specification_results$specification == "primary", ]
primary_output <- primary_output[match(systems, primary_output$system), ]
output <- data.frame(
  system = primary_output$system,
  threshold = config$primary_benchmark,
  raw = primary_output$raw_full,
  raw_se = primary_output$raw_full_se,
  standardized = primary_output$standardized_common,
  standardized_se = primary_output$standardized_common_se,
  standardized_minus_raw = primary_output$standardized_minus_raw_full,
  difference_se = primary_output$standardized_minus_raw_full_se,
  difference_ci_lower = primary_output$total_ci_lower,
  difference_ci_upper = primary_output$total_ci_upper,
  stringsAsFactors = FALSE
)

pairwise <- list()
row <- 1L
for (first in seq_len(length(systems) - 1L)) {
  for (second in (first + 1L):length(systems)) {
    first_name <- component_parameter("primary", "standardized_common", systems[first])
    second_name <- component_parameter("primary", "standardized_common", systems[second])
    contrast <- linear_summary(combined_all, setNames(c(1, -1), c(first_name, second_name)))
    pairwise[[row]] <- data.frame(
      system_1 = systems[first],
      system_2 = systems[second],
      difference = contrast["estimate"],
      standard_error = contrast["standard_error"],
      z = contrast["estimate"] / contrast["standard_error"],
      p_value = 2 * stats::pnorm(-abs(contrast["estimate"] / contrast["standard_error"])),
      stringsAsFactors = FALSE
    )
    row <- row + 1L
  }
}
pairwise <- do.call(rbind, pairwise)
pairwise$p_holm <- stats::p.adjust(pairwise$p_value, method = "holm")

primary_support <- support_diagnostics[support_diagnostics$specification == "primary", ]
primary_model <- model_diagnostics[model_diagnostics$surface == "quadratic_interaction", c(
  "draw", "system", "model_iterations"
)]
primary_support <- merge(primary_support, primary_model, by = c("draw", "system"), all.x = TRUE, sort = FALSE)
replicate_expected <- vapply(draw_results, `[[`, integer(1), "replicate_expected")
replicate_attempted <- vapply(draw_results, `[[`, integer(1), "replicate_attempted")
replicate_completed <- vapply(draw_results, `[[`, integer(1), "replicate_completed")
primary_support$replicate_expected <- replicate_expected[match(primary_support$draw, seq_len(config$pv_count))]
primary_support$replicate_attempted <- replicate_attempted[match(primary_support$draw, seq_len(config$pv_count))]
primary_support$replicate_completed <- replicate_completed[match(primary_support$draw, seq_len(config$pv_count))]
legacy_support <- primary_support[, c(
  "draw", "system", "risk_n", "support_retained_weight", "support_effective_n",
  "math_lower", "math_upper", "science_lower", "science_upper", "model_iterations",
  "replicate_expected", "replicate_attempted", "replicate_completed"
)]
legacy_support <- legacy_support[order(legacy_support$draw, match(legacy_support$system, systems)), ]
rownames(legacy_support) <- NULL

# Validate all inferential objects before replacing canonical artifacts.
primary <- readRDS(file.path(paths$derived, "primary_transition_results.rds"))
primary_raw <- vapply(systems, function(system) primary$g45[[system]]$estimate["BB_to_AA"], numeric(1))
primary_raw_se <- vapply(systems, function(system) primary$g45[[system]]$se["BB_to_AA"], numeric(1))
if (max(abs(primary_raw - output$raw)) > 1e-10) {
  stop("Raw recovery point estimates differ between primary and standardization pipelines")
}
if (max(abs(primary_raw_se - output$raw_se)) > 1e-10) {
  stop("Raw recovery standard errors differ between primary and standardization pipelines")
}
if (!is.null(previous_standardization)) {
  frozen_members <- c("estimate", "within", "between", "variance", "se", "draw_estimates")
  if (is.null(previous_standardization$combined) ||
      !identical(names(previous_standardization$combined), frozen_members)) {
    stop("Prior standardization artifact lacks the exact frozen primary combined contract")
  }
  if (!identical(names(combined_primary), frozen_members)) {
    stop("Current primary standardization combined contract is not exact")
  }
  for (member in frozen_members) {
    previous_value <- previous_standardization$combined[[member]]
    current_value <- combined_primary[[member]]
    structure_matches <- identical(dim(previous_value), dim(current_value)) &&
      identical(names(previous_value), names(current_value)) &&
      identical(dimnames(previous_value), dimnames(current_value))
    if (!structure_matches || length(previous_value) != length(current_value) ||
        any(!is.finite(previous_value)) || any(!is.finite(current_value))) {
      stop("Prior primary standardization member is incompatible: ", member)
    }
    maximum_difference <- if (length(current_value)) {
      max(abs(as.numeric(previous_value) - as.numeric(current_value)))
    } else {
      0
    }
    if (!is.finite(maximum_difference) || maximum_difference > 1e-10) {
      stop("Refactor changed frozen primary standardization member ", member)
    }
  }
}
if (any(!is.finite(combined_all$estimate)) || any(combined_all$estimate < 0 | combined_all$estimate > 1)) {
  stop("A raw, supported, modeled, or standardized recovery estimate is invalid")
}
if (any(!is.finite(combined_all$variance)) ||
    max(abs(combined_all$variance - t(combined_all$variance))) > 1e-10) {
  stop("The combined cross-specification covariance matrix is invalid or asymmetric")
}
minimum_eigenvalue <- min(eigen(
  (combined_all$variance + t(combined_all$variance)) / 2,
  symmetric = TRUE,
  only.values = TRUE
)$values)
if (minimum_eigenvalue < -1e-8) {
  stop("The combined cross-specification covariance matrix is not positive semidefinite")
}
if (max(abs(specification_results$decomposition_residual)) > 1e-12) {
  stop("Standardization component decomposition failed its identity check")
}
model_clipping_share_columns <- c(
  "pooled_clipped_prediction_weighted_share",
  "own_system_clipped_prediction_weighted_share",
  "common_target_clipped_prediction_weighted_share"
)
model_clipping_count_columns <- c(
  "pooled_clipped_prediction_n", "own_system_clipped_prediction_n",
  "common_target_clipped_prediction_n"
)
if (any(!model_diagnostics$converged) ||
    any(model_diagnostics$model_rank != model_diagnostics$design_columns) ||
    any(!is.finite(model_diagnostics$condition_number)) ||
    any(!is.finite(model_diagnostics$maximum_absolute_coefficient)) ||
    any(model_diagnostics$warning_count < 0) ||
    any(as.matrix(model_diagnostics[model_clipping_count_columns]) < 0) ||
    any(!is.finite(as.matrix(model_diagnostics[model_clipping_share_columns]))) ||
    any(as.matrix(model_diagnostics[model_clipping_share_columns]) < 0 |
      as.matrix(model_diagnostics[model_clipping_share_columns]) > 1) ||
    any(model_diagnostics$predicted_probability_minimum < 0 |
      model_diagnostics$predicted_probability_maximum > 1)) {
  stop("At least one full-sample response-surface diagnostic failed")
}
strict_diagnostics <- support_diagnostics[
  support_diagnostics$specification == "strict_common_cells",
]
if (any(strict_diagnostics$grid_cells_retained <= 0) ||
    any(strict_diagnostics$minimum_common_cell_n < config$strict_cell_min_n) ||
    any(strict_diagnostics$minimum_common_cell_ess < config$strict_cell_min_ess) ||
    any(strict_diagnostics$target_mass_retained_before_renormalization <= 0 |
      strict_diagnostics$target_mass_retained_before_renormalization > 1)) {
  stop("Strict common-cell diagnostics failed their frozen criteria")
}
if (any(abs(support_diagnostics$target_weight_sum - 1) > 1e-10) ||
    any(support_diagnostics$support_retained_weight <= 0 |
      support_diagnostics$support_retained_weight > 1) ||
    any(support_diagnostics$support_effective_n <= 0)) {
  stop("A support or target diagnostic is invalid")
}
replicate_execution_complete <- all(
  replicate_diagnostics$stacked_replicates_attempted == replicate_diagnostics$stacked_replicates_expected &
    replicate_diagnostics$stacked_replicates_completed == replicate_diagnostics$stacked_replicates_expected &
    replicate_diagnostics$affected_replicates_attempted == replicate_diagnostics$affected_replicates_expected &
    replicate_diagnostics$affected_replicates_completed == replicate_diagnostics$affected_replicates_expected &
    replicate_diagnostics$support_rebuilds_attempted == replicate_diagnostics$support_rebuilds_expected &
    replicate_diagnostics$support_rebuilds_completed == replicate_diagnostics$support_rebuilds_expected &
    replicate_diagnostics$shared_surface_refits_attempted == replicate_diagnostics$shared_surface_refits_expected &
    replicate_diagnostics$shared_surface_refits_completed == replicate_diagnostics$shared_surface_refits_expected
)
replicate_clipping_share_columns <- c(
  "maximum_pooled_clipped_prediction_weighted_share",
  "maximum_own_system_clipped_prediction_weighted_share",
  "maximum_common_target_clipped_prediction_weighted_share"
)
replicate_count_columns <- setdiff(
  names(replicate_diagnostics),
  c("draw", "specification", "surface", replicate_clipping_share_columns)
)
replicate_diagnostics_valid <-
  all(is.finite(as.matrix(replicate_diagnostics[replicate_count_columns]))) &&
  all(as.matrix(replicate_diagnostics[replicate_count_columns]) >= 0) &&
  all(is.finite(as.matrix(replicate_diagnostics[replicate_clipping_share_columns]))) &&
  all(as.matrix(replicate_diagnostics[replicate_clipping_share_columns]) >= 0 &
    as.matrix(replicate_diagnostics[replicate_clipping_share_columns]) <= 1)
if (!replicate_execution_complete || !replicate_diagnostics_valid) {
  stop("At least one fail-fast stacked standardization stage or diagnostic is incomplete")
}
expected_specification_rows <- sum(vapply(
  names(standardization_specifications),
  function(specification) length(specification_systems(specification, systems)),
  integer(1)
))
if (nrow(specification_results) != expected_specification_rows ||
    nrow(support_diagnostics) != config$pv_count * expected_specification_rows ||
    nrow(model_diagnostics) != config$pv_count * length(systems) * length(unique(vapply(
      standardization_specifications, `[[`, character(1), "surface"
    ))) ||
    nrow(alternative_contrasts) != expected_specification_rows - length(systems) ||
    nrow(pairwise) != choose(length(systems), 2) ||
    nrow(legacy_support) != config$pv_count * length(systems)) {
  stop("Unexpected standardization artifact dimensions")
}

covariance_validation <- data.frame(
  check = c(
    "primary raw estimates reproduce transition pipeline",
    "primary raw SEs reproduce transition pipeline",
    "all component estimates are finite probabilities",
    "cross-specification covariance is finite and symmetric",
    "cross-specification covariance is positive semidefinite",
    "component decomposition identity holds",
    "all response surfaces converged at full rank",
    "strict common cells satisfy n and ESS criteria",
    "all targets sum to one",
    "all fail-fast replicate stages and diagnostics completed",
    "artifact dimensions match frozen specifications"
  ),
  passed = c(
    max(abs(primary_raw - output$raw)) <= 1e-10,
    max(abs(primary_raw_se - output$raw_se)) <= 1e-10,
    all(is.finite(combined_all$estimate) & combined_all$estimate >= 0 & combined_all$estimate <= 1),
    all(is.finite(combined_all$variance)) &&
      max(abs(combined_all$variance - t(combined_all$variance))) <= 1e-10,
    minimum_eigenvalue >= -1e-8,
    max(abs(specification_results$decomposition_residual)) <= 1e-12,
    all(model_diagnostics$converged & model_diagnostics$model_rank == model_diagnostics$design_columns),
    all(strict_diagnostics$minimum_common_cell_n >= config$strict_cell_min_n &
      strict_diagnostics$minimum_common_cell_ess >= config$strict_cell_min_ess),
    all(abs(support_diagnostics$target_weight_sum - 1) <= 1e-10),
    replicate_execution_complete && replicate_diagnostics_valid,
    nrow(specification_results) == expected_specification_rows
  ),
  value = c(
    format(max(abs(primary_raw - output$raw)), scientific = TRUE),
    format(max(abs(primary_raw_se - output$raw_se)), scientific = TRUE),
    paste(length(combined_all$estimate), "parameters"),
    format(max(abs(combined_all$variance - t(combined_all$variance))), scientific = TRUE),
    format(minimum_eigenvalue, scientific = TRUE),
    format(max(abs(specification_results$decomposition_residual)), scientific = TRUE),
    paste(nrow(model_diagnostics), "draw-system-surface fits"),
    paste0(
      min(strict_diagnostics$grid_cells_retained), "-", max(strict_diagnostics$grid_cells_retained),
      " of ", config$strict_cell_grid^2, " cells"
    ),
    format(max(abs(support_diagnostics$target_weight_sum - 1)), scientific = TRUE),
    paste(
      sum(replicate_completed), "of", sum(replicate_expected),
      "stacked replicates completed"
    ),
    paste(expected_specification_rows, "specification-system rows")
  ),
  stringsAsFactors = FALSE
)
if (!all(covariance_validation$passed)) stop("Internal standardization validation table contains a failure")

standardization_results <- list(
  version = "2.1-auditability",
  threshold = config$primary_benchmark,
  systems = systems,
  specifications = standardization_specifications,
  definitions = list(
    matched_pv = "Common terminal PV index across mathematics, science, baseline, and follow-up",
    common_support_probabilities = config$common_support_quantiles,
    fail_fast_replicate_diagnostics = paste(
      "Expected, attempted, and completed counts are retained; any failed stage",
      "prevents canonical artifact replacement"
    ),
    clipping_domains = c("pooled evaluation", "own system", "primary/common specification target"),
    strict_cells = paste0(
      config$strict_cell_grid, "x", config$strict_cell_grid,
      " target-quartile grid; every system n >= ", config$strict_cell_min_n,
      " and Kish ESS >= ", config$strict_cell_min_ess
    ),
    components = c("raw_full", "raw_support", "model_own_support", "standardized_common")
  ),
  combined = combined_primary,
  combined_all = combined_all,
  diagnostics = legacy_support,
  specification_diagnostics = support_diagnostics,
  model_diagnostics = model_diagnostics,
  replicate_diagnostics = replicate_diagnostics,
  validation = covariance_validation
)

(function() {
  stage_root <- tempfile(pattern = ".standardization-stage-", tmpdir = project_root)
  created <- dir.create(stage_root, recursive = TRUE, showWarnings = FALSE)
  if (!created && !dir.exists(stage_root)) {
    stop("Could not create standardization staging directory: ", stage_root)
  }
  on.exit(if (dir.exists(stage_root)) unlink(stage_root, recursive = TRUE, force = TRUE), add = TRUE)

staged <- list(
  results = file.path(stage_root, "standardization_results.rds"),
  primary_table = file.path(stage_root, "raw_vs_standardized_recovery.csv"),
  specification_table = file.path(stage_root, "standardization_specification_results.csv"),
  alternative_table = file.path(stage_root, "standardization_alternative_contrasts.csv"),
  pairwise = file.path(stage_root, "standardized_pairwise_contrasts_exploratory.csv"),
  legacy_support = file.path(stage_root, "standardization_support_by_draw.csv"),
  support = file.path(stage_root, "standardization_specification_diagnostics_by_draw.csv"),
  models = file.path(stage_root, "standardization_model_diagnostics_by_draw.csv"),
  coefficients = file.path(stage_root, "standardization_coefficients_by_draw.csv"),
  replicates = file.path(stage_root, "standardization_replicate_diagnostics.csv"),
  validation = file.path(stage_root, "standardization_covariance_validation.csv")
)
saveRDS(standardization_results, staged$results, compress = "xz")
write.csv(output, staged$primary_table, row.names = FALSE)
write.csv(specification_results, staged$specification_table, row.names = FALSE)
write.csv(alternative_contrasts, staged$alternative_table, row.names = FALSE)
write.csv(pairwise, staged$pairwise, row.names = FALSE)
write.csv(legacy_support, staged$legacy_support, row.names = FALSE)
write.csv(support_diagnostics, staged$support, row.names = FALSE)
write.csv(model_diagnostics, staged$models, row.names = FALSE)
write.csv(coefficient_diagnostics, staged$coefficients, row.names = FALSE)
write.csv(replicate_diagnostics, staged$replicates, row.names = FALSE)
write.csv(covariance_validation, staged$validation, row.names = FALSE)

destinations <- list(
  results = file.path(paths$derived, "standardization_results.rds"),
  primary_table = file.path(paths$tables, "raw_vs_standardized_recovery.csv"),
  specification_table = file.path(paths$tables, "standardization_specification_results.csv"),
  alternative_table = file.path(paths$tables, "standardization_alternative_contrasts.csv"),
  pairwise = file.path(paths$tables, "standardized_pairwise_contrasts_exploratory.csv"),
  legacy_support = file.path(paths$diagnostics, "standardization_support_by_draw.csv"),
  support = file.path(paths$diagnostics, "standardization_specification_diagnostics_by_draw.csv"),
  models = file.path(paths$diagnostics, "standardization_model_diagnostics_by_draw.csv"),
  coefficients = file.path(paths$diagnostics, "standardization_coefficients_by_draw.csv"),
  replicates = file.path(paths$diagnostics, "standardization_replicate_diagnostics.csv"),
  validation = file.path(paths$diagnostics, "standardization_covariance_validation.csv")
)
for (name in names(staged)) replace_file(staged[[name]], destinations[[name]])
message(
  "Standardization complete: ", length(combined_all$estimate),
  " jointly estimated parameters across ", length(standardization_specifications),
  " specifications; strict support retained ",
  paste(sort(unique(strict_diagnostics$grid_cells_retained)), collapse = ", "),
  " of ", config$strict_cell_grid^2, " cells; all ", nrow(covariance_validation),
  " validation checks passed"
)
})()
