#!/usr/bin/env Rscript
# analyse_bare_metal.R
# StarForth LithosAnanke — Multi-Architecture Bare-Metal DoE Analysis
#
# Key insight from data exploration:
#  - tick_interval_ns is constant at 10000 (QEMU deterministic 10µs timer).
#    HR-based phase portraits (as in l8_attractor_map) do not apply here.
#  - The interesting dynamics are in:
#      window_width  — adaptive rolling window BREATHES 256→512→1024→2048→4096→...
#      avg_word_heat_q48 — monotonically rising thermal accumulation
#      hot_word_count — convergent excited-state count
#      time_trust_q48 / variance_q48 — arch-specific timer quality metrics
#  - Rows 1~26251 are IDENTICAL across all three ISAs (deterministic boot/POST).
#  - Divergence starts at row ~26252 in arch-specific timer metrics.
#
# The centrepiece chart: window_width[n] vs window_width[n+1] (Poincaré map).
# The adaptive window orbits a fixed cycle through powers of 2 — this IS the
# compudynamics attractor for bare-metal execution.

suppressPackageStartupMessages({
  library(ggplot2)
  library(svglite)
  library(dplyr)
  library(tidyr)
  library(scales)
  library(patchwork)
  library(viridis)
})

# ── paths ─────────────────────────────────────────────────────────────────────
SCRIPT_DIR <- tryCatch(
  dirname(normalizePath(sys.frames()[[1]]$ofile)),
  error = function(e) getwd()
)
BASE_DIR   <- normalizePath(file.path(SCRIPT_DIR, ".."))
DATA_DIR   <- file.path(BASE_DIR, "latest")
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 — Bare-Metal DoE Analysis\n")
cat("══════════════════════════════════════════════════════════════════\n\n")

Q48 <- 65536.0

load_arch <- function(arch) {
  path <- file.path(DATA_DIR, paste0(arch, ".csv"))
  cat(sprintf("  Loading %s ... ", arch))
  df <- read.csv(path, stringsAsFactors = FALSE)
  cat(sprintf("%d rows\n", nrow(df)))
  df$arch <- arch

  # Guard: remove degenerate rows
  df <- df[!is.na(df$tick_interval_ns) & df$tick_interval_ns > 0, ]

  # Q48.16 decode
  df$heat       <- df$avg_word_heat_q48 / Q48
  df$variance   <- df$variance_q48      / Q48
  df$time_trust <- df$time_trust_q48    / Q48

  # Elapsed in seconds
  df$elapsed_s  <- df$elapsed_ns / 1e9

  # Lag-1 window step for Poincaré map
  df$win_next   <- c(df$window_width[-1], NA_real_)
  df$delta_win  <- c(diff(df$window_width), NA_real_)
  df$delta_heat <- c(diff(df$heat), NA_real_)

  df
}

archs    <- c("amd64", "aarch64", "riscv64")
all_data <- lapply(archs, load_arch)
names(all_data) <- archs
df_all <- bind_rows(all_data)
df_all$arch <- factor(df_all$arch, levels = archs)

cat(sprintf("\n  Total rows: %d\n\n", nrow(df_all)))

# ── summary table ─────────────────────────────────────────────────────────────
cat("Computing per-arch summary...\n")

summary_tbl <- df_all %>%
  group_by(arch) %>%
  summarise(
    n_ticks          = n(),
    elapsed_s        = max(elapsed_s, na.rm = TRUE),
    mean_exec_delta  = mean(word_executions_delta, na.rm = TRUE),
    sd_exec_delta    = sd(word_executions_delta,   na.rm = TRUE),
    mean_hot_words   = mean(hot_word_count,         na.rm = TRUE),
    mean_heat        = mean(heat,                   na.rm = TRUE),
    max_heat         = max(heat,                    na.rm = TRUE),
    mean_win         = mean(window_width,           na.rm = TRUE),
    mean_trust       = mean(time_trust,             na.rm = TRUE),
    mean_variance    = mean(variance,               na.rm = TRUE),
    .groups = "drop"
  )

write.csv(summary_tbl, file.path(OUT_TABLES, "arch_summary.csv"), row.names = FALSE)
cat("  Saved: arch_summary.csv\n")
print(as.data.frame(summary_tbl))

