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

md5_file <- function(path) {
  value <- unname(tools::md5sum(path))
  if (length(value) != 1L || is.na(value)) stop("Could not calculate MD5 for ", path)
  tolower(value)
}

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

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

# Refuse to build derived data from an unreviewed or corrupted source cache.
source_paths <- file.path(paths$raw, raw_files)
missing_sources <- names(source_paths)[!file.exists(source_paths)]
if (length(missing_sources)) {
  stop("Official source files are missing; run analysis/01_acquire_data.R: ", paste(missing_sources, collapse = ", "))
}
observed_source_md5 <- unname(tools::md5sum(source_paths))
names(observed_source_md5) <- names(source_paths)
source_mismatch <- names(observed_source_md5)[
  tolower(observed_source_md5) != tolower(expected_source_md5[names(observed_source_md5)])
]
if (length(source_mismatch)) {
  stop("Official source fingerprint mismatch: ", paste(source_mismatch, collapse = ", "))
}
for (key in c("data_r_g45", "data_r_g89")) {
  directory <- file.path(paths$raw, if (key == "data_r_g45") "r_g45" else "r_g89")
  marker <- file.path(directory, ".extraction_manifest.csv")
  if (!file.exists(marker)) stop("Extraction lineage marker is missing; rerun acquisition: ", marker)
  lineage <- read.csv(marker, stringsAsFactors = FALSE)
  if (!nrow(lineage) || any(tolower(lineage$archive_md5) != tolower(expected_source_md5[[key]]))) {
    stop("Extraction lineage does not match its frozen archive: ", key)
  }
}

xlsx_lineage_rows <- list()
xlsx_lineage_index <- 1L

xlsx_manifest <- function(directory, workbook, workbook_md5, extractor_md5) {
  files <- sort(list.files(directory, pattern = "[.]csv$", recursive = TRUE, full.names = TRUE))
  files <- files[basename(files) != ".cache_manifest.csv"]
  if (!length(files)) stop("Workbook extraction produced no CSV files: ", workbook)
  prefix <- paste0(normalizePath(directory, mustWork = TRUE), .Platform$file.sep)
  relative <- substring(normalizePath(files, mustWork = TRUE), nchar(prefix) + 1L)
  data.frame(
    workbook_filename = basename(workbook),
    workbook_md5 = workbook_md5,
    extractor_filename = "xlsx_extract.py",
    extractor_md5 = extractor_md5,
    relative_path = relative,
    bytes = unname(file.info(files)$size),
    md5 = unname(tools::md5sum(files)),
    stringsAsFactors = FALSE
  )
}

xlsx_cache_valid <- function(directory, workbook, workbook_md5, extractor_md5, expected_sheet) {
  marker <- file.path(directory, ".cache_manifest.csv")
  expected <- file.path(directory, expected_sheet)
  if (!dir.exists(directory) || !file.exists(marker) || !file.exists(expected)) return(FALSE)
  recorded <- tryCatch(read.csv(marker, stringsAsFactors = FALSE), error = function(error) data.frame())
  required <- c(
    "workbook_filename", "workbook_md5", "extractor_filename", "extractor_md5",
    "relative_path", "bytes", "md5"
  )
  if (!nrow(recorded) || !all(required %in% names(recorded))) return(FALSE)
  if (any(tolower(recorded$workbook_md5) != workbook_md5) ||
      any(tolower(recorded$extractor_md5) != extractor_md5)) return(FALSE)
  current <- tryCatch(
    xlsx_manifest(directory, workbook, workbook_md5, extractor_md5),
    error = function(error) data.frame()
  )
  if (!nrow(current) || !identical(recorded$relative_path, current$relative_path)) return(FALSE)
  identical(as.numeric(recorded$bytes), as.numeric(current$bytes)) &&
    identical(tolower(recorded$md5), tolower(current$md5))
}

