#!/usr/bin/env Rscript
# generate_stadium_deepdive.R
# Generates the full per-cell / per-ISA / per-interaction-pair / per-factor
# deep-dive: charts (SVG, light+dark) + LaTeX section fragments, for
# stadium_relaunch_report.tex's expanded edition.
#
# Captain Bob's brief: analyze all 9 cells as a conglomerate Latin square
# (done in analyse_stadium_relaunch_fixed.R), THEN dive into each cell's
# own data with its own charts/tables, THEN cover ISA-internal patterns
# and every factor x factor interaction (including quadratic terms).
# This script produces the "dive into each cell" + "every interaction"
# layer as a generated document, not hand-authored.

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

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", "acl-rwt-20260820-fixed")
OUT_CHARTS <- file.path(SCRIPT_DIR, "charts")
OUT_SEC    <- file.path(SCRIPT_DIR, "report", "sections")
dir.create(OUT_CHARTS, showWarnings = FALSE, recursive = TRUE)
dir.create(OUT_SEC, showWarnings = FALSE, recursive = TRUE)

arch_colours <- c(amd64 = "#E07B39", aarch64 = "#4A90D9", riscv64 = "#50C878")
theme_light_sf <- function(base = 10) {
  theme_minimal(base_size = base) %+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 = base + 1),
    plot.subtitle = element_text(colour = "grey40", size = base - 2),
    legend.position = "bottom", legend.key.size = unit(0.4, "cm"))
}
save_svg <- function(plot, name, w = 9, h = 5.5) {
  path <- file.path(OUT_CHARTS, paste0(name, ".svg"))
  svglite(path, width = w, height = h); print(plot); dev.off()
  invisible(path)
}
texify <- function(x) gsub("_", "\\\\_", x)

col_names <- c("run_id","cfg","rep","ent_in","cv_in","tmp_in","stb_in","l8_mode",
               "win_div","infer_win","infer_dec_q","infer_var_q","early_exit",
               "bc_mean_q","bb_mean_q","fit_q")
archs <- c("amd64","aarch64","riscv64"); seeds <- c("12345","67890","13579")

load_cell <- function(arch, seed) {
  f <- file.path(RUNS_DIR, sprintf("%s-seed%s.csv", arch, seed))
  df <- read.csv(f, header = TRUE); names(df) <- col_names
  df$arch <- arch; df$seed <- seed
  df$ent_f <- factor(ifelse(df$ent_in > 0, "hi", "lo"))
  df$cv_f  <- factor(ifelse(df$cv_in  > 0, "hi", "lo"))
  df$tmp_f <- factor(ifelse(df$tmp_in > 0, "hi", "lo"))
  df$stb_f <- factor(ifelse(df$stb_in > 0, "hi", "lo"))
  df
}
all_data <- bind_rows(lapply(archs, function(a) bind_rows(lapply(seeds, function(s) load_cell(a, s)))))
all_data$arch <- factor(all_data$arch, levels = archs)
all_data$seed <- factor(all_data$seed, levels = seeds)
cat(sprintf("Loaded %d rows across 9 cells.\n", nrow(all_data)))

section_files <- c()

