#!/usr/bin/env Rscript
# analyse_artemis_stress.R
#
# Analysis for the Artemis surface stress test bug hunt (2026-08-02):
# a stale-pointer aliasing bug between the block-subsystem's devblock
# cache and the VM-level window cache in src/word_source/block_words.c,
# found by the first completed run of ART-STRESS-TEST, fixed, and
# re-verified on all three architectures.
#
# Reads experiments/artemis_stress/runs/*.csv, writes SVG charts to
# analysis/charts/ and PNG copies to analysis/img/.

suppressMessages({
  library(ggplot2)
  library(dplyr)
  library(tidyr)
  library(scales)
  library(svglite)
})

run_dir    <- file.path("experiments", "artemis_stress", "runs")
chart_dir  <- file.path("experiments", "artemis_stress", "analysis", "charts")
img_dir    <- file.path("experiments", "artemis_stress", "analysis", "img")
dir.create(chart_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(img_dir,   recursive = TRUE, showWarnings = FALSE)

## ---- palette (matches experiments/bare_metal house style) ----
sfblue   <- "#1F77B4"
sforange <- "#FF7F0E"
sfgreen  <- "#2CA02C"
sfred    <- "#D62728"
sfgray   <- "#505050"

theme_sf <- function() {
  theme_minimal(base_size = 12) +
    theme(
      panel.grid.minor = element_blank(),
      plot.title       = element_text(face = "bold"),
      plot.subtitle    = element_text(color = sfgray),
      legend.position  = "bottom"
    )
}

save_chart <- function(plot, name, width = 7, height = 4.2) {
  svg_path <- file.path(chart_dir, paste0(name, ".svg"))
  png_path <- file.path(img_dir,   paste0(name, ".png"))
  ggsave(svg_path, plot, width = width, height = height, device = svglite)
  ggsave(png_path, plot, width = width, height = height, dpi = 200)
  cat("  wrote", svg_path, "and", png_path, "\n")
}

## ---- load runs ----
runs <- list(
  amd64_buggy    = "amd64 (pre-fix)",
  amd64_fixed    = "amd64 (post-fix)",
  aarch64_fixed  = "aarch64 (post-fix)",
  riscv64_fixed  = "riscv64 (post-fix)"
)

load_run <- function(key, label) {
  df <- read.csv(file.path(run_dir, paste0(key, ".csv")))
  df$run   <- key
  df$label <- label
  df$status <- factor(ifelse(df$result == 1, "PASS", "FAIL"), levels = c("PASS", "FAIL"))
  df
}

all_runs <- bind_rows(Map(load_run, names(runs), runs))

## ---- Chart 1: pass rate before/after, all runs ----
summary_tbl <- all_runs %>%
  group_by(run, label) %>%
  summarise(n = n(), passed = sum(result), .groups = "drop") %>%
  mutate(pass_rate = passed / n,
         phase = ifelse(grepl("pre-fix", label), "pre-fix", "post-fix"))

p1 <- ggplot(summary_tbl, aes(x = reorder(label, pass_rate), y = pass_rate, fill = phase)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = sprintf("%d/%d", passed, n)), vjust = -0.5, size = 4) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1.08)) +
  scale_fill_manual(values = c("pre-fix" = sfred, "post-fix" = sfgreen)) +
  coord_flip() +
  labs(
    title    = "Artemis surface stress test: pass rate before/after fix",
    subtitle = "50 trials per run, unique random blocks across the 22998-block data pool",
    x = NULL, y = "Pass rate", fill = NULL
  ) +
  theme_sf()

save_chart(p1, "pass_rate_comparison", width = 8)

## ---- Chart 2: pass/fail sequence in write/verify order, buggy run ----
## A visual scan of this sequence first suggested serial clustering; the
## Wald-Wolfowitz runs test (H1, below) does not support that reading --
## observed runs (26) exactly match the count expected under a random
## ordering. The chart is kept because it is still the raw sequence the
## tests are run against, not because it shows a confirmed pattern.
buggy <- all_runs %>% filter(run == "amd64_buggy") %>% arrange(trial)

p2 <- ggplot(buggy, aes(x = trial, y = 1, fill = status)) +
  geom_tile(color = "white", linewidth = 0.6, height = 1) +
  scale_fill_manual(values = c(PASS = sfgreen, FAIL = sfred)) +
  scale_x_continuous(breaks = seq(0, 49, 5)) +
  labs(
    title    = "Pass/fail sequence in write order, amd64 pre-fix run",
    subtitle = "25 of 50 trials failed; a formal runs test (H1) finds no\nsignificant serial clustering in this sequence -- see Section 4",
    x = "Trial index (= write/verify order)", y = NULL, fill = NULL
  ) +
  theme_sf() +
  theme(axis.text.y = element_blank(), panel.grid = element_blank())