ensure_report_csv <- function(workbook, sheet_filename) {
  extractor <- file.path(analysis_dir, "xlsx_extract.py")
  workbook_md5 <- md5_file(workbook)
  extractor_md5 <- md5_file(extractor)
  directory <- file.path(paths$derived, "xlsx", tools::file_path_sans_ext(basename(workbook)))
  expected <- file.path(directory, sheet_filename)
  status <- "cached_verified"

  if (!xlsx_cache_valid(directory, workbook, workbook_md5, extractor_md5, sheet_filename)) {
    parent <- dirname(directory)
    dir.create(parent, recursive = TRUE, showWarnings = FALSE)
    staged <- tempfile(pattern = paste0(".", basename(directory), "-"), tmpdir = parent)
    dir.create(staged, recursive = TRUE, showWarnings = FALSE)
    on.exit(if (dir.exists(staged)) unlink(staged, recursive = TRUE, force = TRUE), add = TRUE)
    output <- system2(
      "python3",
      c(shQuote(extractor), shQuote(workbook), shQuote(staged)),
      stdout = TRUE,
      stderr = TRUE
    )
    command_status <- attr(output, "status")
    if (is.null(command_status)) command_status <- 0L
    staged_expected <- file.path(staged, sheet_filename)
    if (command_status != 0L || !file.exists(staged_expected)) {
      stop(
        "Could not extract expected worksheet ", sheet_filename, " from ", basename(workbook),
        if (length(output)) paste0("\n", paste(output, collapse = "\n")) else ""
      )
    }
    manifest <- xlsx_manifest(staged, workbook, workbook_md5, extractor_md5)
    write.csv(manifest, file.path(staged, ".cache_manifest.csv"), row.names = FALSE)
    replace_directory(staged, directory)
    status <- "refreshed_verified"
  }

  xlsx_lineage_rows[[xlsx_lineage_index]] <<- data.frame(
    workbook_filename = basename(workbook),
    workbook_md5 = workbook_md5,
    extractor_md5 = extractor_md5,
    sheet_filename = sheet_filename,
    sheet_md5 = md5_file(expected),
    cache_status = status,
    stringsAsFactors = FALSE
  )
  xlsx_lineage_index <<- xlsx_lineage_index + 1L
  expected
}

mean_report_file <- function(subject, band) {
  key <- paste0("report_", subject, "_", band)
  sheet <- c(
    report_math_g45 = "01_Exhibit_1.1.1a.csv",
    report_science_g45 = "01_Exhibit_2.1.1a.csv",
    report_math_g89 = "01_Exhibit_1.2.1a.csv",
    report_science_g89 = "01_Exhibit_2.2.1a.csv"
  )[[key]]
  if (is.null(sheet)) stop("Unknown official mean report: ", key)
  ensure_report_csv(file.path(paths$raw, raw_files[[key]]), sheet)
}

benchmark_report_file <- function(subject) {
  key <- paste0("report_", subject, "_benchmark_g45")
  sheet <- c(
    report_math_benchmark_g45 = "01_Exhibit_1.1.3.csv",
    report_science_benchmark_g45 = "01_Exhibit_2.1.3.csv"
  )[[key]]
  if (is.null(sheet)) stop("Unknown official benchmark report: ", key)
  ensure_report_csv(file.path(paths$raw, raw_files[[key]]), sheet)
}

parse_mean_report <- function(file, subject, band) {
  x <- read.csv(file, header = FALSE, fill = TRUE, check.names = FALSE, quote = "\"", comment.char = "")
  country <- trimws(x[[3]])
  keep <- country %in% unname(system_labels)
  data.frame(
    system = country[keep],
    subject = subject,
    band = band,
    official_t0 = suppressWarnings(as.numeric(x[[5]][keep])),
    official_t0_se = suppressWarnings(as.numeric(x[[6]][keep])),
    official_t1 = suppressWarnings(as.numeric(x[[8]][keep])),
    official_t1_se = suppressWarnings(as.numeric(x[[9]][keep])),
    official_change = suppressWarnings(as.numeric(x[[11]][keep])),
    official_change_se = suppressWarnings(as.numeric(x[[12]][keep])),
    stringsAsFactors = FALSE
  )
}

