---
name: "Historical Semantic Drift Under Corpus-Balanced Resampling"
description: "Use this skill when you need to test whether an apparent time-series trend — specifically, published claims of 20th-century word-meaning drift (gay, awful, terrific, cell, virus, web, mouse, …) — is a genuine semantic signal or a reporting artifact of shifting genre balance in the source corpus. A negative-control / corpus-stress-test design: the same statistical procedure is applied to two Google Books Ngrams corpora that differ only in genre mix (full English vs fiction-only), with a decade-shuffling permutation null, bootstrap 95% CIs, a synthetic stationary negative control, and sensitivity analysis over alternative decade boundaries."
version: "1.0.1"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "historical-semantics", "semantic-drift", "google-ngrams", "permutation-test", "bootstrap", "negative-control", "corpus-robustness"]
python_version: ">=3.8"
dependencies: []
network_access: "Outbound HTTPS to books.google.com (required for first run only; cached responses reused thereafter)"
disk_space: "~2 MB for cache of 160 JSON responses + results.json + report.md"
expected_runtime: "6–12 min first run; <30 s from cache"
env_vars: "none"
data_source: "Google Books Ngrams Viewer JSON API (https://books.google.com/ngrams/json), corpora en-2019 and en-fiction-2019"
data_revision: "Google Books Ngrams 2019 (20200217) release; cached query responses hashed with SHA256 for integrity"
deterministic: "Yes — seeded (seed=42) RNG, sorted dict iteration, response-content hashing"
---

# Historical Semantic Drift Under Corpus-Balanced Resampling

## When to Use This Skill

Use this skill when you need to test whether published claims of 20th-century lexical drift (e.g., "gay shifted from an adjective meaning 'happy' to a noun denoting sexual identity") are robust to the composition of the corpus they were measured on, or whether they are in part an artifact of shifting genre balance in Google Books. It is a **negative-control / corpus-stress-test** design: the same statistical procedure is applied to two corpora that share a data-generating process except for genre balance, and the question is whether the drift estimate survives that restriction.

## Prerequisites

| Requirement | Value |
|-------------|-------|
| Python version | 3.8+ (standard library only) |
| External libraries | None (no pip, no numpy, no scipy, no pandas) |
| Network access | Outbound HTTPS to `books.google.com` on first run only |
| Disk space | ~2 MB for the cache of 160 JSON responses |
| RAM | < 100 MB |
| Runtime (first run) | 6–12 minutes (160 API calls × ~0.35 s delay + retries) |
| Runtime (cached) | < 30 seconds |
| Environment variables | None required |
| Write permissions | Write access to the workspace directory (for `cache/`, `results.json`, `report.md`, `data_manifest.json`) |
| Random seed | Fixed at 42 in the script |

All network calls go through a single entry-point (`fetch_ngram`) with exponential backoff and a 4-attempt retry; if the 4th attempt fails, the script prints a clear error message to `stderr` and exits with status 1. Cached responses are content-addressed and length-checked before reuse.

## Adaptation Guidance

To retarget this analysis to a different lexical set or a different corpus:

- Change `TARGET_WORDS` (in the `DOMAIN CONFIGURATION` block) to the list of words claimed to have drifted. The published cases currently encoded come from Hamilton et al. (2016) and Kulkarni et al. (2015).
- Change `CORPUS_ALL` and `CORPUS_GENRE` to any valid Google Ngrams corpus identifiers (e.g., `en-us-2019`, `en-gb-2019`, `fre-2019`). The statistical method is corpus-agnostic.
- Change `YEAR_START` / `YEAR_END` to analyze a different historical window; change `DECADE_EARLY` / `DECADE_LATE` to compare different end-point decades.
- Change `POS_TAGS` to widen the syntactic grain (add `PRON`, `DET`, `ADP`, `NUM`).
- Change `N_PERMUTATIONS` and `N_BOOTSTRAPS` to control statistical precision (both default to 1000; the skill is tested at 1000).

The download/parse/verify/report pipeline (`load_data()` → `run_analysis()` → `generate_report()`) stays the same. The statistical core — POS-share vectors per decade, cosine drift, decade-shuffling permutation null, bootstrap CI, and cross-corpus agreement — is a general framework for any categorical-share time series.

---

## Overview

Published word-level semantic-drift estimates (Hamilton et al. 2016, Kulkarni et al. 2015) are typically computed on the full Google Books "eng-all" corpus, whose genre composition is known to drift across the 20th century — fiction's share falls, scientific publishing's share rises. If the observed drift of a word's vector is partly driven by that compositional shift, the drift measured on a single-genre subcorpus (e.g., fiction only) should be markedly smaller.

For each target word we query the Google Books Ngrams JSON API for its yearly frequency disaggregated by part-of-speech tag (`WORD_NOUN`, `WORD_ADJ`, `WORD_VERB`, `WORD_ADV`) in two corpora: the full English corpus (`en-2019`) and the fiction-only subcorpus (`en-fiction-2019`). For each (word, corpus) pair we build POS-share vectors averaged over an early decade (1900–1909) and a late decade (1990–1999), and compute drift as cosine distance between them. We then compare:

- **Drift magnitude** in the full corpus vs. the fiction-only subcorpus
- **Decade-shuffling permutation null** — shuffle the 100 years within each (word, corpus) trajectory and recompute drift 1000 times
- **Bootstrap 95% CI** on drift — resample years within each decade bucket with replacement, 1000 times
- **Cross-corpus agreement** — Spearman rank correlation of per-word drift between the two corpora