save_chart(p2, "failure_clustering", height = 2.6)

## ---- Chart 3: surface coverage map, buggy vs one fixed run ----
cov <- all_runs %>% filter(run %in% c("amd64_buggy", "amd64_fixed"))

p3 <- ggplot(cov, aes(x = lbn, y = trial, color = status)) +
  geom_point(size = 2.2, alpha = 0.85) +
  scale_color_manual(values = c(PASS = sfgreen, FAIL = sfred)) +
  facet_wrap(~label, ncol = 1) +
  labs(
    title    = "Surface coverage: scattered across the full data pool",
    subtitle = "Each point is one trial's block. A Wilcoxon test (H2) finds LBN and\noutcome are not independent, but no sub-region is visually implicated",
    x = "Logical block number (data pool: 3076-26073)", y = "Trial index", color = NULL
  ) +
  theme_sf()

save_chart(p3, "surface_coverage_map", height = 5.5)

## ---- summary table for the report ----
write.csv(summary_tbl, file.path("experiments", "artemis_stress", "analysis", "summary.csv"),
          row.names = FALSE)

cat("\nDone. Summary:\n")
print(as.data.frame(summary_tbl))

## ═══════════════════════════════════════════════════════════════════════
## Hypothesis tests
##
## Three questions, each with an explicit null hypothesis, tested against
## the actual recorded trial data rather than asserted from eyeballing
## charts. Results are written to tables/ for the LaTeX report to quote
## verbatim -- no numbers in the report are hand-typed.
## ═══════════════════════════════════════════════════════════════════════

table_dir <- file.path("experiments", "artemis_stress", "analysis", "tables")
dir.create(table_dir, recursive = TRUE, showWarnings = FALSE)

sink_lines <- c()
say <- function(...) {
  line <- paste0(...)
  sink_lines <<- c(sink_lines, line)
  cat(line, "\n")
}

say("Artemis surface stress test -- hypothesis tests")
say("=================================================")
say("")

## ---- H1: are pre-fix failures independently distributed across trials, ----
## ---- or do they show serial clustering (a runs test)?                  ----
## H0: the pass/fail sequence is a random ordering (no serial correlation).
## Manual Wald-Wolfowitz runs test (base R only, no extra package needed).
runs_test <- function(seq01) {
  n1 <- sum(seq01 == 1); n2 <- sum(seq01 == 0); n <- n1 + n2
  r  <- 1 + sum(diff(seq01) != 0)
  mu <- 1 + 2 * n1 * n2 / n
  var_r <- (2 * n1 * n2 * (2 * n1 * n2 - n)) / (n^2 * (n - 1))
  z <- (r - mu) / sqrt(var_r)
  p <- 2 * pnorm(-abs(z))
  list(n1 = n1, n2 = n2, runs = r, expected = mu, z = z, p = p)
}

buggy_seq <- buggy$result[order(buggy$trial)]
rt <- runs_test(buggy_seq)

say("H1 -- pre-fix failure clustering (Wald-Wolfowitz runs test)")
say("  H0: pass/fail outcomes are independently ordered (no clustering).")
say(sprintf("  Observed runs = %d (expected under H0 = %.2f, sd = %.2f)",
            rt$runs, rt$expected, sqrt((2*rt$n1*rt$n2*(2*rt$n1*rt$n2-(rt$n1+rt$n2))) /
              ((rt$n1+rt$n2)^2*(rt$n1+rt$n2-1)))))
say(sprintf("  z = %.3f, p = %.2e", rt$z, rt$p))
say(sprintf("  -> %s H0 at alpha=0.05: failures %s serially clustered.",
            ifelse(rt$p < 0.05, "REJECT", "fail to reject"),
            ifelse(rt$p < 0.05, "ARE", "are not")))
say("")

## ---- H1b: is outcome associated with trial/write order directly       ----
## ---- (not via the runs test, which only detects consecutive-value     ----
## ---- clustering and can miss a more general order-dependent trend)?   ----
## H0: trial index is drawn from the same distribution for PASS and FAIL.
wt_order <- wilcox.test(trial ~ status, data = buggy)

say("H1b -- pre-fix outcome vs. trial/write order directly (Wilcoxon)")
say("  H0: trial index is drawn from the same distribution for PASS/FAIL")
say("      (complements H1: the runs test only catches consecutive-value")
say("      clustering, not a general rank-order trend).")
say(sprintf("  W = %.1f, p = %.3f", wt_order$statistic, wt_order$p.value))
say(sprintf("  -> %s H0 at alpha=0.05: outcome %s associated with write order.",
            ifelse(wt_order$p.value < 0.05, "REJECT", "fail to reject"),
            ifelse(wt_order$p.value < 0.05, "IS", "is not")))
