#!/usr/bin/env Rscript
# ---------------------------------------------------------------------------
# z93_distribution_charts.R (2026-08-18)
# Publication charts for the /r1a-z93-distribution post, in the site palette.
# Inputs:  runs/z93-ancient-era-clade-2026-08-18.csv (212 ancient males)
#          runs/z93-modern-country-x-stratum-2026-08-18.csv (YFull flags)
# Outputs: SVGs into ~/projects/indoaryan/public/images/blog/z93/
# Style: Pudding/Pew discipline — action titles, direct labels, no legends
# where color is labelled in place, shared scales, ordered panels.
# ---------------------------------------------------------------------------
suppressMessages({ library(ggplot2); library(scales) })

ROOT <- "/home/cdr/projects/sohi-qpadm"
OUT  <- "/home/cdr/projects/indoaryan/public/images/blog/z93"
dir.create(OUT, recursive = TRUE, showWarnings = FALSE)

INK <- "#262019"; MUTED <- "#7a6f5d"; TEAL <- "#356b70"; OCHRE <- "#a76f2b"
BRONZE <- "#6f4a14"; GREY <- "#9a8f7d"; RULE <- "#d9cfba"

theme_paper <- function(base = 12) {
  theme_minimal(base_size = base) +
    theme(
      text = element_text(colour = INK),
      plot.title = element_text(face = "bold", size = base * 1.25, hjust = 0),
      plot.subtitle = element_text(colour = MUTED, size = base * 0.92, hjust = 0),
      plot.caption = element_text(colour = MUTED, size = base * 0.72, hjust = 0),
      axis.text = element_text(colour = MUTED),
      axis.title = element_text(colour = MUTED, size = base * 0.85),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(colour = RULE, linewidth = 0.35),
      strip.text = element_text(colour = INK, face = "bold", size = base * 0.9),
      plot.background = element_rect(fill = "transparent", colour = NA),
      panel.background = element_rect(fill = "transparent", colour = NA),
      legend.position = "none",
      plot.title.position = "plot", plot.caption.position = "plot"
    )
}

anc <- read.csv(file.path(ROOT, "runs/z93-ancient-era-clade-2026-08-18.csv"),
                stringsAsFactors = FALSE)
anc$year <- -as.numeric(anc$year_bce)          # negative = BCE, positive = CE
anc$arm <- with(anc, ifelse(clade == "Z94 Z2124 sibling", "Z2124 arm",
              ifelse(clade == "Z94 Y3 / personal-side", "Y3 arm (toward L657)",
              ifelse(grepl("basal", clade), "basal Z93", "other Z94/Z93 arms"))))
anc$arm <- factor(anc$arm, levels = c("basal Z93", "Z2124 arm",
                                      "Y3 arm (toward L657)", "other Z94/Z93 arms"))
ARM_COL <- c("basal Z93" = GREY, "Z2124 arm" = TEAL,
             "Y3 arm (toward L657)" = OCHRE, "other Z94/Z93 arms" = "#b7a98f")

## 1. HERO — every ancient man as a dot in time -------------------------------
set.seed(7)
anc$row <- ave(anc$year, anc$arm, FUN = function(x) {
  r <- rank(x, ties.method = "first"); ((r - 1) %% 14) + 1 })
band <- c("basal Z93" = 0, "Z2124 arm" = 16, "Y3 arm (toward L657)" = 33,
          "other Z94/Z93 arms" = 39)
anc$ypos <- band[as.character(anc$arm)] + anc$row * 0.9

annot <- data.frame(
  x = c(-2700, -1950, -801, 550),
  y = c(13.5, 31, 36.5, 31.5),
  lab = c("Fatyanovo-region men,\nall basal", "Sintashta era begins:\nZ2124 arm appears",
          "Butkara II, Swat — the lone\npre-CE Y3 man in this dataset", "Hun–Avar era:\n17 men in Hungary"))

