LithosAnanke/experiments/std79-doe/report-20260912/analysis.R

#!/usr/bin/env Rscript
# SPDX-License-Identifier: LicenseRef-Starship-1.0
# Copyright (c) 2023-2026 Robert A. James
# analysis.R -- std79-doe 3x9 factorial analysis, K-conservation centerpiece.
# Run from this directory (contains combined.csv). Produces charts/*.svg
# (light+dark pairs, matching experiments/bare_metal/analysis's own
# convention) and tables/*.csv for the LaTeX report in report/.

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

dir.create("charts", showWarnings = FALSE)
dir.create("tables", showWarnings = FALSE)

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)

df <- df %>%
  arrange(arch, run_id, tick_number) %>%
  group_by(arch, run_id) %>%
  mutate(apic_delta = apic_ticks - lag(apic_ticks)) %>%
  ungroup()

# ── shared theme, matches experiments/bare_metal/analysis/analyse_bare_metal.R ──
theme_light_report <- function() {
  theme_minimal(base_size = 11) %+replace%
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(colour = "grey88"),
      plot.title = element_text(face = "bold", size = 12),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}
theme_dark_report <- function() {
  theme_minimal(base_size = 11) %+replace%
    theme(
      plot.background = element_rect(fill = "#0d0d0d", colour = NA),
      panel.background = element_rect(fill = "#0d0d0d", colour = NA),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(colour = "#1a1a2e"),
      text = element_text(colour = "#aaaaaa"),
      axis.text = element_text(colour = "#aaaaaa"),
      plot.title = element_text(face = "bold", size = 12, colour = "white"),
      axis.text.x = element_text(angle = 45, hjust = 1, colour = "#aaaaaa"),
      legend.text = element_text(colour = "#aaaaaa"),
      legend.background = element_rect(fill = "#0d0d0d")
    )
}
save_svg <- function(p, name, w = 8, h = 6) {
  ggsave(file.path("charts", paste0(name, ".svg")), p, width = w, height = h, device = svglite)
}

# ── K / conservation centerpiece: fleet_k across every tick, every arch ──
cat(sprintf("fleet_k range: [%.10f, %.10f], distinct values: %d\n",
            min(df$fleet_k), max(df$fleet_k), length(unique(df$fleet_k))))
cat(sprintf("fleet_conserved: %d/%d rows TRUE\n", sum(df$fleet_conserved == 1), nrow(df)))

make_k_plot <- function(dark = FALSE) {
  bg <- if (dark) "#0d0d0d" else "white"
  line_col <- if (dark) "#00e5ff" else "#1f77b4"
  ref_col <- if (dark) "#888888" else "grey50"
  p <- ggplot(df, aes(x = tick_number, y = fleet_k, colour = arch)) +
    geom_hline(yintercept = 1.0, linetype = "dashed", colour = ref_col, linewidth = 0.4) +
    geom_point(size = 0.5, alpha = 0.6) +
    facet_wrap(~arch, ncol = 1, scales = "free_x") +
    scale_y_continuous(limits = c(0.9, 1.1)) +
    labs(title = "Fleet conservation invariant K over the whole campaign",
         subtitle = "K = sum(execution_heat_q48) over all live VMs; dashed line = Q48_ONE (perfect conservation)",
         x = "Heartbeat tick number", y = "K") +
    theme(legend.position = "none")
  if (dark) p <- p + theme_dark_report() + theme(legend.position = "none")
  else      p <- p + theme_light_report() + theme(legend.position = "none")
  p
}
save_svg(make_k_plot(FALSE), "fleet_k_light", w = 9, h = 8)
save_svg(make_k_plot(TRUE),  "fleet_k_dark",  w = 9, h = 8)

# ── per-cell (arch x id_label) summary + ANOVA for every other metric ──
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"
)

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")
}
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)
  summary(fit)[[1]]
}
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(legend.position = "none")
  p_light <- p + theme_light_report() + theme(legend.position = "none")
  p_dark  <- p + theme_dark_report() + theme(legend.position = "none")
  save_svg(p_light, paste0(metric, "_light"), w = 7, h = 8)
  save_svg(p_dark,  paste0(metric, "_dark"),  w = 7, h = 8)
  TRUE
}

anova_rows <- list()
cell_rows <- list()
for (metric in metrics) {
  cs <- cell_summary(metric)
  cs$metric <- metric
  cell_rows[[metric]] <- cs

  at <- anova_table(metric)
  if (!is.null(at)) {
    at_df <- as.data.frame(at)
    at_df$term <- trimws(rownames(at))
    at_df$metric <- metric
    anova_rows[[metric]] <- at_df
  }
  make_plot(metric)
}

write.csv(bind_rows(cell_rows), "tables/cell_summary.csv", row.names = FALSE)
write.csv(bind_rows(anova_rows), "tables/anova.csv", row.names = FALSE)

# ── overall dataset summary table (per architecture) ──
overall <- df %>%
  group_by(arch) %>%
  summarise(
    n_rows = n(),
    n_trials = n_distinct(run_id),
    elapsed_s = max(elapsed_ns) / 1e9,
    mean_hot_words = mean(hot_word_count),
    mean_fleet_k = mean(fleet_k),
    sd_fleet_k = sd(fleet_k),
    pct_conserved = 100 * mean(fleet_conserved == 1),
    .groups = "drop"
  )
write.csv(overall, "tables/overall_summary.csv", row.names = FALSE)

cat("\nDone. charts/, tables/ written.\n")
print(as.data.frame(overall))