A large corpus-dependent gap, combined with weak cross-corpus rank agreement, would suggest that published single-corpus drift estimates are confounded by genre mix.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_historical-embedding-semantic-drift-under-corpus-balanced-re/cache
```

**Expected output:** No output (directory created silently).

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_historical-embedding-semantic-drift-under-corpus-balanced-re/script.py
#!/usr/bin/env python3
"""
Historical Semantic Drift Under Corpus-Balanced Resampling.

For each target word claimed in the published literature to exhibit 20th-century
semantic drift, query the Google Books Ngrams JSON API for POS-tagged annual
frequencies in the full English corpus (en-2019) and the fiction-only subcorpus
(en-fiction-2019). Compute POS-share drift between early and late decade,
permutation null by decade shuffling, bootstrap 95% CI, and cross-corpus
rank agreement.

Python 3.8+ standard library only.
"""

import hashlib
import json
import math
import os
import random
import sys
import time
import urllib.error
import urllib.parse
import urllib.request
from collections import defaultdict
from datetime import datetime

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

# Target words: 20 words that Hamilton et al. (2016, ACL) and Kulkarni et al.
# (2015, WWW) identify as prominent 20th-century semantic drifters.
TARGET_WORDS = [
    "gay", "awful", "terrific", "nice", "cell",
    "virus", "broadcast", "mouse", "browser", "web",
    "record", "tape", "radio", "program", "hack",
    "chip", "server", "net", "message", "wicked",
]

# Google Ngrams corpus identifiers. The 2019 release was published 2020-02-17.
CORPUS_ALL   = "en-2019"
CORPUS_GENRE = "en-fiction-2019"

# POS tags supported by the Ngrams Viewer.
POS_TAGS = ["NOUN", "ADJ", "VERB", "ADV"]

# Historical window.
YEAR_START = 1900
YEAR_END   = 1999  # inclusive; 100 years total

# Decade end-points to compare.
DECADE_EARLY = (1900, 1909)
DECADE_LATE  = (1990, 1999)

# Sensitivity-analysis alternative decade pairs.
SENSITIVITY_DECADE_PAIRS = [
    ((1910, 1919), (1980, 1989)),
    ((1920, 1929), (1970, 1979)),
]

# Statistical parameters.
N_PERMUTATIONS = 1000
N_BOOTSTRAPS   = 1000
RANDOM_SEED    = 42

# Network behavior.
API_BASE       = "https://books.google.com/ngrams/json"
USER_AGENT     = "Mozilla/5.0 (compatible; Claw4S-research-bot/1.0)"
REQUEST_SLEEP  = 0.35      # polite inter-request delay
MAX_RETRIES    = 4
TIMEOUT        = 30

# Output filenames.
RESULTS_FILE = "results.json"
REPORT_FILE  = "report.md"

# ═══════════════════════════════════════════════════════════════
# End of DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════


SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(SCRIPT_DIR, "cache")


# ---------- Statistical helpers ----------

def cosine_similarity(u, v):
    dot = sum(a * b for a, b in zip(u, v))
    nu = math.sqrt(sum(a * a for a in u))
    nv = math.sqrt(sum(b * b for b in v))
    if nu == 0 or nv == 0:
        return 0.0
    return dot / (nu * nv)


def cosine_distance(u, v):
    return 1.0 - cosine_similarity(u, v)


def normalize_share(counts):
    """Convert a list of non-negative counts into a probability vector.
    Returns a uniform vector if total is 0 (so cosine is defined)."""
    total = sum(counts)
    if total <= 0:
        return [1.0 / len(counts)] * len(counts)
    return [c / total for c in counts]


def mean(xs):
    return sum(xs) / len(xs) if xs else 0.0


def percentile_sorted(xs_sorted, p):
    if not xs_sorted:
        return 0.0
    k = (len(xs_sorted) - 1) * p / 100.0
    f = int(k)
    c = min(f + 1, len(xs_sorted) - 1)
    return xs_sorted[f] + (k - f) * (xs_sorted[c] - xs_sorted[f])


def spearman_rho(x, y):
    """Spearman rank correlation. Ties broken by average rank."""
    n = len(x)
    if n < 2:
        return 0.0

    def ranks(v):
        indexed = sorted(range(n), key=lambda i: v[i])
        r = [0.0] * n
        i = 0
        while i < n:
            j = i
            while j + 1 < n and v[indexed[j + 1]] == v[indexed[i]]:
                j += 1
            avg_rank = (i + j) / 2.0 + 1.0
            for k in range(i, j + 1):
                r[indexed[k]] = avg_rank
            i = j + 1
        return r

    rx = ranks(x)
    ry = ranks(y)
    mx = sum(rx) / n
    my = sum(ry) / n
    num = sum((rx[i] - mx) * (ry[i] - my) for i in range(n))
    dx = math.sqrt(sum((rx[i] - mx) ** 2 for i in range(n)))
    dy = math.sqrt(sum((ry[i] - my) ** 2 for i in range(n)))
    if dx == 0 or dy == 0:
        return 0.0
    return num / (dx * dy)


def spearman_permutation_p(x, y, n_permutations, rng):
    """Two-sided permutation p-value on Spearman rho by shuffling y."""
    observed = spearman_rho(x, y)
    n = len(y)
    count = 0
    y_shuffled = list(y)
    for _ in range(n_permutations):
        rng.shuffle(y_shuffled)
        rho = spearman_rho(x, y_shuffled)
        if abs(rho) >= abs(observed) - 1e-12:
            count += 1
    p = (count + 1) / (n_permutations + 1)
    return observed, p


def sign_flip_permutation_of_mean(diffs, n_permutations, rng):
    """Sign-flip permutation test for paired differences. Statistic: mean(diff).
    Under the null of symmetric zero-centered diffs, flipping each sign with
    prob 0.5 yields an exact distribution of mean(diff). Returns
    (observed_mean, two_sided_p_value)."""
    observed = mean(diffs)
    n = len(diffs)
    abs_vals = [abs(d) for d in diffs]
    count = 0
    for _ in range(n_permutations):
        s = sum((1 if rng.random() < 0.5 else -1) * a for a in abs_vals) / n
        if abs(s) >= abs(observed) - 1e-12:
            count += 1
    p = (count + 1) / (n_permutations + 1)
    return observed, p


def pearson_r(x, y):
    """Pearson correlation coefficient, stdlib."""
    n = len(x)
    if n < 2:
        return 0.0
    mx = sum(x) / n
    my = sum(y) / n
    num = sum((x[i] - mx) * (y[i] - my) for i in range(n))
    dx = math.sqrt(sum((x[i] - mx) ** 2 for i in range(n)))
    dy = math.sqrt(sum((y[i] - my) ** 2 for i in range(n)))
    if dx == 0 or dy == 0:
        return 0.0
    return num / (dx * dy)


# ---------- SHA256 caching ----------

def sha256_bytes(b):
    h = hashlib.sha256()
    h.update(b)
    return h.hexdigest()


def sha256_file(path):
    h = hashlib.sha256()
    with open(path, "rb") as f:
        while True:
            chunk = f.read(65536)
            if not chunk:
                break
            h.update(chunk)
    return h.hexdigest()


# ---------- Google Ngrams API ----------

def _cache_path(content, corpus):
    key = f"{content}__{corpus}__{YEAR_START}_{YEAR_END}"
    safe = hashlib.sha1(key.encode("utf-8")).hexdigest()
    return os.path.join(CACHE_DIR, f"{safe}.json")


def fetch_ngram(content, corpus):
    """Fetch a single POS-tagged ngram trajectory. Returns list of floats of length
    YEAR_END - YEAR_START + 1. Uses on-disk cache."""
    path = _cache_path(content, corpus)
    if os.path.exists(path):
        with open(path, "r") as f:
            cached = json.load(f)
        ts = cached.get("timeseries", [])
        expected_len = YEAR_END - YEAR_START + 1
        if len(ts) == expected_len:
            return ts

    params = {
        "content":    content,
        "year_start": str(YEAR_START),
        "year_end":   str(YEAR_END),
        "corpus":     corpus,
        "smoothing":  "0",
    }
    url = API_BASE + "?" + urllib.parse.urlencode(params)

    last_err = None
    body = None
    payload = None
    for attempt in range(1, MAX_RETRIES + 1):
        try:
            time.sleep(REQUEST_SLEEP)
            req = urllib.request.Request(url, headers={
                "User-Agent": USER_AGENT,
                "Accept":     "application/json",
            })
            with urllib.request.urlopen(req, timeout=TIMEOUT) as resp:
                body = resp.read()
            payload = json.loads(body.decode("utf-8"))
            break
        except (urllib.error.URLError, urllib.error.HTTPError, json.JSONDecodeError, OSError) as e:
            last_err = e
            back = min(2.0 ** attempt, 16.0)
            print(f"    WARN: fetch {content} from {corpus} attempt {attempt}/{MAX_RETRIES} failed ({e}); "
                  f"retrying in {back:.1f}s", file=sys.stderr)
            time.sleep(back)
    if payload is None:
        raise RuntimeError(
            f"Failed to fetch `{content}` from corpus `{corpus}` after {MAX_RETRIES} retries. "
            f"Last error: {last_err}. Check network access to {API_BASE} or rerun later; "
            f"previously-cached queries persist on disk."
        )

    expected_len = YEAR_END - YEAR_START + 1
    ts = []
    if isinstance(payload, list) and payload:
        entry = payload[0]
        if isinstance(entry, dict) and entry.get("type") == "NGRAM":
            ts = list(entry.get("timeseries", []))
    if len(ts) != expected_len:
        ts = [0.0] * expected_len

    os.makedirs(CACHE_DIR, exist_ok=True)
    with open(path, "w") as f:
        json.dump({
            "content": content, "corpus": corpus,
            "year_start": YEAR_START, "year_end": YEAR_END,
            "timeseries": ts,
            "sha256_response": sha256_bytes(body),
            "fetched_at": datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ"),
        }, f, indent=2)
    return ts


# ---------- Load ----------

def load_data():
    """Populate a nested dict: data[word][corpus][pos] = list of yearly frequencies.

    Uses the Google Books Ngrams Viewer JSON API. All responses are cached on
    disk. Returns the nested dict plus a manifest (cache hashes, fetch dates).
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    data = {}
    manifest = {"fetches": [], "corpora": [CORPUS_ALL, CORPUS_GENRE], "pos_tags": POS_TAGS}

    total_queries = len(TARGET_WORDS) * 2 * len(POS_TAGS)
    done = 0
    t0 = time.time()
    for word in TARGET_WORDS:
        data[word] = {}
        for corpus in [CORPUS_ALL, CORPUS_GENRE]:
            data[word][corpus] = {}
            for pos in POS_TAGS:
                content = f"{word}_{pos}"
                ts = fetch_ngram(content, corpus)
                data[word][corpus][pos] = ts
                done += 1
                # compact progress
                if done % 20 == 0 or done == total_queries:
                    elapsed = time.time() - t0
                    print(f"    fetched {done}/{total_queries} queries "
                          f"(elapsed {elapsed:.1f}s)")
                # Hash the actual time-series data (not the cache file, which
                # contains a non-deterministic `fetched_at` timestamp). This
                # guarantees corpus_sha256 is identical across reruns whenever
                # the underlying API payload is identical.
                ts_bytes = json.dumps(ts, sort_keys=True, separators=(",", ":")).encode("utf-8")
                manifest["fetches"].append({
                    "content": content,
                    "corpus": corpus,
                    "content_sha256": sha256_bytes(ts_bytes),
                })

    # --- Corpus-liveness sanity check ---
    # If more than 50% of fetched trajectories are all-zero, the corpus
    # identifier has almost certainly been renamed or the API changed.
    n_zero = 0
    n_total = 0
    for word in TARGET_WORDS:
        for corpus in [CORPUS_ALL, CORPUS_GENRE]:
            for pos in POS_TAGS:
                n_total += 1
                if all(v == 0.0 for v in data[word][corpus][pos]):
                    n_zero += 1
    if n_zero > n_total * 0.5:
        raise RuntimeError(
            f"Corpus-liveness check failed: {n_zero}/{n_total} trajectories "
            f"are all-zero. The Ngrams corpora `{CORPUS_ALL}` or "
            f"`{CORPUS_GENRE}` may have been renamed or deprecated."
        )
    manifest["n_all_zero_trajectories"] = n_zero
    manifest["n_total_trajectories"] = n_total

    # --- Global canonical corpus SHA256 across all response payloads ---
    #     Uses content-addressed hashes (not file-mtime hashes) so corpus_sha256
    #     is stable across reruns whenever the API returns identical payloads.
    h = hashlib.sha256()
    for entry in sorted(manifest["fetches"], key=lambda d: (d["content"], d["corpus"])):
        h.update(entry["content"].encode("utf-8"))
        h.update(b"|")
        h.update(entry["corpus"].encode("utf-8"))
        h.update(b"|")
        h.update(entry["content_sha256"].encode("utf-8"))
        h.update(b"\n")
    manifest["corpus_sha256"] = h.hexdigest()

    return data, manifest


# ---------- Core analysis ----------

def decade_share_vector(data_word_corpus, decade):
    """Average POS-share vector over years in [decade[0], decade[1]]."""
    a, b = decade
    idx_a = a - YEAR_START
    idx_b = b - YEAR_START
    yearly_totals = []
    for y_idx in range(idx_a, idx_b + 1):
        counts = [data_word_corpus[pos][y_idx] for pos in POS_TAGS]
        yearly_totals.append(counts)
    # Mean over years, then normalize.
    means = [mean([y[i] for y in yearly_totals]) for i in range(len(POS_TAGS))]
    return normalize_share(means)


def drift_for_word_corpus(data_word_corpus, d_early, d_late):
    ve = decade_share_vector(data_word_corpus, d_early)
    vl = decade_share_vector(data_word_corpus, d_late)
    return cosine_distance(ve, vl), ve, vl


def permutation_null_drift(data_word_corpus, d_early, d_late, n_permutations, rng):
    """Shuffle the 100 years within each POS trajectory (same permutation
    applied across POS to preserve per-year composition), recompute drift.
    Returns observed, p (observed or larger), and null distribution quantiles.
    """
    observed, _, _ = drift_for_word_corpus(data_word_corpus, d_early, d_late)

    year_indices = list(range(YEAR_END - YEAR_START + 1))
    permuted_drifts = []

    for _ in range(n_permutations):
        rng.shuffle(year_indices)
        reindexed = {pos: [data_word_corpus[pos][i] for i in year_indices] for pos in POS_TAGS}
        d, _, _ = drift_for_word_corpus(reindexed, d_early, d_late)
        permuted_drifts.append(d)

    count = sum(1 for d in permuted_drifts if d >= observed - 1e-12)
    p = (count + 1) / (n_permutations + 1)
    permuted_sorted = sorted(permuted_drifts)
    return {
        "observed": observed,
        "p_value": p,
        "null_median": percentile_sorted(permuted_sorted, 50),
        "null_p95":    percentile_sorted(permuted_sorted, 95),
        "null_mean":   mean(permuted_drifts),
    }


def block_permutation_null_drift(data_word_corpus, d_early, d_late, n_permutations, rng, block_size=10):
    """Block-permutation null: shuffle 10-year blocks rather than individual
    years, preserving short-range autocorrelation. Returns p-value."""
    observed, _, _ = drift_for_word_corpus(data_word_corpus, d_early, d_late)
    n_years = YEAR_END - YEAR_START + 1
    n_blocks = n_years // block_size
    # Drop any trailing year(s) that don't fit evenly into blocks.
    trimmed_len = n_blocks * block_size
    count = 0
    for _ in range(n_permutations):
        block_order = list(range(n_blocks))
        rng.shuffle(block_order)
        new_indices = []
        for b in block_order:
            new_indices.extend(range(b * block_size, (b + 1) * block_size))
        # Pad trailing years (unchanged position) back in to keep length.
        new_indices.extend(range(trimmed_len, n_years))
        reindexed = {pos: [data_word_corpus[pos][i] for i in new_indices] for pos in POS_TAGS}
        d, _, _ = drift_for_word_corpus(reindexed, d_early, d_late)
        if d >= observed - 1e-12:
            count += 1
    return (count + 1) / (n_permutations + 1)


def benjamini_hochberg(pvals):
    """Benjamini-Hochberg FDR-adjusted p-values, stdlib. Returns list of q-values
    in the original order of pvals. Equivalent to R's p.adjust(..., 'BH')."""
    n = len(pvals)
    if n == 0:
        return []
    ordered = sorted(range(n), key=lambda i: pvals[i])
    q_ordered = [0.0] * n
    min_q = 1.0
    for rank_from_bottom, idx in enumerate(reversed(ordered)):
        rank = n - rank_from_bottom  # rank from largest down
        q = pvals[idx] * n / rank
        if q < min_q:
            min_q = q
        q_ordered[idx] = min(min_q, 1.0)
    return q_ordered


def bootstrap_ci_drift(data_word_corpus, d_early, d_late, n_bootstraps, rng):
    """Bootstrap years within each decade bucket, 95% CI on cosine drift."""
    def sample_decade(decade):
        a, b = decade
        years = list(range(a - YEAR_START, b - YEAR_START + 1))
        return [rng.choice(years) for _ in range(len(years))]

    samples = []
    for _ in range(n_bootstraps):
        idx_early = sample_decade(d_early)
        idx_late  = sample_decade(d_late)
        means_e = [mean([data_word_corpus[pos][i] for i in idx_early]) for pos in POS_TAGS]
        means_l = [mean([data_word_corpus[pos][i] for i in idx_late])  for pos in POS_TAGS]
        ve = normalize_share(means_e)
        vl = normalize_share(means_l)
        samples.append(cosine_distance(ve, vl))

    samples.sort()
    return {
        "mean":    mean(samples),
        "ci_lo":   percentile_sorted(samples, 2.5),
        "ci_hi":   percentile_sorted(samples, 97.5),
        "ci_10":   percentile_sorted(samples, 10),
        "ci_90":   percentile_sorted(samples, 90),
    }


def run_analysis(data, manifest=None):
    rng = random.Random(RANDOM_SEED)

    # 1. Observed drift + permutation p + bootstrap CI for every (word, corpus).
    per_word = {}
    for word in TARGET_WORDS:
        per_word[word] = {}
        for corpus in [CORPUS_ALL, CORPUS_GENRE]:
            d_obs, ve, vl = drift_for_word_corpus(data[word][corpus], DECADE_EARLY, DECADE_LATE)
            perm = permutation_null_drift(data[word][corpus], DECADE_EARLY, DECADE_LATE,
                                          N_PERMUTATIONS, rng)
            block_p = block_permutation_null_drift(data[word][corpus], DECADE_EARLY, DECADE_LATE,
                                                    N_PERMUTATIONS, rng, block_size=10)
            boot = bootstrap_ci_drift(data[word][corpus], DECADE_EARLY, DECADE_LATE,
                                      N_BOOTSTRAPS, rng)
            per_word[word][corpus] = {
                "drift":            round(d_obs, 6),
                "share_early":      [round(x, 6) for x in ve],
                "share_late":       [round(x, 6) for x in vl],
                "permutation_p":    round(perm["p_value"], 6),
                "block_permutation_p": round(block_p, 6),
                "null_median":      round(perm["null_median"], 6),
                "null_p95":         round(perm["null_p95"], 6),
                "null_mean":        round(perm["null_mean"], 6),
                "bootstrap_mean":   round(boot["mean"], 6),
                "bootstrap_ci_lo":  round(boot["ci_lo"], 6),
                "bootstrap_ci_hi":  round(boot["ci_hi"], 6),
            }

    # Benjamini-Hochberg FDR adjustment across the 40 (word, corpus) permutation p-values.
    p_list = []
    cells = []
    for word in TARGET_WORDS:
        for corpus in [CORPUS_ALL, CORPUS_GENRE]:
            p_list.append(per_word[word][corpus]["permutation_p"])
            cells.append((word, corpus))
    q_list = benjamini_hochberg(p_list)
    for (word, corpus), q in zip(cells, q_list):
        per_word[word][corpus]["permutation_q_bh"] = round(q, 6)

    # 2. Cross-corpus rank agreement of drift scores.
    drift_all      = [per_word[w][CORPUS_ALL]["drift"]   for w in TARGET_WORDS]
    drift_fiction  = [per_word[w][CORPUS_GENRE]["drift"] for w in TARGET_WORDS]
    rho_all_vs_fic, p_rho = spearman_permutation_p(drift_all, drift_fiction,
                                                    N_PERMUTATIONS, rng)

    # 3. Paired difference: drift_all(w) - drift_fiction(w).
    diffs = [drift_all[i] - drift_fiction[i] for i in range(len(TARGET_WORDS))]
    diff_mean, diff_p = sign_flip_permutation_of_mean(diffs, N_PERMUTATIONS, rng)
    # Bootstrap CI on the mean of diffs.
    diff_boots = []
    for _ in range(N_BOOTSTRAPS):
        resample = [diffs[rng.randrange(len(diffs))] for _ in range(len(diffs))]
        diff_boots.append(mean(resample))
    diff_boots.sort()
    diff_ci_lo = percentile_sorted(diff_boots, 2.5)
    diff_ci_hi = percentile_sorted(diff_boots, 97.5)

    # 4. Count of words where drift_all "survives" corpus restriction.
    #    Survival criterion 1: permutation-significant in BOTH corpora (p<=0.05)
    #    Survival criterion 2: fiction drift within bootstrap CI of all-genre drift.
    n_sig_all     = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_ALL]["permutation_p"]   <= 0.05)
    n_sig_fiction = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_GENRE]["permutation_p"] <= 0.05)
    n_sig_both    = sum(1 for w in TARGET_WORDS
                        if per_word[w][CORPUS_ALL]["permutation_p"]   <= 0.05
                        and per_word[w][CORPUS_GENRE]["permutation_p"] <= 0.05)
    n_fdr_sig_all     = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_ALL]["permutation_q_bh"]   <= 0.05)
    n_fdr_sig_fiction = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_GENRE]["permutation_q_bh"] <= 0.05)
    n_block_sig_all     = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_ALL]["block_permutation_p"]   <= 0.05)
    n_block_sig_fiction = sum(1 for w in TARGET_WORDS if per_word[w][CORPUS_GENRE]["block_permutation_p"] <= 0.05)

    # Survival criterion 2: does fiction-corpus drift fall inside all-corpus bootstrap CI?
    n_drift_agreeing = 0
    for w in TARGET_WORDS:
        lo = per_word[w][CORPUS_ALL]["bootstrap_ci_lo"]
        hi = per_word[w][CORPUS_ALL]["bootstrap_ci_hi"]
        f  = per_word[w][CORPUS_GENRE]["drift"]
        if lo <= f <= hi:
            n_drift_agreeing += 1

    # 4b. Trimmed-mean paired diff (drop the 1 most positive and 1 most negative).
    diffs_sorted_indices = sorted(range(len(diffs)), key=lambda i: diffs[i])
    trim_indices = set(diffs_sorted_indices[:1] + diffs_sorted_indices[-1:])
    trimmed_diffs = [diffs[i] for i in range(len(diffs)) if i not in trim_indices]
    trimmed_mean_diff = mean(trimmed_diffs)
    trimmed_drift_all     = [drift_all[i]     for i in range(len(drift_all))     if i not in trim_indices]
    trimmed_drift_fiction = [drift_fiction[i] for i in range(len(drift_fiction)) if i not in trim_indices]

    # 4c. Per-word frequency-trajectory Pearson correlation between corpora.
    #     For each word, total freq per year = sum across POS tags. If fiction and
    #     all-English trajectories correlate perfectly, the word is used identically
    #     in both. A low correlation is a direct fingerprint of genre-driven
    #     usage difference and a red flag for corpus-dependent drift estimates.
    trajectory_correlations = {}
    for w in TARGET_WORDS:
        traj_all = [sum(data[w][CORPUS_ALL][pos][y]   for pos in POS_TAGS) for y in range(YEAR_END - YEAR_START + 1)]
        traj_fic = [sum(data[w][CORPUS_GENRE][pos][y] for pos in POS_TAGS) for y in range(YEAR_END - YEAR_START + 1)]
        trajectory_correlations[w] = round(pearson_r(traj_all, traj_fic), 4)
    traj_r_values = list(trajectory_correlations.values())
    n_high_traj_r = sum(1 for r in traj_r_values if r >= 0.7)
    n_low_traj_r  = sum(1 for r in traj_r_values if r <  0.5)

    # 5. Sensitivity analysis: alternative decade boundaries.
    sensitivity = []
    for d_e, d_l in SENSITIVITY_DECADE_PAIRS:
        drift_all_s = [drift_for_word_corpus(data[w][CORPUS_ALL],   d_e, d_l)[0] for w in TARGET_WORDS]
        drift_fic_s = [drift_for_word_corpus(data[w][CORPUS_GENRE], d_e, d_l)[0] for w in TARGET_WORDS]
        rho_s, p_s = spearman_permutation_p(drift_all_s, drift_fic_s, 500, rng)
        diffs_s = [drift_all_s[i] - drift_fic_s[i] for i in range(len(TARGET_WORDS))]
        sensitivity.append({
            "early_decade":    list(d_e),
            "late_decade":     list(d_l),
            "mean_drift_all":     round(mean(drift_all_s), 6),
            "mean_drift_fiction": round(mean(drift_fic_s), 6),
            "mean_diff":          round(mean(diffs_s), 6),
            "spearman_rho":       round(rho_s, 4),
            "spearman_p":         round(p_s, 4),
        })

    # 5b. Negative / falsification control.
    #     Construct a synthetic "stationary" word whose annual POS-share trajectory
    #     is identically distributed in every year (independently resampled from
    #     a fixed Dirichlet-like mean). Under the null of no temporal drift the
    #     observed drift should be small and the permutation p-value should be
    #     uniformly distributed on (0,1). We expect the null-control drift to be
    #     much smaller than the observed mean real-word drift and the permutation
    #     p-value to be above 0.05 on this control in at least one of the two
    #     corpora — if either of those is false, the test machinery itself is
    #     biased toward inventing drift.
    negative_control_rng = random.Random(RANDOM_SEED + 1)
    fixed_mean_share = [0.50, 0.20, 0.20, 0.10]  # a plausible POS mix; not varied
    control_annual_counts = {pos: [] for pos in POS_TAGS}
    n_years = YEAR_END - YEAR_START + 1
    for _ in range(n_years):
        total = 1000
        per_pos = []
        remainder = total
        for p_idx, share in enumerate(fixed_mean_share):
            if p_idx == len(fixed_mean_share) - 1:
                jitter = remainder
            else:
                mu = share * total
                jitter = int(round(negative_control_rng.gauss(mu, max(1.0, mu * 0.08))))
                jitter = max(0, min(jitter, remainder))
                remainder -= jitter
            per_pos.append(jitter)
        for p_idx, pos in enumerate(POS_TAGS):
            control_annual_counts[pos].append(float(per_pos[p_idx]))
    control_drift, _, _ = drift_for_word_corpus(control_annual_counts, DECADE_EARLY, DECADE_LATE)
    control_perm = permutation_null_drift(control_annual_counts, DECADE_EARLY, DECADE_LATE,
                                          N_PERMUTATIONS, rng)
    negative_control = {
        "description": "Synthetic stationary word: POS-shares drawn from a fixed distribution every year; should show near-zero drift under a valid permutation null.",
        "drift":              round(control_drift, 6),
        "permutation_p":      round(control_perm["p_value"], 6),
        "null_mean":          round(control_perm["null_mean"], 6),
        "null_p95":           round(control_perm["null_p95"], 6),
        "ratio_observed_to_null_p95": round(control_drift / control_perm["null_p95"], 6)
                               if control_perm["null_p95"] > 0 else None,
    }

    # 6. Aggregate drift statistics.
    results = {
        "target_words":           list(TARGET_WORDS),
        "corpora":                [CORPUS_ALL, CORPUS_GENRE],
        "pos_tags":               list(POS_TAGS),
        "decade_early":           list(DECADE_EARLY),
        "decade_late":            list(DECADE_LATE),
        "n_permutations":         N_PERMUTATIONS,
        "n_bootstraps":           N_BOOTSTRAPS,
        "random_seed":            RANDOM_SEED,

        "per_word":               per_word,

        "aggregate": {
            "mean_drift_all":        round(mean(drift_all), 6),
            "mean_drift_fiction":    round(mean(drift_fiction), 6),
            "median_drift_all":      round(percentile_sorted(sorted(drift_all), 50), 6),
            "median_drift_fiction":  round(percentile_sorted(sorted(drift_fiction), 50), 6),
            "mean_paired_diff":      round(diff_mean, 6),
            "paired_diff_p":         round(diff_p, 6),
            "paired_diff_ci_lo":     round(diff_ci_lo, 6),
            "paired_diff_ci_hi":     round(diff_ci_hi, 6),
            "spearman_rho_all_vs_fiction":    round(rho_all_vs_fic, 4),
            "spearman_p_all_vs_fiction":      round(p_rho, 6),
            "n_words":                        len(TARGET_WORDS),
            "n_sig_all":                      n_sig_all,
            "n_sig_fiction":                  n_sig_fiction,
            "n_sig_both":                     n_sig_both,
            "n_fiction_in_all_bootstrap_ci":  n_drift_agreeing,
            "n_fdr_sig_all":                  n_fdr_sig_all,
            "n_fdr_sig_fiction":              n_fdr_sig_fiction,
            "n_block_sig_all":                n_block_sig_all,
            "n_block_sig_fiction":            n_block_sig_fiction,
            "trimmed_mean_paired_diff":       round(trimmed_mean_diff, 6),
            "trimmed_mean_drift_all":         round(mean(trimmed_drift_all), 6),
            "trimmed_mean_drift_fiction":     round(mean(trimmed_drift_fiction), 6),
            "n_trajectory_r_ge_0.7":          n_high_traj_r,
            "n_trajectory_r_lt_0.5":          n_low_traj_r,
            "mean_trajectory_r":              round(mean(traj_r_values), 4),
            "median_trajectory_r":            round(percentile_sorted(sorted(traj_r_values), 50), 4),
        },

        "trajectory_correlations":  trajectory_correlations,

        "negative_control":       negative_control,

        "sensitivity":            sensitivity,
        "corpus_sha256":          manifest.get("corpus_sha256") if manifest else None,
        "n_all_zero_trajectories": manifest.get("n_all_zero_trajectories") if manifest else None,
        "n_total_trajectories":    manifest.get("n_total_trajectories") if manifest else None,
        "limitations": [
            "Drift is operationalized as change in POS-share vector, a coarse-grained proxy for semantic drift; it does not detect meaning shifts that leave the POS distribution invariant (e.g., a noun whose referent changes).",
            "The Google Books corpus has known biases: OCR errors, scanning-coverage changes over time, overrepresentation of academic material after 1950, and genre labels that are not vetted by linguists.",
            "The en-fiction-2019 subcorpus is a proxy for 'a single genre'; fiction itself is heterogeneous and its internal composition also shifts over the 20th century.",
            "Permutation by decade-shuffling tests a null of no temporal ordering; it does not control for year-to-year autocorrelation driven by language-independent publishing trends.",
            "POS tagging in the Google Ngrams corpus is produced by an automatic tagger whose error rate is non-zero; rare POS assignments for a word may reflect tagger noise.",
            "Target words are drawn from a specific pair of papers (Hamilton 2016, Kulkarni 2015); the conclusion may not generalize to unannotated words with subtler drift.",
            "The analysis does NOT replicate the embedding-based drift scores in the source papers; cosine distance over POS shares is a cheaper, coarser instrument and a negative result here does not falsify the original embedding findings.",
            "Words with very low overall frequency (total < 100 occurrences per decade) will have unstable POS-share estimates; we do not filter or exclude such words, and in rare cases their bootstrap CIs may therefore be very wide.",
            "The corpus is restricted to English; similar drift questions in non-English corpora (e.g., `fre-2019`, `ger-2019`) would require re-running with changed corpus identifiers and may surface different artifacts.",
            "Two of the target words (`browser`, `server`, `net`) have essentially no attested usage before 1960, so their early-decade POS-share vectors are estimated from near-zero counts and their drift scores may be dominated by the Laplace-like uniform fallback in `normalize_share`.",
            "The negative-control drift is non-zero because the synthetic word has small Gaussian jitter on POS shares; we require only that it be an order of magnitude smaller than real-word drift, not that it be zero — a stricter control would null-pass deterministically.",
        ],
    }
    return results


# ---------- Report ----------

def generate_report(results, workspace):
    # results.json
    results_path = os.path.join(workspace, RESULTS_FILE)
    with open(results_path, "w") as f:
        json.dump(results, f, indent=2, default=str)

    # report.md
    lines = []
    agg = results["aggregate"]
    lines.append("# Historical Semantic Drift Under Corpus-Balanced Resampling — Report\n")
    lines.append(f"**Target words:** {len(results['target_words'])} "
                 f"({', '.join(results['target_words'])})")
    lines.append(f"**Corpora:** `{results['corpora'][0]}` (full English) vs "
                 f"`{results['corpora'][1]}` (fiction only)")
    lines.append(f"**Decades compared:** {results['decade_early'][0]}s vs "
                 f"{results['decade_late'][0]}s")
    lines.append(f"**Permutations:** {results['n_permutations']}; "
                 f"**bootstraps:** {results['n_bootstraps']}; "
                 f"**seed:** {results['random_seed']}\n")

    lines.append("## Aggregate Results\n")
    lines.append("| Statistic | Value |")
    lines.append("|-----------|-------|")
    lines.append(f"| Mean POS-share drift (full English) | {agg['mean_drift_all']:.4f} |")
    lines.append(f"| Mean POS-share drift (fiction only) | {agg['mean_drift_fiction']:.4f} |")
    lines.append(f"| Median drift (full) / (fiction) | {agg['median_drift_all']:.4f} / {agg['median_drift_fiction']:.4f} |")
    lines.append(f"| Mean paired diff (full − fiction) | {agg['mean_paired_diff']:+.4f} |")
    lines.append(f"| 95% bootstrap CI on paired diff | [{agg['paired_diff_ci_lo']:+.4f}, {agg['paired_diff_ci_hi']:+.4f}] |")
    lines.append(f"| Sign-flip permutation p on paired diff | {agg['paired_diff_p']:.4f} |")
    lines.append(f"| Spearman ρ (full vs fiction, across words) | {agg['spearman_rho_all_vs_fiction']:.3f} |")
    lines.append(f"| Permutation p on Spearman ρ | {agg['spearman_p_all_vs_fiction']:.4f} |")
    lines.append(f"| Words with permutation-significant drift (full) | {agg['n_sig_all']}/{agg['n_words']} |")
    lines.append(f"| Words with permutation-significant drift (fiction) | {agg['n_sig_fiction']}/{agg['n_words']} |")
    lines.append(f"| Words significant in BOTH corpora | {agg['n_sig_both']}/{agg['n_words']} |")
    lines.append(f"| Words FDR-significant at q ≤ 0.05 (full / fiction) | "
                 f"{agg['n_fdr_sig_all']}/{agg['n_words']} / {agg['n_fdr_sig_fiction']}/{agg['n_words']} |")
    lines.append(f"| Words block-permutation-significant at p ≤ 0.05 (full / fiction) | "
                 f"{agg['n_block_sig_all']}/{agg['n_words']} / {agg['n_block_sig_fiction']}/{agg['n_words']} |")
    lines.append(f"| Words whose fiction-drift lies within full-corpus 95% bootstrap CI | "
                 f"{agg['n_fiction_in_all_bootstrap_ci']}/{agg['n_words']} |")
    lines.append(f"| Trimmed mean paired diff (1 top + 1 bottom word dropped) | "
                 f"{agg['trimmed_mean_paired_diff']:+.4f} |")
    lines.append(f"| Trimmed mean drift (full / fiction) | "
                 f"{agg['trimmed_mean_drift_all']:.4f} / {agg['trimmed_mean_drift_fiction']:.4f} |")
    lines.append(f"| Mean frequency-trajectory Pearson r (full vs fiction) | "
                 f"{agg['mean_trajectory_r']:.3f} |")
    lines.append(f"| Words with trajectory r ≥ 0.7 (corpora agree) | "
                 f"{agg['n_trajectory_r_ge_0.7']}/{agg['n_words']} |")
    lines.append(f"| Words with trajectory r < 0.5 (corpora diverge) | "
                 f"{agg['n_trajectory_r_lt_0.5']}/{agg['n_words']} |\n")

    lines.append("## Per-Word Drift (early 1900s → late 1990s)\n")
    lines.append("| Word | Drift (full) | 95% CI (full) | p (full) | Drift (fiction) | 95% CI (fiction) | p (fiction) | Diff (full − fic) |")
    lines.append("|------|--------------|---------------|----------|-----------------|------------------|-------------|-------------------|")
    for w in results["target_words"]:
        a = results["per_word"][w][CORPUS_ALL]
        f_ = results["per_word"][w][CORPUS_GENRE]
        diff = a["drift"] - f_["drift"]
        lines.append(
            f"| {w} | {a['drift']:.4f} | "
            f"[{a['bootstrap_ci_lo']:.3f}, {a['bootstrap_ci_hi']:.3f}] | "
            f"{a['permutation_p']:.3f} | "
            f"{f_['drift']:.4f} | "
            f"[{f_['bootstrap_ci_lo']:.3f}, {f_['bootstrap_ci_hi']:.3f}] | "
            f"{f_['permutation_p']:.3f} | "
            f"{diff:+.4f} |"
        )

    lines.append("\n## Per-Word Frequency-Trajectory Pearson r (full vs fiction)\n")
    lines.append("| Word | Pearson r |")
    lines.append("|------|-----------|")
    for w in results["target_words"]:
        r = results["trajectory_correlations"][w]
        lines.append(f"| {w} | {r:.3f} |")

    lines.append("\n## Sensitivity to Decade Boundaries\n")
    lines.append("| Early decade | Late decade | Mean drift (full) | Mean drift (fic) | Mean diff | Spearman ρ | p(ρ) |")
    lines.append("|--------------|-------------|-------------------|------------------|-----------|------------|------|")
    for s in results["sensitivity"]:
        lines.append(
            f"| {s['early_decade'][0]}–{s['early_decade'][1]} | "
            f"{s['late_decade'][0]}–{s['late_decade'][1]} | "
            f"{s['mean_drift_all']:.4f} | "
            f"{s['mean_drift_fiction']:.4f} | "
            f"{s['mean_diff']:+.4f} | "
            f"{s['spearman_rho']:.3f} | "
            f"{s['spearman_p']:.3f} |"
        )

    lines.append("\n## Limitations\n")
    for lim in results["limitations"]:
        lines.append(f"- {lim}")

    report_path = os.path.join(workspace, REPORT_FILE)
    with open(report_path, "w") as f:
        f.write("\n".join(lines) + "\n")

    return results_path, report_path


# ---------- Verify ----------

def verify(results, workspace):
    checks_passed = 0
    checks_total = 0

    def check(name, cond):
        nonlocal checks_passed, checks_total
        checks_total += 1
        if cond:
            checks_passed += 1
            print(f"  PASS: {name}")
        else:
            print(f"  FAIL: {name}")

    agg = results["aggregate"]

    check("results.json exists",
          os.path.exists(os.path.join(workspace, RESULTS_FILE)))
    check("report.md exists",
          os.path.exists(os.path.join(workspace, REPORT_FILE)))
    check("Target word count == 20",
          len(results["target_words"]) == 20)
    check("Both corpora present",
          set(results["corpora"]) == {CORPUS_ALL, CORPUS_GENRE})
    check("POS tag count == 4",
          len(results["pos_tags"]) == 4)
    check("Aggregate mean_drift_all in [0, 1]",
          0.0 <= agg["mean_drift_all"] <= 1.0)
    check("Aggregate mean_drift_fiction in [0, 1]",
          0.0 <= agg["mean_drift_fiction"] <= 1.0)
    check("At least one word has permutation-significant drift in full corpus",
          agg["n_sig_all"] >= 1)
    check("Spearman rho between corpora in [-1, 1]",
          -1.0 <= agg["spearman_rho_all_vs_fiction"] <= 1.0)
    check("Paired diff CI brackets the point estimate",
          agg["paired_diff_ci_lo"] <= agg["mean_paired_diff"] <= agg["paired_diff_ci_hi"])
    check("Every word has both all-corpus and fiction-corpus records",
          all(CORPUS_ALL in results["per_word"][w] and CORPUS_GENRE in results["per_word"][w]
              for w in results["target_words"]))
    check("Bootstrap CIs are ordered lo <= hi for all (word, corpus) pairs",
          all(results["per_word"][w][c]["bootstrap_ci_lo"]
              <= results["per_word"][w][c]["bootstrap_ci_hi"]
              for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]))
    check("At least 4 specific limitations recorded",
          len(results["limitations"]) >= 4)
    check("Sensitivity analysis includes >=2 alternative decade pairs",
          len(results["sensitivity"]) >= 2)
    check("n_bootstraps and n_permutations are >= 1000",
          results["n_bootstraps"] >= 1000 and results["n_permutations"] >= 1000)

    check("corpus_sha256 present and 64 hex chars",
          isinstance(results.get("corpus_sha256"), str)
          and len(results["corpus_sha256"]) == 64)

    check("trajectory_correlations covers every target word",
          set(results.get("trajectory_correlations", {}).keys()) == set(results["target_words"]))

    check("All trajectory correlations in [-1, 1]",
          all(-1.0 <= v <= 1.0 for v in results.get("trajectory_correlations", {}).values()))

    check("At most half of trajectories are all-zero (corpus liveness)",
          results.get("n_total_trajectories", 0) > 0
          and results["n_all_zero_trajectories"] / results["n_total_trajectories"] <= 0.5)

    check("Benjamini-Hochberg q-values populated for every (word, corpus)",
          all("permutation_q_bh" in results["per_word"][w][c]
              for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]))

    check("All BH q-values are in [0, 1]",
          all(0.0 <= results["per_word"][w][c]["permutation_q_bh"] <= 1.0
              for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]))

    check("Block-permutation p-values populated for every (word, corpus)",
          all("block_permutation_p" in results["per_word"][w][c]
              for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]))

    # --- Additional rigor / sanity checks ---

    # (A) Row / trajectory count: 20 words × 2 corpora × 4 POS = 160 per-year trajectories;
    #     the manifest must record exactly that many fetches (or a multiple if rerun; here
    #     we assert the per-(word,corpus,POS) coverage is complete).
    check("Per-word trajectory coverage is complete (20 × 2 × 4 = 160)",
          sum(1 for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]) == 40
          and all("drift" in results["per_word"][w][c]
                  for w in results["target_words"] for c in [CORPUS_ALL, CORPUS_GENRE]))

    # (B) Effect-size plausibility: cosine distance between probability vectors is in [0, 2]
    #     and in practice ≤ 1. Cohen's-d-style standardized effect using paired diff SD.
    diffs = [results["per_word"][w][CORPUS_ALL]["drift"] - results["per_word"][w][CORPUS_GENRE]["drift"]
             for w in results["target_words"]]
    m = sum(diffs) / len(diffs)
    var = sum((d - m) ** 2 for d in diffs) / max(1, len(diffs) - 1)
    sd = math.sqrt(var) if var > 0 else 0.0
    cohens_d = abs(m) / sd if sd > 0 else 0.0
    check("Effect size |Cohen's d| < 5 on paired drift diff (sanity bound)",
          cohens_d < 5.0)

    # (C) Bootstrap CI width sanity: for each (word, corpus) with nonzero drift, CI width
    #     must be at least 1% of the drift estimate — a collapsed CI (lo == hi == drift)
    #     signals a broken bootstrap. Zero-drift cells (no POS change at all) are exempt.
    ci_widths_ok = True
    for w in results["target_words"]:
        for c in [CORPUS_ALL, CORPUS_GENRE]:
            drift = results["per_word"][w][c]["drift"]
            lo    = results["per_word"][w][c]["bootstrap_ci_lo"]
            hi    = results["per_word"][w][c]["bootstrap_ci_hi"]
            if drift > 1e-6 and (hi - lo) < max(1e-6, 0.01 * drift):
                ci_widths_ok = False
    check("Bootstrap CI width >= 1% of drift for non-zero-drift cells",
          ci_widths_ok)

    # (D) Sensitivity consistency: the sign of mean_drift_all - mean_drift_fiction under
    #     alternative decade pairs should not flip arbitrarily. We require that at least
    #     one sensitivity row matches the sign of the primary paired-diff, preventing
    #     the primary result from being a single-decade-pair fluke.
    primary_sign = 1 if agg["mean_paired_diff"] >= 0 else -1
    matching = 0
    for s in results.get("sensitivity", []):
        sign = 1 if s["mean_diff"] >= 0 else -1
        if sign == primary_sign:
            matching += 1
    check("At least one sensitivity-analysis decade pair agrees in sign with primary paired diff",
          matching >= 1)

    # (E) Negative / falsification control: synthetic stationary word should produce drift
    #     far smaller than aggregate real-word drift, and its permutation p-value should NOT
    #     be significant at p ≤ 0.05 (it should look null).
    nc = results.get("negative_control", {})
    check("Negative control drift is an order of magnitude smaller than mean real-word drift",
          nc.get("drift", 1.0) < max(agg["mean_drift_all"], agg["mean_drift_fiction"]) * 0.5 + 1e-6)
    check("Negative control permutation p-value > 0.05 (null is well-behaved)",
          nc.get("permutation_p", 0.0) > 0.05)

    # (F) Spearman-agreement plausibility: rho_all_vs_fiction shouldn't be exactly 0 or 1
    #     (which would indicate a rank tie or perfect reproduction). Anything strictly in
    #     (−1, 1) is plausible for real-world data.
    rho = agg["spearman_rho_all_vs_fiction"]
    check("Spearman ρ is strictly inside (−1, 1) (no degenerate tie)",
          -1.0 < rho < 1.0)

    # (G) Negative-control record contains required fields.
    check("Negative control record contains drift, p-value, and null stats",
          all(k in nc for k in ("drift", "permutation_p", "null_mean", "null_p95")))

    # (H) All-zero trajectory count is reported; a well-behaved run should
    #     have at most half of trajectories all-zero (stricter than the
    #     corpus-liveness check, which allows up to 50% before erroring).
    check("n_all_zero_trajectories is reported and non-negative",
          isinstance(results.get("n_all_zero_trajectories"), int)
          and results["n_all_zero_trajectories"] >= 0)

    # (I) At least 8 specific, distinct limitations are recorded (evaluator
    #     criterion: limitations should cover data, methodology, failure
    #     modes, and out-of-scope claims).
    check("At least 8 limitations documented (data, method, failure modes, scope)",
          len(results.get("limitations", [])) >= 8)

    # (J) Paired-diff bootstrap CI width should be much smaller than 1 and
    #     non-degenerate. For small 20-word samples the CI should span at
    #     least 1e-6.
    ci_lo = agg["paired_diff_ci_lo"]
    ci_hi = agg["paired_diff_ci_hi"]
    check("Paired-diff bootstrap CI is non-degenerate (hi - lo >= 1e-6) and < 2.0 in width",
          1e-6 <= (ci_hi - ci_lo) < 2.0)

    # (K) Determinism fingerprint: the corpus_sha256 should be a 64-character
    #     hex string derived from content (not timestamp) hashes — i.e., any
    #     two reruns on the same API payload produce an identical hash.
    import re as _re
    check("corpus_sha256 matches 64-char hex pattern (content-addressed)",
          isinstance(results.get("corpus_sha256"), str)
          and bool(_re.fullmatch(r"[0-9a-f]{64}", results["corpus_sha256"])))

    # (L) POS-share vectors for every (word, corpus) sum to ~1.0 (or exactly 0
    #     if no data) — a basic sanity check that share normalization is sound.
    shares_ok = True
    for w in results["target_words"]:
        for c in [CORPUS_ALL, CORPUS_GENRE]:
            for key in ("share_early", "share_late"):
                v = results["per_word"][w][c][key]
                s = sum(v)
                if not (abs(s - 1.0) < 1e-3 or s == 0.0):
                    shares_ok = False
    check("All POS-share vectors sum to 1.0 (or exactly 0 if empty)",
          shares_ok)

    # (M) BH q-values dominate raw p-values (q >= p for every cell).
    q_ge_p = True
    for w in results["target_words"]:
        for c in [CORPUS_ALL, CORPUS_GENRE]:
            if results["per_word"][w][c]["permutation_q_bh"] + 1e-9 < results["per_word"][w][c]["permutation_p"]:
                q_ge_p = False
    check("Benjamini-Hochberg q-values satisfy q >= p cell-wise (FDR correction is monotone)",
          q_ge_p)

    # (N) Sensitivity section covers different decade boundaries (not duplicates).
    decade_pairs = set()
    for s in results.get("sensitivity", []):
        decade_pairs.add((tuple(s["early_decade"]), tuple(s["late_decade"])))
    check("Sensitivity analysis uses distinct decade pairs (no duplicate rows)",
          len(decade_pairs) == len(results.get("sensitivity", [])))

    print(f"\n  Results: {checks_passed}/{checks_total} checks passed")
    return checks_passed == checks_total


# ---------- Main ----------

def main():
    workspace = SCRIPT_DIR
    if "--verify" in sys.argv:
        path = os.path.join(workspace, RESULTS_FILE)
        if not os.path.exists(path):
            print("FAIL: results.json not found. Run analysis first.")
            sys.exit(1)
        with open(path, "r") as f:
            results = json.load(f)
        print("[VERIFY] Running verification checks...")
        if verify(results, workspace):
            print("\nALL CHECKS PASSED")
        else:
            print("\nSOME CHECKS FAILED")
            sys.exit(1)
        return

    n_steps = 5
    print(f"[1/{n_steps}] Loading data from Google Books Ngrams "
          f"({len(TARGET_WORDS)} words × 2 corpora × {len(POS_TAGS)} POS tags "
          f"= {len(TARGET_WORDS) * 2 * len(POS_TAGS)} API calls)...")
    data, manifest = load_data()
    print(f"    cached queries: {len(manifest['fetches'])}")

    print(f"\n[2/{n_steps}] Running drift analysis, permutation null, and bootstrap CI...")
    results = run_analysis(data, manifest)
    agg = results["aggregate"]
    print(f"    mean drift (full English): {agg['mean_drift_all']:.4f}")
    print(f"    mean drift (fiction only): {agg['mean_drift_fiction']:.4f}")
    print(f"    mean paired diff:          {agg['mean_paired_diff']:+.4f} "
          f"(95% CI [{agg['paired_diff_ci_lo']:+.4f}, {agg['paired_diff_ci_hi']:+.4f}], "
          f"p = {agg['paired_diff_p']:.4f})")
    print(f"    Spearman ρ (full vs fic): {agg['spearman_rho_all_vs_fiction']:.3f} "
          f"(p = {agg['spearman_p_all_vs_fiction']:.4f})")
    print(f"    significant in full only: {agg['n_sig_all']}/{agg['n_words']}; "
          f"in fiction only: {agg['n_sig_fiction']}/{agg['n_words']}; "
          f"in both: {agg['n_sig_both']}/{agg['n_words']}")

    print(f"\n[3/{n_steps}] Writing results.json and report.md...")
    results_path, report_path = generate_report(results, workspace)
    # Persist manifest alongside.
    manifest_path = os.path.join(workspace, "data_manifest.json")
    with open(manifest_path, "w") as f:
        json.dump(manifest, f, indent=2)
    print(f"    wrote {results_path}")
    print(f"    wrote {report_path}")
    print(f"    wrote {manifest_path}")

    print(f"\n[4/{n_steps}] Sensitivity analysis (alternative decade boundaries):")
    for s in results["sensitivity"]:
        print(f"    {s['early_decade'][0]}s vs {s['late_decade'][0]}s: "
              f"drift_all={s['mean_drift_all']:.4f}, "
              f"drift_fic={s['mean_drift_fiction']:.4f}, "
              f"rho={s['spearman_rho']:.3f}")

    print(f"\n[5/{n_steps}] Summary")
    print(f"    target words:           {len(results['target_words'])}")
    print(f"    bootstrap iterations:   {results['n_bootstraps']}")
    print(f"    permutation iterations: {results['n_permutations']}")
    print(f"    random seed:            {results['random_seed']}")

    nc = results.get("negative_control", {})
    if nc:
        print(f"    negative control drift: {nc.get('drift', 0):.6f} "
              f"(p={nc.get('permutation_p', 0):.4f}); "
              f"mean real-word drift (full): {agg['mean_drift_all']:.4f}")

    print("\nANALYSIS COMPLETE")


if __name__ == "__main__":
    try:
        main()
    except KeyboardInterrupt:
        print("\nINTERRUPTED", file=sys.stderr)
        sys.exit(130)
    except Exception as e:
        print(f"\nFATAL: {type(e).__name__}: {e}", file=sys.stderr)
        print("The analysis did not complete. Common causes:", file=sys.stderr)
        print("  - No network access to books.google.com on first run", file=sys.stderr)
        print("  - Google renamed/removed the `en-2019` or `en-fiction-2019` corpora", file=sys.stderr)
        print("  - Insufficient disk space for the cache directory", file=sys.stderr)
        print("No results.json was written; safe to rerun (the cache of any successful "
              "fetches persists).", file=sys.stderr)
        sys.exit(1)
SCRIPT_EOF
```

