#!/usr/bin/env Rscript
# analysis.R -- 3(architecture) x 9(identity) factorial analysis of the
# std79-doe heartbeat/physics telemetry captured via HB-ON/HB-OFF
# (see ../README.md's results-20260912-with-heartbeat-csv/ entry).
#
# Input: combined.csv (this directory) -- one row per heartbeat tick,
# tagged with which (arch, run_id, id_idx, id_label, rep) trial was
# active when that tick's CSV row was printed, built by
# correlate_doe.py from the three raw logs in
# ../results-20260912-with-heartbeat-csv/. Q48.16 fields are already
# decoded to plain floats (avg_word_heat, time_trust, variance,
# hera/hermes/artemis_heat) and jitter_bits to jitter_ns (IEEE754
# bit-pattern reinterpretation), both done in Python before this script
# runs, since R has no native uint64/bit-punning story worth using here.
#
# Output: ANALYSIS.md (this directory) + svg/*.svg (referenced from it).

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

set.seed(1)
# Run from this directory (contains combined.csv); the caller cd's here.

df <- read.csv("combined.csv", stringsAsFactors = FALSE)

id_levels <- c("zuse", "rajames", "00", "01", "02", "03", "04", "05", "06")
arch_levels <- c("amd64", "aarch64", "riscv64")
df$id_label <- factor(df$id_label, levels = id_levels)
df$arch <- factor(df$arch, levels = arch_levels)

# apic_ticks is a monotonic hardware counter, not directly comparable
# across trials/architectures as a raw value -- the per-tick DELTA
# (hardware ticks consumed per heartbeat tick) is the actual timing
# signal. Computed per (arch, run_id) so a trial boundary never
# contaminates the delta of the row after it.
df <- df %>%
  arrange(arch, run_id, tick_number) %>%
  group_by(arch, run_id) %>%
  mutate(apic_delta = apic_ticks - lag(apic_ticks)) %>%
  ungroup()

metrics <- c(
  "hot_word_count", "avg_word_heat", "window_width", "actual_window_size",
  "jitter_ns", "apic_delta", "time_trust", "variance", "vm_call_depth_max",
  "hera_heat", "hermes_heat", "artemis_heat"
)

metric_labels <- c(
  hot_word_count = "Hot word count",
  avg_word_heat = "Mean word execution heat",
  window_width = "Rolling window width",
  actual_window_size = "Actual analysis window size",
  jitter_ns = "Estimated timer jitter (ns)",
  apic_delta = "APIC ticks per heartbeat tick",
  time_trust = "TIME-TRUST",
  variance = "Timing variance",
  vm_call_depth_max = "Max VM call depth",
  hera_heat = "Hera fleet heat",
  hermes_heat = "Hermes fleet heat",
  artemis_heat = "Artemis fleet heat"
)

dir.create("svg", showWarnings = FALSE)

# ---- Per-cell (arch x id_label) summary table for one metric ----
cell_summary <- function(metric) {
  df %>%
    filter(!is.na(.data[[metric]])) %>%
    group_by(arch, id_label) %>%
    summarise(
      n = n(),
      mean = mean(.data[[metric]]),
      sd = sd(.data[[metric]]),
      cv_pct = ifelse(mean != 0, 100 * sd / abs(mean), NA_real_),
      .groups = "drop"
    )
}

# ---- Two-way ANOVA (metric ~ arch * id_label), safe on zero-variance ----
anova_table <- function(metric) {
  d <- df %>% filter(!is.na(.data[[metric]]))
  if (length(unique(d[[metric]])) <= 1) {
    return(NULL)
  }
  fit <- tryCatch(
    aov(as.formula(paste0(metric, " ~ arch * id_label")), data = d),
    error = function(e) NULL
  )
  if (is.null(fit)) return(NULL)
  s <- summary(fit)[[1]]
  s
}