# Statistical tests on the divergent region (rows after ~26251)
# Use time_trust and variance which differ per arch.
# If run length == 26251 (3-rep standard sweep) df_tail will be empty;
# fall back to the last 10% of each arch's data.
df_tail <- df_all %>%
  group_by(arch) %>%
  filter(row_number() > 26251) %>%
  ungroup()

if (nrow(df_tail) == 0) {
  cat("  NOTE: No rows beyond 26251 — using last 10% per arch as divergent region\n")
  df_tail <- df_all %>%
    group_by(arch) %>%
    slice_tail(prop = 0.10) %>%
    ungroup()
}

n_groups_trust <- df_tail %>% pull(time_trust) %>% unique() %>% length()
n_groups_var   <- df_tail %>% pull(variance)   %>% unique() %>% length()

kw_note <- ""
if (n_groups_trust < 2) {
  kw_note <- "SKIP — all time_trust values identical across architectures (algorithmic variance = 0)"
  cat(sprintf("\nKruskal-Wallis: time_trust_q48 ~ architecture: %s\n", kw_note))
  kw_trust <- list(statistic = 0, parameter = 0, p.value = 1)
} else {
  cat("\nKruskal-Wallis: time_trust_q48 ~ architecture (divergent region)\n")
  kw_trust <- kruskal.test(time_trust ~ arch, data = df_tail)
  cat(sprintf("  H=%.4f  df=%d  p=%.4e\n", kw_trust$statistic, kw_trust$parameter, kw_trust$p.value))
}

if (n_groups_var < 2) {
  kw_var_note <- "SKIP — all variance values identical across architectures (algorithmic variance = 0)"
  cat(sprintf("Kruskal-Wallis: variance_q48 ~ architecture: %s\n\n", kw_var_note))
  kw_var <- list(statistic = 0, parameter = 0, p.value = 1)
} else {
  cat("Kruskal-Wallis: variance_q48 ~ architecture (divergent region)\n")
  kw_var <- kruskal.test(variance ~ arch, data = df_tail)
  cat(sprintf("  H=%.4f  df=%d  p=%.4e\n\n", kw_var$statistic, kw_var$parameter, kw_var$p.value))
}

sink(file.path(OUT_TABLES, "kw_results.txt"))
cat("Kruskal-Wallis tests (divergent region)\n")
cat("════════════════════════════════════════\n\n")
if (nzchar(kw_note)) {
  cat(sprintf("time_trust_q48 ~ architecture: %s\n", kw_note))
} else {
  cat("time_trust_q48 ~ architecture\n"); print(kw_trust)
}
cat("\n")
if (nzchar(kw_note)) {
  cat(sprintf("variance_q48 ~ architecture: %s\n", kw_note))
} else {
  cat("variance_q48 ~ architecture\n"); print(kw_var)
}
sink()
cat("  Saved: kw_results.txt\n\n")

# ── colour themes ─────────────────────────────────────────────────────────────
arch_colours <- c(amd64   = "#E07B39",
                  aarch64 = "#4A90D9",
                  riscv64 = "#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"
    )
}