parse_benchmark_report <- function(file, subject) {
  x <- read.csv(file, header = FALSE, fill = TRUE, check.names = FALSE, quote = "\"", comment.char = "")
  token <- trimws(x[[3]])
  current_system <- NA_character_
  output <- list()
  index <- 1L
  for (row in seq_len(nrow(x))) {
    if (token[row] %in% unname(system_labels)) current_system <- token[row]
    if (token[row] %in% c("Grade 4", "Grade 5") && !is.na(current_system)) {
      output[[index]] <- data.frame(
        system = current_system,
        subject = subject,
        band = "g45",
        wave = ifelse(token[row] == "Grade 4", "t0", "t1"),
        official_at_or_above_400 = suppressWarnings(as.numeric(x[[15]][row])) / 100,
        official_at_or_above_400_se = suppressWarnings(as.numeric(x[[16]][row])) / 100,
        benchmark_se_suppressed = trimws(x[[16]][row]) == "~",
        stringsAsFactors = FALSE
      )
      index <- index + 1L
    }
  }
  do.call(rbind, output)
}

message("Loading and reducing official Grade 4-5 files")
g45 <- load_band_data(file.path(paths$raw, "r_g45"), "g45")
message("Loading and reducing official Grade 8-9 files")
g89 <- load_band_data(file.path(paths$raw, "r_g89"), "g89")

validate_reduced_design <- function(data, label) {
  design <- as.matrix(data[c("system_id", "school_id", "student_id", "weight", "jkzone", "jkrep")])
  if (any(!is.finite(design)) || any(data$weight <= 0)) {
    stop(label, ": nonfinite design values or nonpositive final weights")
  }
  if (any(!(data$jkrep %in% c(0, 1)))) stop(label, ": JKREP is not exactly binary")
  if (any(!(data$participated_2023 %in% c(0L, 1L))) ||
      any(!(data$participated_2024 %in% c(0L, 1L)))) {
    stop(label, ": derived participation indicators are not binary")
  }
  expected_pattern <- ifelse(
    data$participated_2023 == 1 & data$participated_2024 == 1, "both",
    ifelse(data$participated_2023 == 1, "2023_only",
      ifelse(data$participated_2024 == 1, "2024_only", "neither"))
  )
  if (!identical(data$participation_pattern, expected_pattern)) {
    stop(label, ": participation-pattern labels do not match binary wave indicators")
  }
  members <- split(data$jkrep, paste(data$system_id, data$jkzone, sep = ":"))
  if (!all(vapply(members, function(value) setequal(unique(value), c(0, 1)), logical(1)))) {
    stop(label, ": at least one system-zone lacks JKREP member 0 or 1")
  }
  invisible(TRUE)
}
validate_reduced_design(g45, "Grade 4-5")
validate_reduced_design(g89, "Grade 8-9")

make_inventory <- function(data) {
  systems <- split(data, data$system)
  do.call(rbind, lapply(systems, function(d) {
    data.frame(
      band = unique(d$band),
      system = unique(d$system),
      students = nrow(d),
      schools = length(unique(d$school_id)),
      zones = length(unique(d$jkzone)),
      weighted_population = sum(d$weight),
      both_wave_n = sum(d$participation_pattern == "both"),
      baseline_only_n = sum(d$participation_pattern == "2023_only"),
      followup_only_n = sum(d$participation_pattern == "2024_only"),
      neither_n = sum(d$participation_pattern == "neither"),
      both_wave_weighted_pct = 100 * weighted_probability(d$participation_pattern == "both", d$weight),
      pv_missing_cells = sum(!is.finite(as.matrix(d[all_pv_variables(unique(d$band))]))),
      stringsAsFactors = FALSE
    )
  }))
}
inventory <- rbind(make_inventory(g45), make_inventory(g89))
rownames(inventory) <- NULL