# ════════════════════════════════════════════════════════════════════════
# PART A -- PER-CELL DEEP DIVE (9 sections)
# ════════════════════════════════════════════════════════════════════════
cat("\n=== PART A: per-cell deep dive ===\n")
for (a in archs) for (s in seeds) {
  cell <- filter(all_data, arch == a, seed == s)
  tag <- sprintf("%s_seed%s", a, s)
  cat(sprintf("  Cell %s...\n", tag))

  # cfg-level summary (16 rows)
  cfg_summary <- cell %>% group_by(cfg) %>%
    summarise(n = n(), mean_dec = mean(infer_dec_q), sd_dec = sd(infer_dec_q),
              early_exit_rate = mean(early_exit), .groups = "drop") %>%
    arrange(cfg)

  # Chart: histogram + boxplot-by-cfg, side by side
  p1 <- ggplot(cell, aes(x = infer_dec_q)) +
    geom_histogram(bins = 18, fill = arch_colours[[a]], alpha = 0.85, colour = NA) +
    labs(title = "infer\\_dec\\_q distribution", x = "infer_dec_q", y = "count") +
    theme_light_sf(9)
  p2 <- ggplot(cell, aes(x = factor(cfg), y = infer_dec_q)) +
    geom_boxplot(fill = arch_colours[[a]], alpha = 0.75, outlier.size = 0.5) +
    labs(title = "by L8 config", x = "cfg", y = "infer_dec_q") +
    theme_light_sf(9)
  combined <- p1 + p2
  save_svg(combined, sprintf("cell_%s_detail_light", tag), w = 11, h = 4.2)

  chart_path <- sprintf("cell_%s_detail_light", tag)

  tbl <- paste0(
    "\\begin{tabular}{rrrrr}\n\\toprule\n",
    "cfg & n & mean infer\\_dec\\_q & sd & early-exit rate \\\\\n\\midrule\n",
    paste(sprintf("%d & %d & %.2f & %.2f & %.3f \\\\", cfg_summary$cfg, cfg_summary$n,
                   cfg_summary$mean_dec, cfg_summary$sd_dec, cfg_summary$early_exit_rate),
          collapse = "\n"),
    "\n\\bottomrule\n\\end{tabular}\n")

  overall_mean <- mean(cell$infer_dec_q); overall_sd <- sd(cell$infer_dec_q)
  overall_ee   <- mean(cell$early_exit)
  q <- quantile(cell$infer_dec_q, probs = c(0.25, 0.5, 0.75))

  sec <- sprintf("\\subsection{Cell: %s, seed %s}\n\\label{sec:cell_%s}\n\n",
                  a, s, tag)
  sec <- paste0(sec, sprintf(
    "480 rows captured, all 16 \\texttt{cfg} values represented exactly 30 times, zero errors. Overall \\texttt{infer\\_dec\\_q}: mean %.2f, sd %.2f, median %.1f (IQR %.1f--%.1f). Early-exit rate: %.3f.\n\n",
    overall_mean, overall_sd, q[2], q[1], q[3], overall_ee))
  sec <- paste0(sec, sprintf(
    "\\begin{figure}[h]\n\\centering\n\\includegraphics[width=\\textwidth]{%s.pdf}\n\\caption{Cell %s/seed %s: \\texttt{infer\\_dec\\_q} distribution (left) and per-config breakdown (right).}\n\\end{figure}\n\n",
    chart_path, texify(a), s))
  sec <- paste0(sec, "\\begin{table}[h]\n\\centering\n\\caption{Per-config summary, cell ", texify(a), "/seed ", s, ".}\n", tbl, "\\end{table}\n\n\\clearpage\n")

  # Rep-trajectory chart: infer_dec_q across the 480 sequential runs,
  # visualizing order-dependence (why seed matters, architecture doesn't).
  cell_ord <- cell %>% arrange(run_id)
  p3 <- ggplot(cell_ord, aes(x = run_id, y = infer_dec_q)) +
    geom_line(colour = arch_colours[[a]], alpha = 0.5, linewidth = 0.3) +
    geom_point(aes(colour = factor(cfg)), size = 0.9, show.legend = FALSE) +
    geom_smooth(se = FALSE, colour = "black", linewidth = 0.6, method = "loess", span = 0.15) +
    scale_colour_viridis_d(option = "turbo") +
    labs(title = sprintf("Rep-order trajectory of infer_dec_q -- %s / seed %s", a, s),
         subtitle = "x-axis is run_id (sequential execution order through the shuffled matrix), not cfg or rep",
         x = "run_id (execution order)", y = "infer_dec_q") +
    theme_light_sf(9)
  save_svg(p3, sprintf("cell_%s_trajectory_light", tag), w = 11, h = 4)

  sec2 <- sprintf("\\subsubsection{Execution-order trajectory}\n\n")
  sec2 <- paste0(sec2, sprintf(
    "The per-config summary above pools all 30 replications of each config regardless of when they ran. Plotting \\texttt{infer\\_dec\\_q} against \\texttt{run\\_id} (the sequential execution order through the shuffled 480-row matrix, not the config or replication index) shows whatever path-dependence exists in the underlying rolling-window/decay-slope estimator directly.\n\n"))
  sec2 <- paste0(sec2, sprintf(
    "\\begin{figure}[h]\n\\centering\n\\includegraphics[width=\\textwidth]{cell_%s_trajectory_light.pdf}\n\\caption{Cell %s/seed %s: infer\\_dec\\_q by execution order, coloured by L8 config, with a loess trend.}\n\\end{figure}\n\n\\clearpage\n",
    tag, texify(a), s))

  fname <- file.path(OUT_SEC, sprintf("cell_%s.tex", tag))
  writeLines(paste0(sec, sec2), fname)
  section_files <- c(section_files, sprintf("cell_%s", tag))
}