say("")

## ---- H2: is failure status independent of disk location (logical      ----
## ---- block number)?                                                    ----
## H0: LBN is drawn from the same distribution for PASS and FAIL trials.
## Wilcoxon rank-sum test (non-parametric, no normality assumption on LBN).
wt <- wilcox.test(lbn ~ status, data = buggy)

say("H2 -- pre-fix failure vs. disk location (Wilcoxon rank-sum test)")
say("  H0: logical block number is drawn from the same distribution")
say("      for PASS and FAIL trials (failure does not depend on location).")
say(sprintf("  W = %.1f, p = %.3f", wt$statistic, wt$p.value))
say(sprintf("  -> %s H0 at alpha=0.05: failure %s associated with LBN.",
            ifelse(wt$p.value < 0.05, "REJECT", "fail to reject"),
            ifelse(wt$p.value < 0.05, "IS", "is not")))
say("")

## Secondary check: sibling position within devblock (0/1/2).
buggy$sib_pos <- (buggy$lbn - 3076) %% 3
pos_tbl <- table(buggy$sib_pos, buggy$status)
ct <- suppressWarnings(chisq.test(pos_tbl))
say("  Secondary check -- sibling position (0/1/2) vs. outcome (chi-squared)")
say(sprintf("  chi-sq = %.3f, df = %d, p = %.3f",
            ct$statistic, ct$parameter, ct$p.value))
say(sprintf("  -> %s H0: outcome %s associated with sibling position.",
            ifelse(ct$p.value < 0.05, "REJECT", "fail to reject"),
            ifelse(ct$p.value < 0.05, "IS", "is not")))
say("")

## ---- H3: did the fix change the pass rate, or is 25/50 vs 150/150     ----
## ---- explainable by chance?                                            ----
## H0: pre-fix and post-fix runs share the same true pass probability.
post_passed <- sum(all_runs$run != "amd64_buggy" & all_runs$result == 1)
post_n      <- sum(all_runs$run != "amd64_buggy")
pre_passed  <- sum(all_runs$run == "amd64_buggy" & all_runs$result == 1)
pre_n       <- sum(all_runs$run == "amd64_buggy")

pt <- prop.test(c(pre_passed, post_passed), c(pre_n, post_n), correct = TRUE)

say("H3 -- did the fix change the pass rate (two-proportion test)")
say("  H0: pre-fix and post-fix trials share the same true pass rate.")
say(sprintf("  Pre-fix: %d/%d = %.1f%%   Post-fix (pooled, 3 arches): %d/%d = %.1f%%",
            pre_passed, pre_n, 100*pre_passed/pre_n,
            post_passed, post_n, 100*post_passed/post_n))
say(sprintf("  chi-sq = %.2f, df = %d, p = %.2e", pt$statistic, 1, pt$p.value))
say(sprintf("  95%% CI on the difference in pass rate: [%.3f, %.3f]",
            pt$conf.int[1], pt$conf.int[2]))
say(sprintf("  -> %s H0 at alpha=0.05: the fix %s the pass rate.",
            ifelse(pt$p.value < 0.05, "REJECT", "fail to reject"),
            ifelse(pt$p.value < 0.05, "DID CHANGE", "did not measurably change")))

writeLines(sink_lines, file.path(table_dir, "hypothesis_tests.txt"))

## Machine-readable version for \newcommand substitution in the report.
hyp_df <- data.frame(
  key   = c("h1_runs_obs","h1_runs_exp","h1_z","h1_p",
            "h1b_W","h1b_p",
            "h2_W","h2_p","h2b_chisq","h2b_df","h2b_p",
            "h3_pre_n","h3_pre_pass","h3_post_n","h3_post_pass",
            "h3_chisq","h3_p","h3_ci_lo","h3_ci_hi"),
  value = c(rt$runs, rt$expected, rt$z, rt$p,
            unname(wt_order$statistic), wt_order$p.value,
            unname(wt$statistic), wt$p.value, unname(ct$statistic), unname(ct$parameter), ct$p.value,
            pre_n, pre_passed, post_n, post_passed,
            unname(pt$statistic), pt$p.value, pt$conf.int[1], pt$conf.int[2])
)
write.csv(hyp_df, file.path(table_dir, "hypothesis_tests.csv"), row.names = FALSE)
cat("\nWrote", file.path(table_dir, "hypothesis_tests.txt"), "and .csv\n")