# ---- Boxplot: metric by id_label, faceted by arch ----
make_plot <- function(metric) {
  d <- df %>% filter(!is.na(.data[[metric]]))
  if (length(unique(d[[metric]])) <= 1) return(NULL)
  p <- ggplot(d, aes(x = id_label, y = .data[[metric]], fill = id_label)) +
    geom_boxplot(outlier.size = 0.6, alpha = 0.85) +
    facet_wrap(~arch, ncol = 1) +
    labs(
      title = metric_labels[[metric]],
      x = "Identity", y = metric_labels[[metric]]
    ) +
    theme_bw(base_size = 11) +
    theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1))
  fname <- file.path("svg", paste0(metric, ".svg"))
  ggsave(fname, p, width = 7, height = 8, device = svglite)
  fname
}

# ---- Run everything, collect markdown fragments ----
md <- c()
md <- c(md, "# std79-doe 3×9 Factorial Analysis — Heartbeat/Physics Telemetry")
md <- c(md, "")
md <- c(md, sprintf("Generated %s. Source: `combined.csv` (%d rows), built by", format(Sys.time(), "%Y-%m-%d %H:%M"), nrow(df)))
md <- c(md, "`correlate_doe.py` from the three raw logs in")
md <- c(md, "`../results-20260912-with-heartbeat-csv/`. See that file's own docstring for the")
md <- c(md, "correlation methodology (the heartbeat tick's async CSV printer and the DoE trial")
md <- c(md, "loop's own console output share one serial line with no locking, so rows can land")
md <- c(md, "mid-token in the raw logs; the underlying FORTH execution and values are")
md <- c(md, "unaffected, only reconstructing which trial owns which row needed care).")
md <- c(md, "")
md <- c(md, "**Design:** 3 architectures (amd64, aarch64, riscv64) × 9 identities (zuse,")
md <- c(md, "rajames, 00-06) × 3 replicates, Fisher-Yates-shuffled run order, same seed")
md <- c(md, "(12345) on every architecture (FABRIC-3.md §XV). Each of the 27 trials per")
md <- c(md, "architecture runs the same 24-word FORTH-79 exerciser; every heartbeat tick during")
md <- c(md, "the whole `EXEC-STD79-DOE` run (bracketed by `HB-ON`/`HB-OFF`) emits one telemetry")
md <- c(md, "row (`doe_log.c`), correlated back to whichever trial was active when it printed.")
md <- c(md, "")
md <- c(md, "**Known limitation:** aarch64's raw log lost 1 of 27 run_id markers to interleaving")
md <- c(md, "beyond recovery (trial 12, identity `rajames` rep 1) -- its rows are folded into")
md <- c(md, "trial 11 (identity `rajames` rep 0) in this dataset, so aarch64's `rajames` cell")
md <- c(md, "for rep-sensitive metrics is not perfectly separable for those two replicates")
md <- c(md, "specifically. Every other cell on every architecture is unaffected.")
md <- c(md, "")

# ---- Pass 1: compute everything, so the "Key findings" summary (which
# needs every metric's ANOVA result) can be written before the detailed
# per-metric sections that follow it. ----
results <- list()
for (metric in metrics) {
  results[[metric]] <- list(
    cs = cell_summary(metric),
    at = anova_table(metric),
    plot_file = make_plot(metric)
  )
}

sig <- function(p) !is.na(p) && p < 0.05

arch_sig <- c(); id_sig <- c(); inter_sig <- c(); none_sig <- c()
for (metric in metrics) {
  at <- results[[metric]]$at
  if (is.null(at)) next
  p_arch <- if ("arch" %in% trimws(rownames(at))) at[trimws(rownames(at)) == "arch", "Pr(>F)"] else NA
  p_id <- if ("id_label" %in% trimws(rownames(at))) at[trimws(rownames(at)) == "id_label", "Pr(>F)"] else NA
  p_int <- if ("arch:id_label" %in% trimws(rownames(at))) at[trimws(rownames(at)) == "arch:id_label", "Pr(>F)"] else NA
  label <- metric_labels[[metric]]
  any_sig <- FALSE
  if (sig(p_arch)) { arch_sig <- c(arch_sig, label); any_sig <- TRUE }
  if (sig(p_id)) { id_sig <- c(id_sig, label); any_sig <- TRUE }
  if (sig(p_int)) { inter_sig <- c(inter_sig, label); any_sig <- TRUE }
  if (!any_sig) none_sig <- c(none_sig, label)
}