# ════════════════════════════════════════════════════════════════════════
# PART B -- PER-ISA DEEP DIVE (3 sections)
# ════════════════════════════════════════════════════════════════════════
cat("\n=== PART B: per-ISA deep dive ===\n")
isa_section_files <- c()
for (a in archs) {
  sub <- filter(all_data, arch == a)
  tag <- a
  cat(sprintf("  ISA %s...\n", a))

  seed_summary <- sub %>% group_by(seed) %>%
    summarise(mean_dec = mean(infer_dec_q), sd_dec = sd(infer_dec_q),
              median_dec = median(infer_dec_q), early_exit_rate = mean(early_exit),
              .groups = "drop")

  kw <- kruskal.test(infer_dec_q ~ seed, data = sub)

  # per-factor means within this ISA
  factor_means <- sub %>%
    summarise(
      ent_hi = mean(infer_dec_q[ent_f=="hi"]), ent_lo = mean(infer_dec_q[ent_f=="lo"]),
      cv_hi  = mean(infer_dec_q[cv_f=="hi"]),  cv_lo  = mean(infer_dec_q[cv_f=="lo"]),
      tmp_hi = mean(infer_dec_q[tmp_f=="hi"]), tmp_lo = mean(infer_dec_q[tmp_f=="lo"]),
      stb_hi = mean(infer_dec_q[stb_f=="hi"]), stb_lo = mean(infer_dec_q[stb_f=="lo"]))

  p <- ggplot(sub, aes(x = seed, y = infer_dec_q, fill = seed)) +
    geom_violin(alpha = 0.8, colour = NA) +
    geom_boxplot(width = 0.15, fill = "white", outlier.size = 0.4) +
    scale_fill_manual(values = c("12345"="#D62728","67890"="#FF7F0E","13579"="#2CA02C"), guide = "none") +
    labs(title = sprintf("%s -- infer_dec_q by seed (violin + box)", a),
         subtitle = sprintf("Kruskal-Wallis H=%.2f, df=%d, p=%.3g", kw$statistic, kw$parameter, kw$p.value),
         x = "seed", y = "infer_dec_q") +
    theme_light_sf(10)
  save_svg(p, sprintf("isa_%s_seed_violin_light", tag), w = 8, h = 5)

  tbl1 <- paste0("\\begin{tabular}{lrrrr}\n\\toprule\nseed & mean & sd & median & early-exit rate \\\\\n\\midrule\n",
    paste(sprintf("%s & %.2f & %.2f & %.1f & %.3f \\\\", seed_summary$seed, seed_summary$mean_dec,
                   seed_summary$sd_dec, seed_summary$median_dec, seed_summary$early_exit_rate), collapse="\n"),
    "\n\\bottomrule\n\\end{tabular}\n")

  tbl2 <- paste0("\\begin{tabular}{lrrrr}\n\\toprule\nfactor & mean (hi) & mean (lo) & difference \\\\\n\\midrule\n",
    sprintf("entropy & %.2f & %.2f & %.2f \\\\\n", factor_means$ent_hi, factor_means$ent_lo, factor_means$ent_hi-factor_means$ent_lo),
    sprintf("coeff.\\ of variation & %.2f & %.2f & %.2f \\\\\n", factor_means$cv_hi, factor_means$cv_lo, factor_means$cv_hi-factor_means$cv_lo),
    sprintf("temporal decay & %.2f & %.2f & %.2f \\\\\n", factor_means$tmp_hi, factor_means$tmp_lo, factor_means$tmp_hi-factor_means$tmp_lo),
    sprintf("stability & %.2f & %.2f & %.2f \\\\\n", factor_means$stb_hi, factor_means$stb_lo, factor_means$stb_hi-factor_means$stb_lo),
    "\\bottomrule\n\\end{tabular}\n")

  sec <- sprintf("\\subsection{Architecture: %s}\n\\label{sec:isa_%s}\n\n", a, tag)
  sec <- paste0(sec, sprintf("All three seeds' cells for %s combined: 1{,}440 rows, zero errors. Seed remains the dominant source of variation within this architecture alone (Kruskal-Wallis $H=%.2f$, $p=%.3g$) -- identical in shape to the other two architectures (see \\S\\ref{sec:isa_amd64}--\\ref{sec:isa_riscv64}).\n\n",
    texify(a), kw$statistic, kw$p.value))
  sec <- paste0(sec, sprintf("\\begin{figure}[h]\n\\centering\n\\includegraphics[width=0.75\\textwidth]{isa_%s_seed_violin_light.pdf}\n\\caption{%s: infer\\_dec\\_q distribution by seed (violin + inner boxplot).}\n\\end{figure}\n\n", tag, texify(a)))
  sec <- paste0(sec, "\\begin{table}[h]\n\\centering\n\\caption{Per-seed summary, ", texify(a), ".}\n", tbl1, "\\end{table}\n\n")
  sec <- paste0(sec, "\\begin{table}[h]\n\\centering\n\\caption{Per-factor main-effect means, ", texify(a), " (all seeds pooled).}\n", tbl2, "\\end{table}\n\n\\clearpage\n")

  fname <- file.path(OUT_SEC, sprintf("isa_%s.tex", tag))
  writeLines(sec, fname)
  isa_section_files <- c(isa_section_files, sprintf("isa_%s", tag))
}