p1 <- ggplot(anc, aes(year, ypos, colour = arm)) +
  annotate("rect", xmin = -2100, xmax = -1550, ymin = -1, ymax = 46,
           fill = TEAL, alpha = 0.06) +
  geom_point(size = 1.9, alpha = 0.85) +
  geom_text(data = annot, aes(x, y, label = lab), inherit.aes = FALSE,
            colour = INK, size = 3.0, hjust = 0, lineheight = 0.95) +
  annotate("text", x = -3050, y = band + 5.5, label = names(band),
           colour = ARM_COL[names(band)], fontface = "bold", size = 3.4, hjust = 0) +
  scale_colour_manual(values = ARM_COL) +
  scale_x_continuous(breaks = seq(-3000, 1500, 500),
                     labels = function(v) ifelse(v < 0, paste0(abs(v), " BCE"),
                                          ifelse(v == 0, "0", paste0(v, " CE")))) +
  scale_y_continuous(NULL, breaks = NULL, limits = c(-1, 47)) +
  labs(title = "Every sampled ancient R1a-Z93 man, placed in time",
       subtitle = "212 males from the AADR v66 audit; one dot per man, grouped by tree arm. Shaded band: the Sintashta–Andronovo window.",
       x = NULL,
       caption = "Data: AADR v66 × YFull v14.04 project audit · dot rows within an arm are layout only") +
  theme_paper()
svg(file.path(OUT, "hero-ancient-timeline.svg"), width = 10.4, height = 6.4, bg = "transparent")
print(p1); invisible(dev.off())

## 2. Composition by era ------------------------------------------------------
tab <- as.data.frame(table(anc$era, anc$arm))
names(tab) <- c("era", "arm", "n")
tab$era <- factor(tab$era, levels = sort(unique(as.character(tab$era))))
era_n <- aggregate(n ~ era, tab, sum)

p2 <- ggplot(tab, aes(era, n, fill = arm)) +
  geom_col(position = "fill", width = 0.62) +
  geom_text(data = era_n, aes(era, 1.045, label = paste0("n=", n)),
            inherit.aes = FALSE, colour = MUTED, size = 3) +
  scale_fill_manual(values = ARM_COL) +
  scale_y_continuous(labels = percent_format(), breaks = seq(0, 1, 0.25), limits = c(0, 1.09)) +
  scale_x_discrete(labels = function(v) gsub("^\\d ", "", gsub(" ", "\n", v, fixed = FALSE))) +
  labs(title = "The Z2124 arm takes over the sampled record",
       subtitle = "Share of each tree arm among sampled ancient Z93 males, by era — direction robust, exact rates carry the sampling biases discussed in the text.",
       x = NULL, y = NULL,
       caption = "teal = Z2124 arm · grey = basal Z93 · ochre = Y3 arm · sand = other arms") +
  theme_paper()
svg(file.path(OUT, "era-composition.svg"), width = 10.4, height = 5.4, bg = "transparent")
print(p2); invisible(dev.off())

## 3. Small multiples — the geographic sweep ----------------------------------
top_c <- names(sort(table(anc$entity), decreasing = TRUE))[1:6]
sub <- anc[anc$entity %in% top_c, ]
sub$entity <- factor(sub$entity, levels = top_c)
tab3 <- as.data.frame(table(sub$era, sub$entity))
names(tab3) <- c("era", "entity", "n")
tab3$eranum <- as.integer(substr(as.character(tab3$era), 1, 1))

p3 <- ggplot(tab3, aes(eranum, n)) +
  geom_col(fill = TEAL, width = 0.7) +
  facet_wrap(~entity, nrow = 1) +
  scale_x_continuous(breaks = 1:6, labels = c("pre-\nSint.", "Sint.–\nAndr.", "LBA–\nIron", "Scyth.–\nXiongnu", "Migr.\nera", "Med.+")) +
  labs(title = "Sampling sweeps east, then boomerangs west",
       subtitle = "Ancient Z93 men per era in the six best-sampled countries — Russia leads early, Kazakhstan–China–Mongolia peak in the Scythian era, Hungary spikes with the Hun–Avar return.",
       x = NULL, y = "men",
       caption = "Shared scale across panels; countries ordered by total sample count") +
  theme_paper(11) +
  theme(axis.text.x = element_text(size = 7))
svg(file.path(OUT, "geographic-sweep.svg"), width = 10.4, height = 4.6, bg = "transparent")
print(p3); invisible(dev.off())

