Files

512 lines
23 KiB
R
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
#!/usr/bin/env Rscript
# analyse_multirep.R
# StarForth LithosAnanke — Multi-Replication DoE Analysis
#
# Pools all 9 outer DoE cells (3 reps × 3 arches) and tests:
# 1. Replication main effect on attractor metrics (should be nil)
# 2. Seed → shuffle independence (independent permutations = independent reps)
# 3. Inter-rep consistency of execution rate invariance
# 4. Per-rep Kruskal-Wallis to verify arch divergence in timer quality is
# reproducible across seeds
#
# Rep 1: seed 12345 (amd64, aarch64, riscv64)
# Rep 2: seed 67890 (riscv64, aarch64, amd64)
# Rep 3: seed 13579 (aarch64, amd64, riscv64)
suppressPackageStartupMessages({
library(ggplot2)
library(svglite)
library(dplyr)
library(tidyr)
library(scales)
library(patchwork)
library(viridis)
})
SCRIPT_DIR <- tryCatch(
dirname(normalizePath(sys.frames()[[1]]$ofile)),
error = function(e) getwd()
)
BASE_DIR <- normalizePath(file.path(SCRIPT_DIR, ".."))
RUNS_DIR <- file.path(BASE_DIR, "runs")
OUT_CHARTS <- file.path(SCRIPT_DIR, "charts")
OUT_TABLES <- file.path(SCRIPT_DIR, "tables")
dir.create(OUT_CHARTS, showWarnings = FALSE, recursive = TRUE)
dir.create(OUT_TABLES, showWarnings = FALSE, recursive = TRUE)
cat("══════════════════════════════════════════════════════════════════\n")
cat(" StarForth LithosAnanke — Multi-Replication DoE Analysis\n")
cat(" 9 outer cells (3 reps × 3 arches), ~783k rows each rep\n")
cat("══════════════════════════════════════════════════════════════════\n\n")
Q48 <- 65536.0
# ── manifest of all 9 canonical runs ─────────────────────────────────────────
manifest <- list(
list(rep=1, seed=12345, arch="amd64", file="doe-amd64-20260612-110212.csv"),
list(rep=1, seed=12345, arch="aarch64", file="doe-aarch64-20260612-132412.csv"),
list(rep=1, seed=12345, arch="riscv64", file="doe-riscv64-20260612-135612.csv"),
list(rep=2, seed=67890, arch="riscv64", file="doe-riscv64-20260612-141129.csv"),
list(rep=2, seed=67890, arch="aarch64", file="doe-aarch64-20260612-141722.csv"),
list(rep=2, seed=67890, arch="amd64", file="doe-amd64-20260612-142323.csv"),
list(rep=3, seed=13579, arch="aarch64", file="doe-aarch64-20260612-160132.csv"),
list(rep=3, seed=13579, arch="amd64", file="doe-amd64-20260612-160656.csv"),
list(rep=3, seed=13579, arch="riscv64", file="doe-riscv64-20260612-174414.csv")
)
load_run <- function(m) {
path <- file.path(RUNS_DIR, m$file)
cat(sprintf(" Rep%d seed=%-6d %-8s ... ", m$rep, m$seed, m$arch))
df <- read.csv(path, stringsAsFactors = FALSE)
cat(sprintf("%d rows\n", nrow(df)))
df$rep <- m$rep
df$seed <- m$seed
df$arch <- m$arch
df <- df[!is.na(df$tick_interval_ns) & df$tick_interval_ns > 0, ]
df$heat <- df$avg_word_heat_q48 / Q48
df$variance <- df$variance_q48 / Q48
df$time_trust <- df$time_trust_q48 / Q48
df$elapsed_s <- df$elapsed_ns / 1e9
df
}
all_runs <- lapply(manifest, load_run)
df_all <- bind_rows(all_runs)
df_all$arch <- factor(df_all$arch, levels = c("amd64","aarch64","riscv64"))
df_all$rep <- factor(df_all$rep)
cat(sprintf("\n Total rows across all 9 cells: %s\n\n",
formatC(nrow(df_all), format="d", big.mark=",")))
# ── 1. Per-cell summary ───────────────────────────────────────────────────────
cat("Computing per-cell summary (rep × arch)...\n")
cell_tbl <- df_all %>%
group_by(rep, seed, arch) %>%
summarise(
n_ticks = n(),
mean_exec_delta = mean(word_executions_delta, na.rm=TRUE),
sd_exec_delta = sd(word_executions_delta, na.rm=TRUE),
mean_heat = mean(heat, na.rm=TRUE),
max_heat = max(heat, na.rm=TRUE),
mean_trust = mean(time_trust, na.rm=TRUE),
mean_variance = mean(variance, na.rm=TRUE),
mean_win = mean(window_width, na.rm=TRUE),
.groups = "drop"
)
write.csv(cell_tbl, file.path(OUT_TABLES, "cell_summary.csv"), row.names=FALSE)
cat(" Saved: cell_summary.csv\n")
print(as.data.frame(cell_tbl))
cat("\n")
# ── 2. Replication main effect on exec_delta (should be F≈0, p≈1) ───────────
cat("ANOVA: exec_delta ~ rep (should show nil replication effect)...\n")
fit_rep <- aov(mean_exec_delta ~ rep, data = cell_tbl)
print(summary(fit_rep))
cat("\n")
# ── 3. Per-rep Kruskal-Wallis: timer quality divergence reproducible? ─────────
cat("Per-rep Kruskal-Wallis: time_trust ~ arch (divergent region only)\n")
cat("──────────────────────────────────────────────────────────────────\n")
kw_results <- list()
for (r in 1:3) {
df_r <- df_all %>%
filter(rep == r) %>%
group_by(arch) %>%
filter(row_number() > 26251) %>%
ungroup()
kw <- kruskal.test(time_trust ~ arch, data = df_r)
kw_results[[r]] <- data.frame(
rep = r,
seed = unique(df_all$seed[df_all$rep == r]),
H = round(kw$statistic, 1),
df = kw$parameter,
p_value = kw$p.value,
decision = if (kw$p.value < 2.2e-16) "REJECT H0 (arch divergence confirmed)" else "fail to reject"
)
cat(sprintf(" Rep %d (seed %5d): H=%.1f df=%d p=%.2e → %s\n",
r, kw_results[[r]]$seed, kw$statistic, kw$parameter, kw$p.value,
kw_results[[r]]$decision))
}
kw_tbl <- bind_rows(kw_results)
write.csv(kw_tbl, file.path(OUT_TABLES, "kw_per_rep.csv"), row.names=FALSE)
cat(" Saved: kw_per_rep.csv\n\n")
# ── 4. Exec rate: stable across reps? ────────────────────────────────────────
cat("Exec rate stability across reps (per arch):\n")
exec_stability <- cell_tbl %>%
group_by(arch) %>%
summarise(
rep_mean_of_means = mean(mean_exec_delta),
rep_sd_of_means = sd(mean_exec_delta),
cv_across_reps = sd(mean_exec_delta) / mean(mean_exec_delta) * 100,
.groups = "drop"
)
print(as.data.frame(exec_stability))
write.csv(exec_stability, file.path(OUT_TABLES, "exec_stability.csv"), row.names=FALSE)
cat(" Saved: exec_stability.csv\n\n")
# ── colour themes ─────────────────────────────────────────────────────────────
arch_colours <- c(amd64 = "#E07B39",
aarch64 = "#4A90D9",
riscv64 = "#50C878")
rep_colours <- c("1" = "#E07B39", "2" = "#4A90D9", "3" = "#50C878")
theme_light_report <- function() {
theme_minimal(base_size = 11) %+replace%
theme(
panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "grey90"),
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = 12),
plot.subtitle = element_text(colour = "grey40", size = 9),
legend.position = "bottom"
)
}
save_svg <- function(plot, name, w = 12, h = 7) {
path <- file.path(OUT_CHARTS, paste0(name, ".svg"))
svglite(path, width = w, height = h)
print(plot)
dev.off()
cat(sprintf(" Saved: %s.svg\n", name))
invisible(path)
}
cat("Generating multi-rep charts...\n\n")
# ── CHART MR-1: Exec rate per cell (3×3 grid, rep×arch) ──────────────────────
cat("[MR-1] Per-cell execution rate comparison...\n")
p_exec_cells <- ggplot(cell_tbl,
aes(x = arch, y = mean_exec_delta, colour = rep, shape = rep)) +
geom_point(size = 4, alpha = 0.85) +
geom_errorbar(aes(ymin = mean_exec_delta - sd_exec_delta,
ymax = mean_exec_delta + sd_exec_delta),
width = 0.15, alpha = 0.5) +
scale_colour_manual(values = rep_colours, name = "Replication") +
scale_shape_manual(values = c("1"=16,"2"=17,"3"=15), name = "Replication") +
scale_y_continuous(limits = c(253, 259)) +
labs(
title = "Execution Rate per Cell (mean ± SD)",
subtitle = "3 replications × 3 architectures — rate invariant across ISA and seed",
x = "Architecture",
y = "Words per 10 µs tick"
) +
theme_light_report()
save_svg(p_exec_cells, "multirep_exec_rate", w = 8, h = 6)
# ── CHART MR-2: Heat mean per cell ────────────────────────────────────────────
cat("[MR-2] Per-cell heat comparison...\n")
p_heat_cells <- ggplot(cell_tbl,
aes(x = rep, y = mean_heat, fill = arch, colour = arch)) +
geom_col(position = "dodge", alpha = 0.8, width = 0.6) +
scale_fill_manual(values = arch_colours) +
scale_colour_manual(values = arch_colours) +
labs(
title = "Mean Word Heat by Replication × Architecture",
subtitle = "amd64 heat ~8× higher due to APIC residual accumulation; consistent across reps",
x = "Replication",
y = "Mean avg word heat (Q48.16 decoded)"
) +
theme_light_report()
save_svg(p_heat_cells, "multirep_heat", w = 8, h = 6)
# ── CHART MR-3: Trust mean per cell ───────────────────────────────────────────
cat("[MR-3] Per-cell time-trust comparison...\n")
p_trust_cells <- ggplot(cell_tbl,
aes(x = rep, y = mean_trust, fill = arch, colour = arch)) +
geom_col(position = "dodge", alpha = 0.8, width = 0.6) +
scale_fill_manual(values = arch_colours) +
scale_colour_manual(values = arch_colours) +
coord_cartesian(ylim = c(0.99, 1.001)) +
labs(
title = "Mean Time-Trust by Replication × Architecture",
subtitle = "amd64 trust ~0.993; aarch64/riscv64 = 1.000 exactly — reproducible across all 3 seeds",
x = "Replication",
y = "Mean time_trust (Q48.16 decoded)"
) +
theme_light_report()
save_svg(p_trust_cells, "multirep_trust", w = 8, h = 6)
# ── CHART MR-4: KW H-statistic across reps (robustness check) ────────────────
cat("[MR-4] KW H-statistic reproducibility across reps...\n")
p_kw <- ggplot(kw_tbl, aes(x = factor(rep), y = H, fill = factor(rep))) +
geom_col(alpha = 0.85, width = 0.5) +
geom_text(aes(label = formatC(H, format="f", digits=0)),
vjust = -0.3, size = 4, fontface = "bold") +
scale_fill_manual(values = rep_colours, guide = "none") +
scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.15))) +
labs(
title = "Kruskal-Wallis H-Statistic: Timer Quality Divergence per Replication",
subtitle = "H ≈ 668,000 in all three reps — arch divergence is fully reproducible across seeds",
x = "Replication",
y = "K-W H statistic"
) +
theme_light_report() +
theme(legend.position = "none")
save_svg(p_kw, "multirep_kw_h", w = 7, h = 6)
# ── CHART MR-5: Window Poincaré map — overlay all 3 reps by arch ─────────────
cat("[MR-5] Multi-rep Poincaré map overlay (all 9 cells)...\n")
df_poin <- df_all %>%
mutate(
win_next = ave(window_width, interaction(rep, arch),
FUN = function(x) c(x[-1], NA_real_))
) %>%
filter(!is.na(win_next), window_width > 0, win_next > 0)
p_poin_mr <- ggplot(df_poin,
aes(x = log2(window_width), y = log2(win_next), colour = rep)) +
stat_bin_2d(bins = 40, aes(fill = after_stat(count)), alpha = 0.7) +
scale_fill_viridis_c(option = "plasma", trans = "sqrt",
name = "Density", labels = scales::comma) +
scale_colour_manual(values = rep_colours, guide = "none") +
scale_x_continuous(breaks = log2(c(256,512,1024,2048,4096)),
labels = c("256","512","1024","2048","4096")) +
scale_y_continuous(breaks = log2(c(256,512,1024,2048,4096)),
labels = c("256","512","1024","2048","4096")) +
geom_abline(slope = 1, intercept = 0, colour = "grey75",
linewidth = 0.4, linetype = "dashed") +
facet_wrap(~ arch, ncol = 3) +
labs(
title = "Multi-Replication Window Poincaré Map (Reps 13 pooled)",
subtitle = "3 seeds × 3 arches — attractor orbit structure identical across all replications",
x = expression(W[n]~"(log"[2]*" scale)"),
y = expression(W[n+1]~"(log"[2]*" scale)")
) +
theme_light_report() +
theme(legend.position = "right")
save_svg(p_poin_mr, "multirep_poincare", w = 14, h = 6)
# ══════════════════════════════════════════════════════════════════════════════
# CHART MR-6+7: 3×3 Thermal Phase Portrait — "the attractor grid"
# Rows = replication (Rep 1/2/3), Columns = architecture
# Styled after L8 capsule mode portraits: bold facet labels, dark + light
# Shows: same attractor across ALL 9 outer DoE cells
# ══════════════════════════════════════════════════════════════════════════════
cat("[MR-6+7] 3×3 attractor grid — thermal portrait (light + dark)...\n")
df_th9 <- df_all %>%
mutate(
delta_heat = ave(heat, interaction(rep, arch),
FUN = function(x) c(diff(x), NA_real_))
) %>%
filter(!is.na(delta_heat)) %>%
group_by(arch, rep) %>%
filter(abs(delta_heat) < quantile(abs(delta_heat), 0.995, na.rm = TRUE)) %>%
ungroup() %>%
mutate(
arch_label = factor(
case_when(
arch == "amd64" ~ "x86-64 (amd64)",
arch == "aarch64" ~ "ARMv8-A (aarch64)",
arch == "riscv64" ~ "RISC-V (riscv64)"
),
levels = c("x86-64 (amd64)", "ARMv8-A (aarch64)", "RISC-V (riscv64)")
),
rep_label = factor(
case_when(
rep == "1" ~ "REP 1 · seed 12345",
rep == "2" ~ "REP 2 · seed 67890",
rep == "3" ~ "REP 3 · seed 13579"
),
levels = c("REP 1 · seed 12345",
"REP 2 · seed 67890",
"REP 3 · seed 13579")
)
)
make_attractor_grid <- function(dark = FALSE) {
bg <- if (dark) "#0d0d0d" else "white"
lo <- if (dark) "#0d0d0d" else "white"
mid <- if (dark) "#003366" else "#ffdd00"
hi <- if (dark) "#00e5ff" else "#cc2200"
vhigh <- if (dark) "white" else "#1a0000"
grid_col <- if (dark) "#1a1a2e" else "grey88"
txt_col <- if (dark) "#aaaaaa" else "grey25"
ttl_col <- if (dark) "white" else "black"
stl_col <- if (dark) "#666666" else "grey45"
ggplot(df_th9, aes(x = heat, y = delta_heat)) +
stat_bin_2d(bins = 120, aes(fill = after_stat(count))) +
scale_fill_gradientn(
colours = c(lo, mid, hi, vhigh),
values = scales::rescale(c(0, 0.05, 0.30, 1)),
name = "Density",
trans = "sqrt",
labels = scales::comma
) +
geom_hline(yintercept = 0,
colour = if (dark) "#2a3a4a" else "grey78",
linewidth = 0.4, linetype = "dashed") +
facet_grid(rep_label ~ arch_label) +
labs(
title = "Compudynamic Attractor — Thermal Phase Space (All 9 Cells)",
subtitle = expression(paste("Heat"[n], " vs ", Delta, "Heat"[n],
" · 3 replications × 3 architectures · attractor topology invariant")),
x = "Avg word heat (Q48.16 decoded)",
y = expression(Delta * "Heat per tick")
) +
theme_minimal(base_size = 10) %+replace%
theme(
panel.background = element_rect(fill = bg, colour = NA),
plot.background = element_rect(fill = bg, colour = NA),
panel.grid.major = element_line(colour = grid_col),
panel.grid.minor = element_blank(),
axis.text = element_text(colour = txt_col, size = 7),
axis.title = element_text(colour = txt_col, size = 9),
# Column labels (architecture) — bold, prominent like the reference
strip.text.x = element_text(colour = ttl_col, face = "bold",
size = 10, margin = margin(b = 4)),
# Row labels (replication)
strip.text.y = element_text(colour = ttl_col, face = "bold",
size = 8, angle = 270,
margin = margin(l = 4)),
strip.background = element_rect(fill = if(dark) "#1a1a2e" else "grey95",
colour = NA),
plot.title = element_text(colour = ttl_col, face = "bold",
size = 12, margin = margin(b = 4)),
plot.subtitle = element_text(colour = stl_col, size = 8,
margin = margin(b = 8)),
legend.position = "right",
legend.background = element_rect(fill = bg, colour = NA),
legend.text = element_text(colour = txt_col, size = 7),
legend.title = element_text(colour = txt_col, size = 8),
panel.spacing = unit(0.6, "lines")
)
}
save_svg(make_attractor_grid(FALSE), "attractor_grid_light", w = 14, h = 10)
save_svg(make_attractor_grid(TRUE), "attractor_grid_dark", w = 14, h = 10)
# ══════════════════════════════════════════════════════════════════════════════
# CHART MR-8+9: 3×3 Poincaré Map Grid — window attractor, all 9 cells
# ══════════════════════════════════════════════════════════════════════════════
cat("[MR-8+9] 3×3 attractor grid — Poincaré map (light + dark)...\n")
df_po9 <- df_all %>%
mutate(
win_next = ave(window_width, interaction(rep, arch),
FUN = function(x) c(x[-1], NA_real_))
) %>%
filter(!is.na(win_next), window_width > 0, win_next > 0) %>%
mutate(
arch_label = factor(
case_when(
arch == "amd64" ~ "x86-64 (amd64)",
arch == "aarch64" ~ "ARMv8-A (aarch64)",
arch == "riscv64" ~ "RISC-V (riscv64)"
),
levels = c("x86-64 (amd64)", "ARMv8-A (aarch64)", "RISC-V (riscv64)")
),
rep_label = factor(
case_when(
rep == "1" ~ "REP 1 · seed 12345",
rep == "2" ~ "REP 2 · seed 67890",
rep == "3" ~ "REP 3 · seed 13579"
),
levels = c("REP 1 · seed 12345",
"REP 2 · seed 67890",
"REP 3 · seed 13579")
)
)
make_poincare_grid <- function(dark = FALSE) {
bg <- if (dark) "#0d0d0d" else "white"
lo <- if (dark) "#0d0d0d" else "white"
mid <- if (dark) "#003366" else "#ffdd00"
hi <- if (dark) "#00e5ff" else "#cc2200"
vhigh <- if (dark) "white" else "#1a0000"
grid_col <- if (dark) "#1a1a2e" else "grey88"
txt_col <- if (dark) "#aaaaaa" else "grey25"
ttl_col <- if (dark) "white" else "black"
stl_col <- if (dark) "#666666" else "grey45"
ggplot(df_po9, aes(x = log2(window_width), y = log2(win_next))) +
stat_bin_2d(bins = 40, aes(fill = after_stat(count))) +
scale_fill_gradientn(
colours = c(lo, mid, hi, vhigh),
values = scales::rescale(c(0, 0.05, 0.30, 1)),
name = "Density",
trans = "sqrt",
labels = scales::comma
) +
geom_abline(slope = 1, intercept = 0,
colour = if (dark) "#2a3a4a" else "grey78",
linewidth = 0.4, linetype = "dashed") +
scale_x_continuous(
breaks = log2(c(256, 512, 1024, 2048, 4096)),
labels = c("256","512","1024","2048","4096")
) +
scale_y_continuous(
breaks = log2(c(256, 512, 1024, 2048, 4096)),
labels = c("256","512","1024","2048","4096")
) +
facet_grid(rep_label ~ arch_label) +
labs(
title = "Compudynamic Attractor — Window Poincaré Map (All 9 Cells)",
subtitle = expression(paste(W[n], " vs ", W[n+1],
" (log"[2], " scale) · breathing cycle 256 ↔ 4096 preserved across all cells")),
x = expression(W[n]~"(words, log"[2]*" scale)"),
y = expression(W[n+1]~"(words, log"[2]*" scale)")
) +
theme_minimal(base_size = 10) %+replace%
theme(
panel.background = element_rect(fill = bg, colour = NA),
plot.background = element_rect(fill = bg, colour = NA),
panel.grid.major = element_line(colour = grid_col),
panel.grid.minor = element_blank(),
axis.text = element_text(colour = txt_col, size = 6),
axis.title = element_text(colour = txt_col, size = 9),
strip.text.x = element_text(colour = ttl_col, face = "bold",
size = 10, margin = margin(b = 4)),
strip.text.y = element_text(colour = ttl_col, face = "bold",
size = 8, angle = 270,
margin = margin(l = 4)),
strip.background = element_rect(fill = if(dark) "#1a1a2e" else "grey95",
colour = NA),
plot.title = element_text(colour = ttl_col, face = "bold",
size = 12, margin = margin(b = 4)),
plot.subtitle = element_text(colour = stl_col, size = 8,
margin = margin(b = 8)),
legend.position = "right",
legend.background = element_rect(fill = bg, colour = NA),
legend.text = element_text(colour = txt_col, size = 7),
legend.title = element_text(colour = txt_col, size = 8),
panel.spacing = unit(0.6, "lines")
)
}
save_svg(make_poincare_grid(FALSE), "poincare_grid_light", w = 14, h = 10)
save_svg(make_poincare_grid(TRUE), "poincare_grid_dark", w = 14, h = 10)
# ── summary ───────────────────────────────────────────────────────────────────
svg_new <- list.files(OUT_CHARTS, pattern = "^multirep.*\\.svg$")
cat(sprintf("\n══════════════════════════════════════════════════════════════════\n"))
cat(sprintf(" Multi-rep analysis complete. %d new charts.\n", length(svg_new)))
cat(sprintf("══════════════════════════════════════════════════════════════════\n\n"))
cat("Key findings:\n")
cat(sprintf(" Exec rate CV across reps: amd64=%.4f%% / aarch64=%.4f%% / riscv64=%.4f%%\n",
exec_stability$cv_across_reps[exec_stability$arch=="amd64"],
exec_stability$cv_across_reps[exec_stability$arch=="aarch64"],
exec_stability$cv_across_reps[exec_stability$arch=="riscv64"]))
cat(sprintf(" KW H range: %.0f %.0f (all p < 2.2e-16)\n",
min(kw_tbl$H), max(kw_tbl$H)))
cat(" Conclusion: attractor metrics stable across seeds; timer divergence reproducible.\n\n")