# ════════════════════════════════════════════════════════════════════════
# PART C -- ALL SIX PAIRWISE FACTOR INTERACTIONS, FACETED BY ARCH
# ════════════════════════════════════════════════════════════════════════
cat("\n=== PART C: pairwise factor interactions ===\n")
factor_pairs <- list(
  c("ent_f","cv_f","entropy","coefficient of variation"),
  c("ent_f","tmp_f","entropy","temporal decay"),
  c("ent_f","stb_f","entropy","stability"),
  c("cv_f","tmp_f","coefficient of variation","temporal decay"),
  c("cv_f","stb_f","coefficient of variation","stability"),
  c("tmp_f","stb_f","temporal decay","stability")
)
interaction_section_files <- c()
for (fp in factor_pairs) {
  f1 <- fp[1]; f2 <- fp[2]; f1name <- fp[3]; f2name <- fp[4]
  tag <- paste0(f1, "_", f2)
  cat(sprintf("  Pair %s x %s...\n", f1, f2))

  df <- all_data %>% group_by(arch, .data[[f1]], .data[[f2]]) %>%
    summarise(mean_dec = mean(infer_dec_q), se = sd(infer_dec_q)/sqrt(n()), .groups = "drop")
  names(df)[2:3] <- c("F1","F2")

  # two-way ANOVA within this pair, per arch (does the interaction differ by arch?)
  formula_str <- sprintf("infer_dec_q ~ %s * %s * arch", f1, f2)
  fit <- lm(as.formula(formula_str), data = all_data)
  at <- anova(fit)
  interact_row <- at[grepl(":", rownames(at)) & grepl("arch", rownames(at)) &
                       grepl(f1, rownames(at)) & grepl(f2, rownames(at)), ]

  p <- ggplot(df, aes(x = F2, y = mean_dec, colour = F1, group = F1)) +
    geom_line(linewidth = 1) + geom_point(size = 2.5) +
    geom_errorbar(aes(ymin = mean_dec - se, ymax = mean_dec + se), width = 0.08) +
    facet_wrap(~arch, nrow = 1) +
    scale_colour_manual(values = c(hi = "#D62728", lo = "#4A90D9"), name = f1name) +
    labs(title = sprintf("%s x %s interaction on infer_dec_q", f1name, f2name),
         x = f2name, y = "mean infer_dec_q") +
    theme_light_sf(9)
  save_svg(p, sprintf("interact_%s_light", tag), w = 10, h = 4)

  # Companion: same interaction, early_exit rate as the response
  df_ee <- all_data %>% group_by(arch, .data[[f1]], .data[[f2]]) %>%
    summarise(rate = mean(early_exit), .groups = "drop")
  names(df_ee)[2:3] <- c("F1","F2")
  p_ee <- ggplot(df_ee, aes(x = F2, y = rate, colour = F1, group = F1)) +
    geom_line(linewidth = 1) + geom_point(size = 2.5) +
    facet_wrap(~arch, nrow = 1) +
    scale_colour_manual(values = c(hi = "#D62728", lo = "#4A90D9"), name = f1name) +
    labs(title = sprintf("%s x %s interaction on early_exit rate", f1name, f2name),
         x = f2name, y = "early_exit rate") +
    theme_light_sf(9)
  save_svg(p_ee, sprintf("interact_%s_earlyexit_light", tag), w = 10, h = 4)

  tbl <- paste0("\\begin{tabular}{lllrr}\n\\toprule\narch & ", texify(f1), " & ", texify(f2), " & mean & se \\\\\n\\midrule\n",
    paste(sprintf("%s & %s & %s & %.2f & %.2f \\\\", df$arch, df$F1, df$F2, df$mean_dec, df$se), collapse="\n"),
    "\n\\bottomrule\n\\end{tabular}\n")

  three_way_p <- if(nrow(interact_row) > 0) interact_row[1,"Pr(>F)"] else NA

  sec <- sprintf("\\subsection{%s $\\times$ %s}\n\\label{sec:interact_%s}\n\n", f1name, f2name, tag)
  sec <- paste0(sec, sprintf("Three-way interaction term (%s:%s:arch) in the full model: $p=%.3g$ -- %s.\n\n",
    texify(f1), texify(f2), three_way_p,
    ifelse(is.na(three_way_p) || three_way_p > 0.05, "no evidence the interaction itself depends on architecture", "flagged for follow-up")))
  sec <- paste0(sec, sprintf("\\begin{figure}[h]\n\\centering\n\\includegraphics[width=\\textwidth]{interact_%s_light.pdf}\n\\caption{%s $\\times$ %s interaction on mean infer\\_dec\\_q, faceted by architecture.}\n\\end{figure}\n\n", tag, f1name, f2name))
  sec <- paste0(sec, "\\begin{table}[h]\n\\centering\n\\caption{Cell means, ", f1name, " x ", f2name, " x architecture.}\n", tbl, "\\end{table}\n\n\\clearpage\n")
  sec <- paste0(sec, sprintf("\\begin{figure}[h]\n\\centering\n\\includegraphics[width=\\textwidth]{interact_%s_earlyexit_light.pdf}\n\\caption{%s $\\times$ %s interaction on \\texttt{early\\_exit} rate, faceted by architecture -- same factor pair, binary-outcome response.}\n\\end{figure}\n\n\\clearpage\n", tag, f1name, f2name))

  fname <- file.path(OUT_SEC, sprintf("interact_%s.tex", tag))
  writeLines(sec, fname)
  interaction_section_files <- c(interaction_section_files, sprintf("interact_%s", tag))
}