calculate_summary <- function(data, band) {
  systems <- split(data, data$system)
  rows <- list()
  row <- 1L
  for (system in names(systems)) {
    d <- systems[[system]]
    draw_results <- lapply(seq_len(config$pv_count), function(draw) {
      map <- pv_variable_map(band, draw)
      estimator <- function(dataset, weight) {
        c(
          math_t0 = weighted_mean(dataset[[map$math_t0]], weight),
          math_t1 = weighted_mean(dataset[[map$math_t1]], weight),
          math_change = weighted_mean(dataset[[map$math_change]], weight),
          math_benchmark_t0 = weighted_probability(dataset[[map$math_t0]] >= config$low_benchmark, weight),
          math_benchmark_t1 = weighted_probability(dataset[[map$math_t1]] >= config$low_benchmark, weight),
          science_t0 = weighted_mean(dataset[[map$science_t0]], weight),
          science_t1 = weighted_mean(dataset[[map$science_t1]], weight),
          science_change = weighted_mean(dataset[[map$science_change]], weight),
          science_benchmark_t0 = weighted_probability(dataset[[map$science_t0]] >= config$low_benchmark, weight),
          science_benchmark_t1 = weighted_probability(dataset[[map$science_t1]] >= config$low_benchmark, weight)
        )
      }
      result <- jrr_estimate(d, estimator)
      result$replicates <- NULL
      result
    })
    combined <- combine_pv_results(draw_results)
    for (subject in c("math", "science")) {
      value <- function(suffix) unname(combined$estimate[paste0(subject, "_", suffix)])
      standard_error <- function(suffix) unname(combined$se[paste0(subject, "_", suffix)])
      rows[[row]] <- data.frame(
        system = system,
        subject = subject,
        band = band,
        estimated_t0 = value("t0"),
        estimated_t0_se = standard_error("t0"),
        estimated_t1 = value("t1"),
        estimated_t1_se = standard_error("t1"),
        estimated_change = value("change"),
        estimated_change_se = standard_error("change"),
        estimated_at_or_above_400_t0 = value("benchmark_t0"),
        estimated_at_or_above_400_t0_se = standard_error("benchmark_t0"),
        estimated_at_or_above_400_t1 = value("benchmark_t1"),
        estimated_at_or_above_400_t1_se = standard_error("benchmark_t1"),
        stringsAsFactors = FALSE
      )
      row <- row + 1L
    }
  }
  do.call(rbind, rows)
}

message("Estimating official marginal points and SEs with matched PVs and JRR")
calculated <- rbind(calculate_summary(g45, "g45"), calculate_summary(g89, "g89"))

official_means <- rbind(
  parse_mean_report(mean_report_file("math", "g45"), "math", "g45"),
  parse_mean_report(mean_report_file("science", "g45"), "science", "g45"),
  parse_mean_report(mean_report_file("math", "g89"), "math", "g89"),
  parse_mean_report(mean_report_file("science", "g89"), "science", "g89")
)
validation <- merge(calculated, official_means, by = c("system", "subject", "band"), all.x = TRUE, sort = FALSE)

benchmarks <- rbind(
  parse_benchmark_report(benchmark_report_file("math"), "math"),
  parse_benchmark_report(benchmark_report_file("science"), "science")
)
for (wave in c("t0", "t1")) {
  wave_rows <- benchmarks[benchmarks$wave == wave, c(
    "system", "subject", "band", "official_at_or_above_400",
    "official_at_or_above_400_se", "benchmark_se_suppressed"
  )]
  names(wave_rows)[4:6] <- paste0(
    c("official_at_or_above_400_", "official_at_or_above_400_se_", "benchmark_se_suppressed_"),
    wave
  )
  validation <- merge(validation, wave_rows, by = c("system", "subject", "band"), all.x = TRUE, sort = FALSE)
}

validation$diff_t0 <- validation$estimated_t0 - validation$official_t0
validation$diff_t0_se <- validation$estimated_t0_se - validation$official_t0_se
validation$diff_t1 <- validation$estimated_t1 - validation$official_t1
validation$diff_t1_se <- validation$estimated_t1_se - validation$official_t1_se
validation$diff_change <- validation$estimated_change - validation$official_change
validation$diff_change_se <- validation$estimated_change_se - validation$official_change_se
validation$diff_benchmark_t0_pp <- 100 * (
  validation$estimated_at_or_above_400_t0 - validation$official_at_or_above_400_t0
)
validation$diff_benchmark_t0_se_pp <- 100 * (
  validation$estimated_at_or_above_400_t0_se - validation$official_at_or_above_400_se_t0
)
validation$diff_benchmark_t1_pp <- 100 * (
  validation$estimated_at_or_above_400_t1 - validation$official_at_or_above_400_t1
)
validation$diff_benchmark_t1_se_pp <- 100 * (
  validation$estimated_at_or_above_400_t1_se - validation$official_at_or_above_400_se_t1
)