theme_dark_report <- function() {
  theme_minimal(base_size = 11) %+replace%
    theme(
      panel.background = element_rect(fill = "#0d0d0d", colour = NA),
      plot.background  = element_rect(fill = "#0d0d0d", colour = NA),
      panel.grid.major = element_line(colour = "#1e1e1e"),
      panel.grid.minor = element_blank(),
      axis.text        = element_text(colour = "#aaaaaa"),
      axis.title       = element_text(colour = "#cccccc"),
      strip.text       = element_text(colour = "white", face = "bold"),
      plot.title       = element_text(colour = "white", face = "bold", size = 12),
      plot.subtitle    = element_text(colour = "#666666", size = 9),
      legend.text      = element_text(colour = "#aaaaaa"),
      legend.title     = element_text(colour = "#cccccc"),
      legend.background = element_rect(fill = "#0d0d0d", colour = NA),
      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 charts...\n\n")

# ══════════════════════════════════════════════════════════════════════════════
# CHART 1+2: Window Poincaré Map — THE CENTREPIECE
# window_width[n] vs window_width[n+1], density-coloured
# Shows the adaptive window's fixed-cycle attractor
# ══════════════════════════════════════════════════════════════════════════════
cat("[1+2] Window Poincaré map (centrepiece — light + dark)...\n")

df_poin <- df_all %>%
  filter(!is.na(win_next), window_width > 0, win_next > 0)

NBINS <- 60   # coarser bins since window_width has only 5 discrete levels

make_poincare <- function(dark = FALSE) {
  bg       <- if (dark) "#0d0d0d" else "white"
  lo       <- if (dark) "#0d0d0d"  else "white"
  mid      <- if (dark) "#003366"  else "#ffcc00"
  hi       <- if (dark) "#00e5ff"  else "#cc0000"
  vhigh    <- if (dark) "white"    else "#1a0000"
  grid_col <- if (dark) "#1a1a2e"  else "grey88"
  txt_col  <- if (dark) "#aaaaaa"  else "grey30"
  ttl_col  <- if (dark) "white"    else "black"

  # Log-scale window for better visual separation of powers-of-2
  ggplot(df_poin, aes(x = log2(window_width), y = log2(win_next))) +
    stat_bin_2d(bins = NBINS, aes(fill = after_stat(count))) +
    scale_fill_gradientn(
      colours = c(lo, mid, hi, vhigh),
      values  = scales::rescale(c(0, 0.05, 0.3, 1)),
      name    = "Density",
      trans   = "sqrt",
      labels  = scales::comma
    ) +
    geom_abline(slope = 1, intercept = 0,
                colour = if (dark) "#334455" else "grey75",
                linewidth = 0.5, 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_wrap(~ arch, ncol = 3) +
    labs(
      title    = "Adaptive Window Poincaré Map — Bare-Metal Multi-Architecture",
      subtitle = "window_width[n] vs window_width[n+1]; diagonal = identity; orbit = adaptive breathing cycle",
      x        = expression(W[n]~"(words, log"[2]*" scale)"),
      y        = expression(W[n+1]~"(words, log"[2]*" scale)")
    ) +
    theme_minimal(base_size = 11) %+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 = 8),
      axis.title        = element_text(colour = txt_col),
      strip.text        = element_text(colour = ttl_col, face = "bold"),
      plot.title        = element_text(colour = ttl_col, face = "bold", size = 13),
      plot.subtitle     = element_text(colour = if(dark) "#888888" else "grey40", size = 9),
      legend.position   = "right",
      legend.background = element_rect(fill = bg, colour = NA),
      legend.text       = element_text(colour = txt_col),
      legend.title      = element_text(colour = txt_col)
    )
}

save_svg(make_poincare(dark = FALSE), "poincare_light", w = 14, h = 6)
save_svg(make_poincare(dark = TRUE),  "poincare_dark",  w = 14, h = 6)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 3+4: Window width time series — shows the breathing oscillation
# ══════════════════════════════════════════════════════════════════════════════
cat("[3+4] Window width time series (light + dark)...\n")

df_win_ts <- df_all %>%
  group_by(arch) %>%
  slice(seq(1, n(), by = 3)) %>%
  ungroup()

make_window_ts <- function(dark = FALSE) {
  base <- ggplot(df_win_ts, aes(x = elapsed_s, y = window_width, colour = arch)) +
    geom_line(alpha = 0.65, linewidth = 0.35) +
    scale_colour_manual(values = arch_colours) +
    scale_y_log10(
      breaks = c(256, 512, 1024, 2048, 4096),
      labels = c("256","512","1024","2048","4096")
    ) +
    scale_x_continuous(labels = scales::comma) +
    facet_wrap(~ arch, ncol = 1, scales = "free_x") +
    labs(
      title    = "Adaptive Window Width Over Session",
      subtitle = "Log₂ scale; breathing cycle: 256→4096→256 repeating — the adaptive attractor",
      x        = "Elapsed time (s)",
      y        = "Window width (words)"
    )
  if (dark) base + theme_dark_report() + theme(legend.position = "none")
  else      base + theme_light_report() + theme(legend.position = "none")
}

save_svg(make_window_ts(FALSE), "window_ts_light", w = 12, h = 9)
save_svg(make_window_ts(TRUE),  "window_ts_dark",  w = 12, h = 9)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 5+6: Thermal accumulation — avg_word_heat time series
# ══════════════════════════════════════════════════════════════════════════════
cat("[5+6] Thermal accumulation time series (light + dark)...\n")

df_heat_ts <- df_all %>%
  group_by(arch) %>%
  slice(seq(1, n(), by = 5)) %>%
  ungroup()

make_heat_ts <- function(dark = FALSE) {
  base <- ggplot(df_heat_ts, aes(x = elapsed_s, y = heat, colour = arch)) +
    geom_line(alpha = 0.6, linewidth = 0.35) +
    scale_colour_manual(values = arch_colours) +
    scale_x_continuous(labels = scales::comma) +
    facet_wrap(~ arch, ncol = 1, scales = "free_x") +
    labs(
      title    = "Average Word Heat Accumulation Over Session",
      subtitle = "Q48.16 decoded; monotonic rise during DOE-WORK; arch-specific divergence in later phase",
      x        = "Elapsed time (s)",
      y        = "Avg word heat"
    )
  if (dark) base + theme_dark_report() + theme(legend.position = "none")
  else      base + theme_light_report() + theme(legend.position = "none")
}

save_svg(make_heat_ts(FALSE), "heat_ts_light", w = 12, h = 9)
save_svg(make_heat_ts(TRUE),  "heat_ts_dark",  w = 12, h = 9)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 7+8: Execution rate distribution (word_executions_delta)
# ══════════════════════════════════════════════════════════════════════════════
cat("[7+8] Execution rate distribution (light + dark)...\n")

df_exec <- df_all %>% filter(word_executions_delta > 0)

make_exec_hist <- function(dark = FALSE) {
  base <- ggplot(df_exec, aes(x = word_executions_delta, fill = arch)) +
    geom_histogram(binwidth = 8, alpha = 0.75, position = "identity") +
    scale_fill_manual(values = arch_colours) +
    scale_y_continuous(labels = scales::comma) +
    facet_wrap(~ arch, ncol = 1) +
    labs(
      title    = "FORTH Execution Rate Distribution",
      subtitle = "Words executed per 10 µs tick; ISA-independent rate convergence expected",
      x        = "Words per tick",
      y        = "Count"
    )
  if (dark) base + theme_dark_report() + theme(legend.position = "none")
  else      base + theme_light_report() + theme(legend.position = "none")
}

save_svg(make_exec_hist(FALSE), "exec_rate_light", w = 10, h = 9)
save_svg(make_exec_hist(TRUE),  "exec_rate_dark",  w = 10, h = 9)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 9+10: time_trust & variance by architecture (divergent region)
# ══════════════════════════════════════════════════════════════════════════════
cat("[9+10] Time-trust & variance density (light + dark)...\n")

df_trust <- df_tail %>% filter(time_trust < 1.01, time_trust > 0.9)

make_trust <- function(dark = FALSE) {
  p_trust <- ggplot(df_trust, aes(x = time_trust, fill = arch, colour = arch)) +
    geom_density(alpha = 0.4, linewidth = 0.6) +
    scale_fill_manual(values = arch_colours) +
    scale_colour_manual(values = arch_colours) +
    labs(title = "Time-Trust Distribution (divergent region)",
         subtitle = "Q48.16 decoded; timer quality diverges per ISA after boot phase",
         x = "Time trust", y = "Density")

  df_vq <- df_tail %>% filter(variance > 0, variance < quantile(variance, 0.99, na.rm=TRUE))
  p_var  <- ggplot(df_vq, aes(x = variance, fill = arch, colour = arch)) +
    geom_density(alpha = 0.4, linewidth = 0.6) +
    scale_fill_manual(values = arch_colours) +
    scale_colour_manual(values = arch_colours) +
    labs(title = "Timing Variance Distribution (divergent region)",
         subtitle = "Q48.16 decoded; lower variance = more stable adaptive timing",
         x = "Variance", y = "Density")

  if (dark) {
    p_trust <- p_trust + theme_dark_report()
    p_var   <- p_var   + theme_dark_report()
  } else {
    p_trust <- p_trust + theme_light_report()
    p_var   <- p_var   + theme_light_report()
  }
  p_trust / p_var
}

save_svg(make_trust(FALSE), "trust_variance_light", w = 10, h = 10)
save_svg(make_trust(TRUE),  "trust_variance_dark",  w = 10, h = 10)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 11+12: Hot word count evolution
# ══════════════════════════════════════════════════════════════════════════════
cat("[11+12] Hot word count evolution (light + dark)...\n")

df_hot_ts <- df_all %>%
  group_by(arch) %>%
  slice(seq(1, n(), by = 5)) %>%
  ungroup()

make_hot_ts <- function(dark = FALSE) {
  base <- ggplot(df_hot_ts, aes(x = elapsed_s, y = hot_word_count, colour = arch)) +
    geom_line(alpha = 0.55, linewidth = 0.35) +
    scale_colour_manual(values = arch_colours) +
    scale_x_continuous(labels = scales::comma) +
    facet_wrap(~ arch, ncol = 1, scales = "free_x") +
    labs(
      title    = "Hot Word Count Over Session",
      subtitle = "Words above heat threshold; converges then differentiates per ISA in DoE phase",
      x        = "Elapsed time (s)",
      y        = "Hot word count"
    )
  if (dark) base + theme_dark_report() + theme(legend.position = "none")
  else      base + theme_light_report() + theme(legend.position = "none")
}

save_svg(make_hot_ts(FALSE), "hot_words_ts_light", w = 12, h = 9)
save_svg(make_hot_ts(TRUE),  "hot_words_ts_dark",  w = 12, h = 9)

# ══════════════════════════════════════════════════════════════════════════════
# CHART 13+14: Execution phase portrait (heat[n] vs delta_heat[n])
# Thermal dynamics — analogous to HR vs ΔHR in l8_attractor_map
# ══════════════════════════════════════════════════════════════════════════════
cat("[13+14] Thermal phase portrait (light + dark)...\n")

df_thermal <- df_all %>%
  filter(!is.na(delta_heat)) %>%
  filter(abs(delta_heat) < quantile(abs(delta_heat), 0.995, na.rm=TRUE))

make_thermal_portrait <- function(dark = FALSE) {
  bg       <- if (dark) "#0d0d0d" else "white"
  lo       <- if (dark) "#0d0d0d"  else "white"
  mid      <- if (dark) "#003366"  else "#ffcc00"
  hi       <- if (dark) "#00e5ff"  else "#cc0000"
  vhigh    <- if (dark) "white"    else "#1a0000"
  grid_col <- if (dark) "#1a1a2e"  else "grey88"
  txt_col  <- if (dark) "#aaaaaa"  else "grey30"
  ttl_col  <- if (dark) "white"    else "black"

  ggplot(df_thermal, aes(x = heat, y = delta_heat)) +
    stat_bin_2d(bins = 200, aes(fill = after_stat(count))) +
    scale_fill_gradientn(
      colours = c(lo, mid, hi, vhigh),
      values  = scales::rescale(c(0, 0.05, 0.3, 1)),
      name    = "Density",
      trans   = "sqrt",
      labels  = scales::comma
    ) +
    geom_hline(yintercept = 0,
               colour = if (dark) "#334455" else "grey75",
               linewidth = 0.4, linetype = "dashed") +
    facet_wrap(~ arch, ncol = 3) +
    labs(
      title    = "Thermal Phase Portrait — Heat[n] vs ΔHeat[n]",
      subtitle = "Thermal analogue of the l8 HR phase portrait; attractor at ΔHeat≈0 as system reaches steady state",
      x        = "Avg word heat (Q48.16 decoded)",
      y        = "ΔHeat per tick"
    ) +
    theme_minimal(base_size = 11) %+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 = 8),
      axis.title        = element_text(colour = txt_col),
      strip.text        = element_text(colour = ttl_col, face = "bold"),
      plot.title        = element_text(colour = ttl_col, face = "bold", size = 13),
      plot.subtitle     = element_text(colour = if(dark) "#888888" else "grey40", size = 9),
      legend.position   = "right",
      legend.background = element_rect(fill = bg, colour = NA),
      legend.text       = element_text(colour = txt_col),
      legend.title      = element_text(colour = txt_col)
    )
}

save_svg(make_thermal_portrait(FALSE), "thermal_portrait_light", w = 14, h = 6)
save_svg(make_thermal_portrait(TRUE),  "thermal_portrait_dark",  w = 14, h = 6)

# ══════════════════════════════════════════════════════════════════════════════
# done
# ══════════════════════════════════════════════════════════════════════════════
svg_files <- list.files(OUT_CHARTS, pattern = "\\.svg$")
cat(sprintf("\n══════════════════════════════════════════════════════════════════\n"))
cat(sprintf("  Analysis complete. %d SVG charts in %s\n", length(svg_files), OUT_CHARTS))
cat(sprintf("══════════════════════════════════════════════════════════════════\n\n"))