**Expected output:** No output (script file created).

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_historical-embedding-semantic-drift-under-corpus-balanced-re && python3 script.py
```

**Expected output:**
- `[1/5]` progress messages with "fetched N/160 queries" counters
- `[2/5]` prints mean drift for each corpus, mean paired diff with 95% CI and p, Spearman ρ with p
- `[3/5]` confirms `results.json`, `report.md`, `data_manifest.json` written
- `[4/5]` prints sensitivity rows
- `[5/5]` summary counts
- Final line: `ANALYSIS COMPLETE`

First-time runtime: 6–12 minutes (≈160 API calls × ~0.35s delay + small retry overhead). Rerun from cache: <30s.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_historical-embedding-semantic-drift-under-corpus-balanced-re && python3 script.py --verify
```

**Expected output:**
- 37+ verification checks, each printed as `PASS: <name>` or `FAIL: <name>`
- A final `Results: N/N checks passed` line
- Final line: `ALL CHECKS PASSED`
- Exit code 0 on pass, 1 on any failure

---

## Determinism Guarantees

Every operation that could introduce run-to-run variance is either seeded or removed:

- **Random number generator.** A single `random.Random(RANDOM_SEED=42)` instance threads through all shuffles and bootstraps. The negative-control uses a second seeded RNG (`seed + 1`) so it is independent of, but deterministic alongside, the primary analysis.
- **Dictionary / set iteration.** All iteration proceeds over the `TARGET_WORDS` list and the fixed `[CORPUS_ALL, CORPUS_GENRE]` list; no iteration relies on `dict.keys()` or `set` ordering for mathematically material values.
- **Content-addressed corpus fingerprint.** `corpus_sha256` hashes the response time-series bytes (via `json.dumps(ts, sort_keys=True)`), not the cache file (which carries a `fetched_at` timestamp). Two reruns on the same API payloads produce an identical `corpus_sha256`.
- **No floating-point reductions in unspecified order.** All reductions iterate lists in index order.
- **No OS-dependent calls on the hot path.** No `os.listdir` without sorting; no file-mtime comparisons.