point_tolerance_score <- 0.55
point_tolerance_pp <- 0.55
se_tolerance_score <- 0.055
se_tolerance_pp <- 0.055
validation$passes_point_rounding_check <- with(validation,
  abs(diff_t0) <= point_tolerance_score &
    abs(diff_t1) <= point_tolerance_score &
    abs(diff_change) <= point_tolerance_score &
    (band != "g45" | (
      abs(diff_benchmark_t0_pp) <= point_tolerance_pp &
        abs(diff_benchmark_t1_pp) <= point_tolerance_pp
    ))
)
validation$passes_se_rounding_check <- with(validation,
  abs(diff_t0_se) <= se_tolerance_score &
    abs(diff_t1_se) <= se_tolerance_score &
    abs(diff_change_se) <= se_tolerance_score &
    (band != "g45" | (
      (is.na(official_at_or_above_400_se_t0) | abs(diff_benchmark_t0_se_pp) <= se_tolerance_pp) &
        (is.na(official_at_or_above_400_se_t1) | abs(diff_benchmark_t1_se_pp) <= se_tolerance_pp)
    ))
)
validation$passes_rounding_check <- validation$passes_point_rounding_check & validation$passes_se_rounding_check

validation_specs <- list(
  list(value_type = "point", statistic = "mean_t0", estimated = "estimated_t0", official = "official_t0", units = "TIMSS points", tolerance = point_tolerance_score, applies = function(row) TRUE),
  list(value_type = "point", statistic = "mean_t1", estimated = "estimated_t1", official = "official_t1", units = "TIMSS points", tolerance = point_tolerance_score, applies = function(row) TRUE),
  list(value_type = "point", statistic = "change", estimated = "estimated_change", official = "official_change", units = "TIMSS points", tolerance = point_tolerance_score, applies = function(row) TRUE),
  list(value_type = "point", statistic = "at_or_above_400_t0", estimated = "estimated_at_or_above_400_t0", official = "official_at_or_above_400_t0", units = "percentage points", tolerance = point_tolerance_pp, multiplier = 100, applies = function(row) row$band == "g45"),
  list(value_type = "point", statistic = "at_or_above_400_t1", estimated = "estimated_at_or_above_400_t1", official = "official_at_or_above_400_t1", units = "percentage points", tolerance = point_tolerance_pp, multiplier = 100, applies = function(row) row$band == "g45"),
  list(value_type = "standard_error", statistic = "mean_t0_se", estimated = "estimated_t0_se", official = "official_t0_se", units = "TIMSS points", tolerance = se_tolerance_score, applies = function(row) TRUE),
  list(value_type = "standard_error", statistic = "mean_t1_se", estimated = "estimated_t1_se", official = "official_t1_se", units = "TIMSS points", tolerance = se_tolerance_score, applies = function(row) TRUE),
  list(value_type = "standard_error", statistic = "change_se", estimated = "estimated_change_se", official = "official_change_se", units = "TIMSS points", tolerance = se_tolerance_score, applies = function(row) TRUE),
  list(value_type = "standard_error", statistic = "at_or_above_400_t0_se", estimated = "estimated_at_or_above_400_t0_se", official = "official_at_or_above_400_se_t0", units = "percentage points", tolerance = se_tolerance_pp, multiplier = 100, applies = function(row) row$band == "g45"),
  list(value_type = "standard_error", statistic = "at_or_above_400_t1_se", estimated = "estimated_at_or_above_400_t1_se", official = "official_at_or_above_400_se_t1", units = "percentage points", tolerance = se_tolerance_pp, multiplier = 100, applies = function(row) row$band == "g45")
)
validation_long_rows <- list()
long_row <- 1L
for (index in seq_len(nrow(validation))) {
  row <- validation[index, , drop = FALSE]
  for (specification in validation_specs) {
    if (!specification$applies(row)) next
    multiplier <- if (is.null(specification$multiplier)) 1 else specification$multiplier
    estimated <- row[[specification$estimated]] * multiplier
    official <- row[[specification$official]] * multiplier
    reported <- is.finite(official)
    difference <- if (reported) estimated - official else NA_real_
    validation_long_rows[[long_row]] <- data.frame(
      system = row$system,
      subject = row$subject,
      band = row$band,
      value_type = specification$value_type,
      statistic = specification$statistic,
      units = specification$units,
      estimated = estimated,
      official = official,
      difference = difference,
      absolute_tolerance = specification$tolerance,
      official_reported = reported,
      included_in_gate = reported,
      passed = if (reported) abs(difference) <= specification$tolerance else NA,
      stringsAsFactors = FALSE
    )
    long_row <- long_row + 1L
  }
}
validation_long <- do.call(rbind, validation_long_rows)
rownames(validation_long) <- NULL