## 4. Modern map — FTDNA haplotree census, country x arm -----------------------
mod <- read.csv(file.path(ROOT, "runs/ftdna-modern-country-x-arm-2026-08-18.csv"),
                check.names = FALSE)
mod <- mod[order(-mod$total), ][1:12, ]
arm_cols <- c("Y3 arm (toward L657)", "Z2124 arm", "basal / other Z93-Z94")
long <- data.frame(country = rep(mod$country, length(arm_cols)),
                   arm = rep(arm_cols, each = nrow(mod)),
                   n = c(mod[["Y3 arm (toward L657)"]], mod[["Z2124 arm"]],
                         mod[["basal / other Z93-Z94"]]))
long <- long[long$n > 0, ]
long$arm <- factor(long$arm, levels = arm_cols)
long$country <- factor(long$country, levels = rev(mod$country))
p4 <- ggplot(long, aes(n, country, fill = arm)) +
  geom_col(width = 0.66) +
  scale_fill_manual(values = c("Y3 arm (toward L657)" = OCHRE,
                               "Z2124 arm" = TEAL,
                               "basal / other Z93-Z94" = GREY),
                    name = NULL) +
  labs(title = "The modern map inverts the ancient one",
       subtitle = "FTDNA haplotree census, 1,867 Z93 subclades: kits by self-reported country (top 12; 3,071 country-tagged of 8,417).\nEastern Europe reports zero Y3-arm kits; India runs 135 of 193 the other way.",
       x = "tester kits (FTDNA public haplotree, read 2026-08-18)", y = NULL,
       caption = "Tester geography records who tests, not where clades arose · ~63% of kits report no country and are excluded · YFull cross-check in the text") +
  theme_paper() +
  theme(legend.position = "top", legend.justification = "left")
svg(file.path(OUT, "modern-map.svg"), width = 10.4, height = 5.6, bg = "transparent")
print(p4); invisible(dev.off())

cat("wrote 4 SVGs to", OUT, "\n")

## 5. Map — located ancient Z93 men, colored by arm ---------------------------
suppressMessages(library(maps))
mp <- read.csv(file.path(ROOT, "runs/discover-ftdna-2026-08-18/chart-ancient-map.csv"),
               stringsAsFactors = FALSE)
mp <- mp[mp$year > -4000, ]                       # drop nothing else; all Z93-era
mp$armf <- factor(mp$arm, levels = c("basal/other Z93-Z94", "Z2124 arm", "Y3 arm"))
world <- map_data("world")
set.seed(11)
mp$jlat <- mp$lat + runif(nrow(mp), -0.4, 0.4)    # legibility jitter, disclosed
mp$jlon <- mp$lon + runif(nrow(mp), -0.4, 0.4)
p5 <- ggplot() +
  geom_polygon(data = world, aes(long, lat, group = group),
               fill = "#e9e1d0", colour = "#d9cfba", linewidth = 0.15) +
  geom_point(data = mp, aes(jlon, jlat, colour = armf), size = 1.7, alpha = 0.8) +
  annotate("text", x = 40, y = 59.5, label = "Middle Don &\nVolga forest-steppe",
           colour = INK, size = 2.9, lineheight = 0.9) +
  annotate("text", x = 69, y = 39.5, label = "Swat (Butkara II)", colour = INK, size = 2.9) +
  annotate("text", x = 17, y = 43.5, label = "Carpathian Basin\n(Hun-Avar era)", colour = INK, size = 2.9, lineheight = 0.9) +
  annotate("text", x = 99, y = 52.5, label = "Xiongnu-era\ncemeteries", colour = INK, size = 2.9, lineheight = 0.9) +
  scale_colour_manual(values = c("basal/other Z93-Z94" = GREY, "Z2124 arm" = TEAL,
                                 "Y3 arm" = OCHRE), name = NULL) +
  coord_quickmap(xlim = c(5, 120), ylim = c(20, 62), expand = FALSE) +
  labs(title = "Two hundred located ancient Z93 men",
       subtitle = "Every register sample with coordinates and a date, colored by tree arm. The 46-56°N steppe corridor holds nearly the entire record;\nthe anonymized Akbari-release samples cannot appear here.",
       x = NULL, y = NULL,
       caption = "FTDNA Discover register, n=200 located of 363 · positions jittered ±0.4° for legibility · sites, dates and sources in the published CSV") +
  theme_paper() +
  theme(legend.position = "top", legend.justification = "left",
        axis.text = element_blank(), panel.grid.major = element_blank())