# ════════════════════════════════════════════════════════════════════════
# PART D -- QUADRATIC / CONTINUOUS-FACTOR RESPONSE (4 factors)
# ════════════════════════════════════════════════════════════════════════
cat("\n=== PART D: quadratic response per factor ===\n")
quad_factors <- list(
  c("ent_in","entropy", "49152"), c("cv_in","coefficient of variation","9830"),
  c("tmp_in","temporal decay","32768"), c("stb_in","stability","32768"))
quad_section_files <- c()
for (qf in quad_factors) {
  fcol <- qf[1]; fname <- qf[2]; fmax <- as.numeric(qf[3])
  tag <- fcol
  cat(sprintf("  Factor %s...\n", fcol))

  all_data$xnorm <- all_data[[fcol]] / fmax
  df <- all_data %>% group_by(arch, xnorm) %>%
    summarise(mean_dec = mean(infer_dec_q), .groups = "drop")

  p <- ggplot(df, aes(x = xnorm, y = mean_dec, colour = arch)) +
    geom_point(size = 3) + geom_line(aes(group = arch), linewidth = 0.8) +
    scale_colour_manual(values = arch_colours, name = "ISA") +
    scale_x_continuous(breaks = c(0,1), labels = c("lo","hi")) +
    labs(title = sprintf("infer_dec_q response to %s (normalized)", fname),
         x = fname, y = "mean infer_dec_q") +
    theme_light_sf(10)
  save_svg(p, sprintf("quad_%s_light", tag), w = 7, h = 5)

  fit <- lm(infer_dec_q ~ xnorm * arch, data = all_data)
  at <- anova(fit)
  lin_p <- at["xnorm", "Pr(>F)"]

  sec <- sprintf("\\subsection{Response to %s}\n\\label{sec:quad_%s}\n\n", fname, tag)
  sec <- paste0(sec, sprintf("Only two levels of %s were sampled (hi/lo) by the L8 factorial design, so a genuine quadratic term is not identifiable from this dataset -- a third, intermediate level would be needed to separate curvature from the linear effect. Reported here as the linear main-effect model instead: linear term $p=%.3g$.\n\n",
    fname, lin_p))
  sec <- paste0(sec, sprintf("\\begin{figure}[h]\n\\centering\n\\includegraphics[width=0.75\\textwidth]{quad_%s_light.pdf}\n\\caption{Mean infer\\_dec\\_q vs.\\ %s level, by architecture.}\n\\end{figure}\n\n\\clearpage\n", tag, fname))

  fname_out <- file.path(OUT_SEC, sprintf("quad_%s.tex", tag))
  writeLines(sec, fname_out)
  quad_section_files <- c(quad_section_files, sprintf("quad_%s", tag))
}