baseline_state_audit <- do.call(rbind, lapply(split(g45, g45$system), function(d) {
  do.call(rbind, lapply(seq_len(config$pv_count), function(draw) {
    map <- pv_variable_map("g45", draw)
    state <- achievement_state(d[[map$math_t0]], d[[map$science_t0]], config$low_benchmark)
    do.call(rbind, lapply(levels(state), function(level) {
      data.frame(
        system = unique(d$system), draw = draw, state = level,
        unweighted_n = sum(state == level, na.rm = TRUE),
        weighted_pct = 100 * weighted_probability(state == level, d$weight),
        stringsAsFactors = FALSE
      )
    }))
  }))
}))

variable_map <- do.call(rbind, lapply(c("g45", "g89"), function(band) {
  do.call(rbind, lapply(seq_len(config$pv_count), function(draw) {
    data.frame(
      band = band,
      draw = draw,
      role = names(pv_variable_map(band, draw)),
      variable = unlist(pv_variable_map(band, draw)),
      row.names = NULL
    )
  }))
}))

xlsx_lineage <- do.call(rbind, xlsx_lineage_rows)
rownames(xlsx_lineage) <- NULL
xlsx_lineage <- xlsx_lineage[!duplicated(xlsx_lineage[, c("workbook_filename", "sheet_filename")]), ]

# Validate every object before any canonical audit or reduced-data artifact is replaced.
if (nrow(validation) != 24L) stop("Official validation must contain exactly 24 system-subject-band rows")
if (nrow(validation_long[validation_long$value_type == "point", ]) != 108L) {
  stop("Official validation must contain exactly 108 scalar point comparisons")
}
if (nrow(validation_long[validation_long$value_type == "standard_error", ]) != 108L) {
  stop("Official validation must contain exactly 108 scalar SE fields")
}
if (sum(validation_long$value_type == "standard_error" & validation_long$included_in_gate) != 102L) {
  stop("Expected exactly 102 numerically reported official SE comparisons")
}
if (sum(validation_long$value_type == "standard_error" & !validation_long$official_reported) != 6L) {
  stop("Expected exactly six suppressed benchmark SE cells")
}
if (any(inventory$pv_missing_cells != 0)) stop("Unexpected missing plausible values in reduced analysis data")
if (any(!validation$passes_rounding_check)) {
  print(validation[!validation$passes_rounding_check, ])
  stop("Official point/SE reproduction failed rounding tolerance")
}
gated_rows <- validation_long[validation_long$included_in_gate, ]
if (any(!gated_rows$passed) || any(!is.finite(gated_rows$estimated))) {
  stop("At least one scalar official-result validation failed")
}
if (nrow(inventory) != 12L || nrow(variable_map) != 60L || nrow(baseline_state_audit) != 180L) {
  stop("Unexpected audit artifact dimensions")
}

