Add 3x9 factorial analysis of std79-doe heartbeat/physics telemetry
Correlates every HB-ON/HB-OFF heartbeat-tick CSV row (results-20260912- with-heartbeat-csv/*-doe-raw.log) back to which trial (run_id, id_idx, id_label, rep) was active when it printed, despite the async tick printer splicing rows mid-token -- including mid a DOE-RUN marker itself -- into the trial loop's own console output on the shared serial line. Pipeline (analysis-20260912/, see its own README.md): - correlate_doe.py: two-pass reconstruction per architecture (remove atomic CSV-row spans to rebuild the clean trial-output stream, map each removed row's offset back to the nearest preceding run_id marker); identity/rep looked up from a known-clean prior run's run_id mapping rather than re-parsed, since one aarch64 marker (trial 12, rajames rep 1) lost its id_idx digit to a zero-separator collision with an adjacent CSV field and is unrecoverable from that log alone -- its rows fold into trial 11 instead, documented as a known limitation. - combine.py: merges all three architectures into combined.csv (768 rows), decoding Q48.16 fields to floats and jitter_bits' IEEE754 bit pattern to real jitter_ns. - analysis.R: per-cell (architecture x identity) means/SD and two-way ANOVA for each of 12 telemetry metrics, one boxplot SVG per metric, written up as ANALYSIS.md. Key findings: identity significantly affects word-heat/window-sizing metrics (expected -- different identities execute different word sets), architecture significantly affects timing metrics (APIC ticks/tick, timing variance, fleet heat -- expected, different QEMU targets), zero architecture x identity interaction on any metric. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EXieurDfDSsDFdnSyusuWo
This commit is contained in:
Robert Allan James
co-authored by
Claude Sonnet 5
Claude Sonnet 5
parent
403a7639e1
commit
934be5a257
@@ -0,0 +1,243 @@
|
||||
#!/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("", 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")
|
||||
Reference in New Issue
Block a user