Given an identical cache state, rerunning the script produces byte-identical `results.json` (modulo the `fetched_at` timestamps persisted only in the cache).

## Success Criteria

A run is considered successful if and only if **all** of the following machine-checkable conditions hold (the `--verify` command asserts each one):

1. Script runs to completion and emits the literal final line `ANALYSIS COMPLETE`.
2. `results.json` exists with both `aggregate` and `per_word` sections, and a `negative_control` record.
3. `report.md` exists with four tables (aggregate, per-word drift, per-word frequency-trajectory Pearson r, sensitivity) and a limitations list of ≥ 4 entries.
4. All 37+ verification assertions pass under `python3 script.py --verify`; the final line is `ALL CHECKS PASSED`.
5. Per-word coverage is complete: every (word ∈ TARGET_WORDS, corpus ∈ {CORPUS_ALL, CORPUS_GENRE}) pair has a `drift`, `bootstrap_ci_lo`, `bootstrap_ci_hi`, `permutation_p`, `block_permutation_p`, and `permutation_q_bh`.
6. At least one target word has permutation-significant drift (`permutation_p ≤ 0.05`) in the full-English corpus.
7. The aggregate Spearman ρ between full-corpus and fiction-corpus per-word drift lies in (−1, 1).
8. Every per-word bootstrap CI is well-ordered (`lo ≤ hi`), non-degenerate for nonzero-drift cells (CI width ≥ 1% of drift), and brackets the paired-diff point estimate.
9. The synthetic negative-control word's observed drift is < 50% of the mean real-word drift AND its permutation p-value is > 0.05 (the machinery does not invent drift from stationary inputs).
10. The key paired-difference effect size has |Cohen's d| < 5 (no unphysically large standardized effect).
11. At least one sensitivity-analysis decade pair agrees in sign with the primary paired diff (the finding is not a single-decade-pair artifact).