point_rows <- validation_long[validation_long$value_type == "point" & validation_long$included_in_gate, ]
se_rows <- validation_long[validation_long$value_type == "standard_error" & validation_long$included_in_gate, ]
score_point_rows <- point_rows[point_rows$units == "TIMSS points", ]
benchmark_point_rows <- point_rows[point_rows$units == "percentage points", ]
score_se_rows <- se_rows[se_rows$units == "TIMSS points", ]
benchmark_se_rows <- se_rows[se_rows$units == "percentage points", ]
summary_lines <- c(
  "TIMSS 2023 Longitudinal data and official-result audit: PASS",
  "Design gates: finite identifiers/weights/zones/members, positive weights, JKREP in {0,1}, and both members per system-zone: PASS",
  "Participation gates: nonmissing source codes in {1,2}, code 1 used for participation, and derived indicators binary: PASS",
  paste("Grade 4-5 students:", nrow(g45), "across", length(unique(g45$system)), "systems"),
  paste("Grade 8-9 students:", nrow(g89), "across", length(unique(g89$system)), "systems"),
  paste("Official validation rows:", nrow(validation)),
  paste("Published scalar point comparisons passed:", sum(point_rows$passed), "of", nrow(point_rows)),
  paste("Published numeric SE comparisons passed:", sum(se_rows$passed), "of", nrow(se_rows)),
  paste("Suppressed benchmark SE cells recorded but excluded from numeric gate:", sum(!validation_long$official_reported)),
  paste("Maximum absolute mean/growth point discrepancy:", round(max(abs(score_point_rows$difference)), 4), "TIMSS points"),
  paste("Maximum absolute benchmark point discrepancy:", round(max(abs(benchmark_point_rows$difference)), 4), "percentage points"),
  paste("Maximum absolute mean/growth SE discrepancy:", round(max(abs(score_se_rows$difference)), 4), "TIMSS points"),
  paste("Maximum absolute benchmark SE discrepancy:", round(max(abs(benchmark_se_rows$difference)), 4), "percentage points")
)

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

staged <- list(
  g45 = file.path(stage_root, "timss_l_g45_analysis.rds"),
  g89 = file.path(stage_root, "timss_l_g89_analysis.rds"),
  inventory = file.path(stage_root, "data_inventory.csv"),
  validation = file.path(stage_root, "official_result_validation.csv"),
  validation_long = file.path(stage_root, "official_result_validation_long.csv"),
  baseline = file.path(stage_root, "baseline_state_draws.csv"),
  variable_map = file.path(stage_root, "variable_map.csv"),
  xlsx_lineage = file.path(stage_root, "xlsx_cache_lineage.csv"),
  summary = file.path(stage_root, "audit_summary.txt")
)
saveRDS(g45, staged$g45, compress = "xz")
saveRDS(g89, staged$g89, compress = "xz")
write.csv(inventory, staged$inventory, row.names = FALSE)
write.csv(validation, staged$validation, row.names = FALSE)
write.csv(validation_long, staged$validation_long, row.names = FALSE)
write.csv(baseline_state_audit, staged$baseline, row.names = FALSE)
write.csv(variable_map, staged$variable_map, row.names = FALSE)
write.csv(xlsx_lineage, staged$xlsx_lineage, row.names = FALSE)
writeLines(summary_lines, staged$summary)

destinations <- list(
  g45 = file.path(paths$derived, "timss_l_g45_analysis.rds"),
  g89 = file.path(paths$derived, "timss_l_g89_analysis.rds"),
  inventory = file.path(paths$diagnostics, "data_inventory.csv"),
  validation = file.path(paths$diagnostics, "official_result_validation.csv"),
  validation_long = file.path(paths$diagnostics, "official_result_validation_long.csv"),
  baseline = file.path(paths$diagnostics, "baseline_state_draws.csv"),
  variable_map = file.path(paths$diagnostics, "variable_map.csv"),
  xlsx_lineage = file.path(paths$diagnostics, "xlsx_cache_lineage.csv"),
  summary = file.path(paths$diagnostics, "audit_summary.txt")
)
for (name in names(staged)) replace_file(staged[[name]], destinations[[name]])
message(paste(summary_lines, collapse = "\n"))
})()