svg(file.path(OUT, "ancient-map.svg"), width = 10.4, height = 6.2, bg = "transparent")
print(p5); invisible(dev.off())

## 6. Branch formation timeline — the radiation burst --------------------------
bf <- read.csv(file.path(ROOT, "runs/discover-ftdna-2026-08-18/chart-branch-formations.csv"),
               stringsAsFactors = FALSE)
p6 <- ggplot(bf, aes(formed)) +
  annotate("rect", xmin = -2500, xmax = -1500, ymin = 0, ymax = Inf, fill = TEAL, alpha = 0.07) +
  geom_histogram(breaks = seq(-3250, 1500, 250), fill = BRONZE, colour = "transparent") +
  annotate("text", x = -1350, y = 21, label = "65 of 106 dated branches\nform in 2500-1500 BCE",
           colour = INK, size = 3.2, lineheight = 0.95, hjust = 0) +
  annotate("text", x = -680, y = 9, label = "4 new surviving branches\nin 1000 BCE - 1 CE",
           colour = MUTED, size = 3.0, lineheight = 0.95, hjust = 0) +
  scale_x_continuous(breaks = seq(-3000, 1500, 500),
                     labels = function(v) ifelse(v < 0, paste0(abs(v), " BCE"),
                                          ifelse(v == 0, "0", paste0(v, " CE")))) +
  labs(title = "Branch formation concentrates in one Bronze Age pulse",
       subtitle = "Formation dates of the 106 dated Z93 subclades in the crawl set (FTDNA sequence-divergence estimates). The trough after 1500 BCE\nsurvives an ascertainment check: median kit counts are identical for pre- and post-trough formations.",
       x = NULL, y = "branches formed per 250 years",
       caption = "Dated clades are the crawl set (larger branches); a formal size-bias null remains open — see the analysis receipt") +
  theme_paper()
svg(file.path(OUT, "branch-formation.svg"), width = 10.4, height = 4.8, bg = "transparent")
print(p6); invisible(dev.off())

## 7. Attestation lag — coalescence vs first physical attestation --------------
lg <- read.csv(file.path(ROOT, "runs/discover-ftdna-2026-08-18/chart-attestation-lag.csv"),
               stringsAsFactors = FALSE)
lg <- lg[order(lg$lag), ]
lg$clade <- factor(lg$clade, levels = lg$clade)
lg$dir <- ifelse(lg$lag > 500, "long gap", "prompt")
p7 <- ggplot(lg) +
  geom_segment(aes(x = tmrca, xend = oldest_sample, y = clade, yend = clade, colour = dir),
               linewidth = 1.1) +
  geom_point(aes(tmrca, clade), colour = INK, size = 2.2) +
  geom_point(aes(oldest_sample, clade), colour = OCHRE, size = 2.2) +
  geom_text(data = lg[lg$lag > 1500, ],
            aes(x = (tmrca + oldest_sample) / 2, y = clade, label = paste0(round(lag), " y")),
            vjust = -0.8, colour = MUTED, size = 2.8) +
  scale_colour_manual(values = c("long gap" = OCHRE, "prompt" = GREY), guide = "none") +
  scale_x_continuous(breaks = seq(-3000, 1500, 500),
                     labels = function(v) ifelse(v < 0, paste0(abs(v), " BCE"),
                                          ifelse(v == 0, "0", paste0(v, " CE")))) +
  labs(title = "Attestation lag: when a lineage coalesced vs. when we first see a body",
       subtitle = "Black dot: clade coalescence (FTDNA TMRCA). Ochre dot: oldest register sample on that clade. The steppe clades are attested\nalmost immediately; L657 waits 1,906 years — the archaeological invisibility of the cremating zone, measured.",
       x = NULL, y = NULL,
       caption = "14 clades: the named backbone plus the five largest lags · negative gaps are stem placements within dating uncertainty") +
  theme_paper()
svg(file.path(OUT, "attestation-lag.svg"), width = 10.4, height = 5.6, bg = "transparent")
print(p7); invisible(dev.off())
cat("charts 5-7 written\n")