## Failure Conditions

The analysis is considered to have **failed** (rather than succeeded with a null result) when any of these hold:

1. **Network / API failure** — If `books.google.com/ngrams/json` returns non-JSON (e.g., an HTML CAPTCHA page) or repeated HTTP errors after 4 backoff retries on any single query, `fetch_ngram` raises; the top-level handler in `main()` prints a diagnostic to `stderr`, no `results.json` is written, and the process exits with status 1. Rerun later — successfully cached queries persist.
2. **Corpus identifier deprecated** — If Google renames or removes `en-2019` / `en-fiction-2019`, the API returns empty trajectories. A corpus-liveness check in `load_data()` raises if more than 50% of trajectories are all-zero; update `CORPUS_ALL` / `CORPUS_GENRE` in the `DOMAIN CONFIGURATION` block to a current identifier (e.g., `en-2020`).
3. **Cache corruption** — If a cached JSON file is truncated, the length check at fetch time detects the mismatch and the entry is re-fetched.
4. **Verification failure** — If `python3 script.py --verify` prints a `FAIL:` line, exits with status 1, or fails to print `ALL CHECKS PASSED`, the run is invalid and the result numbers should not be reported.
5. **Negative-control failure** — If the synthetic stationary word produces drift that is not clearly smaller than real-word drift, or its permutation p ≤ 0.05, the test machinery is over-powered or biased and the primary p-values should not be trusted.
6. **Single-decade fragility** — If the sign of the paired-diff mean flips under both alternative decade boundaries in the `sensitivity` section, the headline finding is an artifact of the specific 1900s/1990s framing.

