#!/usr/bin/env python3
"""Recompute "the one letter" from the published vectors alone.

Input: adams-jefferson-17d.npz (vectors_17d: 1103 x 17 float64; letter_id per row;
schema_names). No letter text, no embeddings. numpy + standard library only.

What this reproduces: the declared permutation test (RIFF_two-clocks.md, DECLARED
RUN 1, Read 3). For each of 102 reading orders (true, reversed, 100 random with
seed 5) the cumulative covariance's eigenvalue trajectory is classified by
habitat's own jump rule, reproduced here verbatim from
src/habitat/core/field/overlay_eventuality.py (classify_eventualities, the
ACHIEVEMENT branch), so that this file has no dependency on the habitat package.

Output: per letter, the fraction of random orders in which it hosts a jump onset,
its status in the true/reversed orders, and the size-expected fraction under the
null. Same seed, same bytes, same numbers, on any machine.

Run:  python recompute_one_letter.py adams-jefferson-17d.npz
"""
from __future__ import annotations

import sys
from collections import defaultdict

import numpy as np

K = 100
SEED = 5

# ── habitat's jump rule, verbatim in effect (overlay_eventuality.classify_eventualities) ──
_REGIME_MIN = {"gradient": 1e-4, "compressed": 1e-5, "categorical": 1e-6}


def _regime_for_dim(d: int) -> str:
    return "gradient" if d <= 4 else "compressed" if d <= 9 else "categorical"


def achievement_dims(E: np.ndarray) -> set:
    """Dimensions (columns) classified ACHIEVEMENT on the trajectory matrix E (n_steps x n_dims)."""
    if E.ndim != 2 or E.shape[0] < 3:
        return set()
    n_steps, n_dims = E.shape
    out = set()
    for d in range(n_dims):
        traj = E[:, d]
        regime = _regime_for_dim(d)
        dim_std = float(np.std(traj))
        diffs = np.diff(traj)
        max_step = float(np.max(np.abs(diffs))) if len(diffs) else 0.0
        regime_min = _REGIME_MIN[regime]
        jump_threshold = dim_std * (1.5 + 1.0 / np.sqrt(n_steps)) if dim_std > 1e-6 else regime_min * 10
        if len(diffs) > 0 and max_step > jump_threshold:
            last_step = abs(float(diffs[-1]))
            if last_step > jump_threshold * 0.8:
                out.add(d)
    return out


def onset_steps(V: np.ndarray) -> set:
    N = V.shape[0]
    E = np.zeros((N, V.shape[1]))
    for t in range(2, N):
        E[t] = np.linalg.eigvalsh(np.cov(V[: t + 1], rowvar=False))
    steps, prev = set(), set()
    for t in range(3, N + 1):
        now = achievement_dims(E[:t])
        if now - prev:
            steps.add(t)
        prev = now
    return steps


def main(path: str) -> None:
    z = np.load(path, allow_pickle=True)
    V = np.asarray(z["vectors_17d"], dtype=np.float64)
    letters = np.array([str(s) for s in z["letter_id"]])
    N = V.shape[0]
    uniq = sorted(set(letters))
    rng = np.random.default_rng(SEED)
    orders = {"true": np.arange(N), "reversed": np.arange(N)[::-1]}
    for i in range(K):
        orders[f"rand{i}"] = rng.permutation(N)
    hosts = defaultdict(lambda: defaultdict(int))
    n_onsets = {}
    for name, perm in orders.items():
        st = onset_steps(V[perm])
        n_onsets[name] = len(st)
        for t in st:
            hosts[letters[perm[t - 1]]][name] += 1
    rand = [n for n in orders if n.startswith("rand")]
    base = float(np.mean([n_onsets[n] for n in rand])) / N
    print(f"{N} compositions, {len(uniq)} letters, {len(orders)} orders (seed {SEED})")
    print(f"onsets: true={n_onsets['true']} reversed={n_onsets['reversed']} random mean={np.mean([n_onsets[n] for n in rand]):.1f}")
    rows = []
    for L in uniq:
        frac = float(np.mean([hosts[L][n] > 0 for n in rand]))
        fw, rv = hosts[L]["true"] > 0, hosts[L]["reversed"] > 0
        status = "both" if fw and rv else "forward-only" if fw else "reverse-only" if rv else "neither"
        k = int((letters == L).sum())
        rows.append((L, k, frac, 1 - (1 - base) ** k, status))
    print("\nLetters jumping in BOTH the true and the reversed order:")
    for L, k, frac, exp, status in rows:
        if status == "both":
            print(f"  {L.split('/')[-1]:34s} k={k:3d}  random-order fraction={frac:.2f}  size-expected={exp:.2f}")
    print("\nTop 10 by random-order fraction over size-expected:")
    for L, k, frac, exp, status in sorted(rows, key=lambda r: -(r[2] - r[3]))[:10]:
        print(f"  {L.split('/')[-1]:34s} k={k:3d}  frac={frac:.2f}  expected={exp:.2f}  {status}")


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else "adams-jefferson-17d.npz")