md <- c(md, "## Key findings (p < 0.05, two-way ANOVA)")
md <- c(md, "")
md <- c(md, sprintf("- **Architecture main effect:** %s",
  if (length(arch_sig)) paste(arch_sig, collapse = ", ") else "none"))
md <- c(md, sprintf("- **Identity main effect:** %s",
  if (length(id_sig)) paste(id_sig, collapse = ", ") else "none"))
md <- c(md, sprintf("- **Architecture × identity interaction:** %s",
  if (length(inter_sig)) paste(inter_sig, collapse = ", ") else "none"))
md <- c(md, sprintf("- **No significant effect of either factor:** %s",
  if (length(none_sig)) paste(none_sig, collapse = ", ") else "none"))
md <- c(md, "")
md <- c(md, "TIME-TRUST's architecture effect is a numerical artifact worth reading correctly,")
md <- c(md, "not a substantive result: it is essentially constant *within* each architecture")
md <- c(md, "(residual variance ~1e-31, i.e. floating-point noise) and differs *between*")
md <- c(md, "architectures — so the F statistic is enormous simply because the within-group")
md <- c(md, "denominator is near zero, not because TIME-TRUST is meaningfully more variable")
md <- c(md, "across architectures than the other metrics here. It is architecture-determined")
md <- c(md, "and identity-independent, which is itself the interesting part.")
md <- c(md, "")

for (metric in metrics) {
  cs <- results[[metric]]$cs
  at <- results[[metric]]$at
  plot_file <- results[[metric]]$plot_file

  md <- c(md, sprintf("## %s (`%s`)", metric_labels[[metric]], metric))
  md <- c(md, "")

  if (is.null(plot_file)) {
    md <- c(md, "_Constant across every trial and architecture — no variance, excluded from ANOVA/plot._")
    md <- c(md, "")
    next
  }

  md <- c(md, sprintf("![%s](%s)", metric_labels[[metric]], plot_file))
  md <- c(md, "")
  md <- c(md, "| Architecture | Identity | n | mean | sd | CV% |")
  md <- c(md, "|---|---|---:|---:|---:|---:|")
  for (i in seq_len(nrow(cs))) {
    r <- cs[i, ]
    md <- c(md, sprintf("| %s | %s | %d | %.4g | %.4g | %s |",
      r$arch, r$id_label, r$n, r$mean, r$sd,
      ifelse(is.na(r$cv_pct), "—", sprintf("%.2f", r$cv_pct))))
  }
  md <- c(md, "")

  if (!is.null(at)) {
    md <- c(md, "Two-way ANOVA (`metric ~ arch * id_label`):")
    md <- c(md, "")
    md <- c(md, "| Term | Df | Sum Sq | Mean Sq | F value | Pr(>F) |")
    md <- c(md, "|---|---:|---:|---:|---:|---:|")
    terms <- rownames(at)
    for (i in seq_len(nrow(at))) {
      row <- at[i, ]
      # %g-style: near-zero within-cell residual variance (e.g. a metric
      # that's essentially constant within an architecture, only varying
      # BETWEEN architectures -- TIME-TRUST does exactly this) can blow
      # the F statistic up to absurd magnitudes; format generally rather
      # than with a fixed decimal count so that renders sanely too.
      fval <- if (!is.na(row[["F value"]])) formatC(row[["F value"]], format = "g", digits = 4) else "—"
      pval <- if (!is.na(row[["Pr(>F)"]])) formatC(row[["Pr(>F)"]], format = "g", digits = 4) else "—"
      md <- c(md, sprintf("| %s | %d | %.4g | %.4g | %s | %s |",
        trimws(terms[i]), row[["Df"]], row[["Sum Sq"]], row[["Mean Sq"]], fval, pval))
    }
    md <- c(md, "")
  } else {
    md <- c(md, "_ANOVA not computed (insufficient variance)._")
    md <- c(md, "")
  }
}

writeLines(md, "ANALYSIS.md")
cat("Wrote ANALYSIS.md and svg/*.svg\n")