# ════════════════════════════════════════════════════════════════════════
# PART E -- RAW DATA APPENDIX, per cell, full 480-row longtable
# ════════════════════════════════════════════════════════════════════════
cat("\n=== PART E: raw data appendix ===\n")
appendix_section_files <- c()
for (a in archs) for (s in seeds) {
  cell <- filter(all_data, arch == a, seed == s) %>% arrange(run_id)
  tag <- sprintf("%s_seed%s", a, s)
  cat(sprintf("  Appendix cell %s...\n", tag))

  rows_tex <- sprintf("%d & %d & %d & %d & %.2f & %d \\\\",
                       cell$run_id, cell$cfg, cell$rep, cell$ent_in > 0, cell$infer_dec_q, cell$early_exit)
  body <- paste(rows_tex, collapse = "\n")

  sec <- sprintf("\\subsection{Raw data: %s, seed %s}\n\\label{sec:appendix_%s}\n\n", a, s, tag)
  sec <- paste0(sec, sprintf("All 480 rows, execution order (\\texttt{run\\_id}), for cell %s / seed %s.\n\n", texify(a), s))
  sec <- paste0(sec,
    "\\begin{longtable}{rrrrrr}\n",
    "\\toprule\nrun\\_id & cfg & rep & entropy=hi & infer\\_dec\\_q & early\\_exit \\\\\n\\midrule\n",
    "\\endfirsthead\n",
    "\\toprule\nrun\\_id & cfg & rep & entropy=hi & infer\\_dec\\_q & early\\_exit \\\\\n\\midrule\n",
    "\\endhead\n",
    body, "\n",
    "\\bottomrule\n",
    "\\end{longtable}\n\n\\clearpage\n")

  fname <- file.path(OUT_SEC, sprintf("appendix_%s.tex", tag))
  writeLines(sec, fname)
  appendix_section_files <- c(appendix_section_files, sprintf("appendix_%s", tag))
}

# ════════════════════════════════════════════════════════════════════════
# Write the master include list
# ════════════════════════════════════════════════════════════════════════
writeLines(sprintf("\\input{sections/%s}", section_files), file.path(OUT_SEC, "_include_cells.tex"))
writeLines(sprintf("\\input{sections/%s}", isa_section_files), file.path(OUT_SEC, "_include_isa.tex"))
writeLines(sprintf("\\input{sections/%s}", interaction_section_files), file.path(OUT_SEC, "_include_interactions.tex"))
writeLines(sprintf("\\input{sections/%s}", quad_section_files), file.path(OUT_SEC, "_include_quad.tex"))
writeLines(sprintf("\\input{sections/%s}", appendix_section_files), file.path(OUT_SEC, "_include_appendix.tex"))

cat(sprintf("\nGenerated %d cell sections, %d ISA sections, %d interaction sections, %d quadratic sections, %d appendix tables.\n",
            length(section_files), length(isa_section_files), length(interaction_section_files), length(quad_section_files), length(appendix_section_files)))
cat("Done.\n")