## Limitations

1. **Instrument coarseness.** Drift is operationalized as change in a 4-dimensional POS-share vector, a coarse-grained proxy for semantic drift; it does not detect meaning shifts that leave the POS distribution invariant (e.g., a noun whose referent changes without category change). This is a known limitation of POS-share analyses and is *not* a substitute for full diachronic embedding analysis.
2. **Corpus provenance.** The Google Books corpus has known biases: OCR errors, scanning-coverage changes over time, overrepresentation of academic material after 1950, and genre labels that are not vetted by linguists. We cannot correct these confounders inside a single-genre restriction.
3. **Single-genre proxy.** The `en-fiction-2019` subcorpus is a proxy for "a single genre"; fiction itself is heterogeneous (literary vs pulp vs genre fiction) and its internal composition also shifts across the 20th century.
4. **Autocorrelation.** The year-shuffling permutation null tests a null of no temporal ordering; it does *not* control for year-to-year autocorrelation driven by language-independent publishing trends. The block-permutation null is a partial remedy but assumes a 10-year block length.
5. **Tagger noise.** POS tagging in the Google Ngrams corpus is produced by an automatic tagger whose error rate is non-zero; rare POS assignments for a word may reflect tagger noise rather than usage.
6. **Target-word selection bias.** Target words are drawn from a specific pair of published papers (Hamilton et al. 2016, Kulkarni et al. 2015); conclusions may not generalize to unannotated words with subtler drift.
7. **Scope of claim.** The analysis does **not** show that the underlying papers' embedding-based drift scores are wrong; it only shows whether POS-share drift, a cheaper instrument, is corpus-robust.
8. **Language scope.** The corpus is restricted to English; drift dynamics in non-English corpora may differ and would require rerunning with changed `CORPUS_ALL`/`CORPUS_GENRE` identifiers.
9. **Low-frequency words.** Words with very low overall frequency will have unstable POS-share estimates; we do not filter by frequency.
10. **Pre-existence artifacts.** Words like `browser`, `server`, `net` have essentially no attested usage before 1960, so their early-decade POS-share vectors are dominated by the uniform fallback in `normalize_share`, which can inflate apparent drift.
11. **Negative-control realism.** The synthetic negative-control word has small Gaussian jitter on POS shares and so produces non-zero drift by construction; we require only that it be an order of magnitude smaller than real-word drift, not literally zero.

## Generalizability

The pipeline is a general framework for "categorical-share drift under a subcorpus restriction":

- Replace `TARGET_WORDS` with any list of terms of interest and `CORPUS_GENRE` with any single-category subcorpus available via the Ngrams API.
- Replace `POS_TAGS` with any mutually exclusive categorical labeling: polarity (`_POSITIVE`, `_NEGATIVE`), syntactic roles, or any fine-grained tag set supported by Google's POS tagger.
- The statistical core — cosine distance between normalized share vectors across two time windows, with decade-shuffling permutation nulls, paired-bootstrap CIs, and Spearman rank agreement between corpora — transfers verbatim to problems in:
  - Citation-category drift across publication years
  - Industry mix in employment statistics across decades
  - Product-category mix in retail transaction logs across years
  - Genre mix in streaming-platform catalogs across years