"""Inductive methods. Each takes pages (texts.py format) and returns (metrics, artifacts):
metrics are scalars comparable across the manuscript and the null texts; artifacts are tables for
inspection (topics, clusters, break points...).  Nothing here tests a hypothesis: these look for
structure, and what they find is turned into hypotheses afterwards.
"""
import os, sys
_here = os.path.dirname(os.path.abspath(__file__))  # data/analysis, or a subfolder of it
sys.path.insert(0, _here if os.path.exists(os.path.join(_here, "paths.py")) else os.path.dirname(_here))
from paths import *  # noqa: repository locations
import collections, math, zlib, random
import numpy as np
import pandas as pd
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.decomposition import LatentDirichletAllocation, TruncatedSVD
from sklearn.metrics import normalized_mutual_info_score as nmi
from sklearn.cluster import KMeans, AgglomerativeClustering
from texts import page_words
import sys, os
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
from corpus import tokenize_glyphs

META = ["section", "lang", "hand", "quire"]


def _docs(pages): return [" ".join(page_words(p)) for p in pages]


def _pair_z(sim, pages, rel, n_perm=500, seed=0):
    """mean similarity of page pairs related by `rel` vs random pairs drawn from the same two sections (z)"""
    idx = [(i, j) for i in range(len(pages)) for j in range(i + 1, len(pages)) if rel(pages[i], pages[j])]
    if not idx: return None
    obs = np.mean([sim[i, j] for i, j in idx]); rng = random.Random(seed); bysec = collections.defaultdict(list)
    for k, p in enumerate(pages): bysec[p["section"]].append(k)
    null = []
    for _ in range(n_perm):
        v = []
        for i, j in idx:
            a = rng.choice(bysec[pages[i]["section"]]); b = rng.choice(bysec[pages[j]["section"]])
            while b == a: b = rng.choice(bysec[pages[j]["section"]])
            v.append(sim[a, b])
        null.append(np.mean(v))
    return round(float((obs - np.mean(null)) / (np.std(null) + 1e-12)), 2)


# ---------------------------------------------------------------- 1. topic model
K_FIXED = 10  # the perplexity-best k is reported; the topic tables use a fixed 10 for detail


def topics(pages, k_grid=(4, 6, 8, 12, 16), seed=0):
    docs = _docs(pages); cv = CountVectorizer(token_pattern=r"\S+", min_df=3); X = cv.fit_transform(docs)
    rng = np.random.RandomState(seed); te = rng.rand(len(docs)) < 0.2
    best = None
    for k in k_grid:
        m = LatentDirichletAllocation(n_components=k, random_state=seed, learning_method="batch", max_iter=60).fit(X[~te])
        pp = m.perplexity(X[te])
        if best is None or pp < best[0]: best = (pp, k)
    k = K_FIXED or best[1]; m = LatentDirichletAllocation(n_components=k, random_state=seed, learning_method="batch", max_iter=100).fit(X)
    th = m.transform(X); top = th.argmax(1); vocab = np.array(cv.get_feature_names_out())
    met = dict(k=k, k_best_perplexity=best[1], heldout_perplexity=round(best[0], 1))
    for f in META: met[f"nmi_{f}"] = round(nmi([p[f] for p in pages], top), 3)
    sim = th @ th.T / (np.linalg.norm(th, axis=1)[:, None] * np.linalg.norm(th, axis=1)[None, :] + 1e-12)
    met["bifolio_z"] = _pair_z(sim, pages, lambda a, b: a["bifolio"] == b["bifolio"] and a["bifolio"] is not None)
    met["mean_max_theta"] = round(float(th.max(1).mean()), 3)
    tw = [dict(topic=t, top_words=" ".join(vocab[np.argsort(-m.components_[t])[:12]]),
               pages=int((top == t).sum()), sections=dict(collections.Counter(p["section"] for p, tt in zip(pages, top) if tt == t).most_common())) for t in range(k)]
    pt = pd.DataFrame(th, columns=[f"t{t}" for t in range(k)]); pt.insert(0, "page", [p["page"] for p in pages]); pt.insert(1, "section", [p["section"] for p in pages])
    pt.insert(2, "lang", [p["lang"] for p in pages]); pt.insert(3, "top", top)
    return met, dict(topics=pd.DataFrame(tw), page_topics=pt)


# ---------------------------------------------------------------- 2. distributional word classes (PPMI-SVD embeddings)
def word_classes(pages, min_count=8, dim=40, k=12, seed=0):
    seqs = [l for p in pages for l in p["lines"]]; cnt = collections.Counter(w for s in seqs for w in s)
    vocab = [w for w, c in cnt.most_common() if c >= min_count]; vi = {w: i for i, w in enumerate(vocab)}; V = len(vocab)
    C = np.zeros((V, 2 * V))
    for s in seqs:
        for i, w in enumerate(s):
            if w not in vi: continue
            if i > 0 and s[i - 1] in vi: C[vi[w], vi[s[i - 1]]] += 1
            if i + 1 < len(s) and s[i + 1] in vi: C[vi[w], V + vi[s[i + 1]]] += 1
    tot = C.sum(); pr = C.sum(1, keepdims=True) / tot; pc = C.sum(0, keepdims=True) / tot
    with np.errstate(divide="ignore", invalid="ignore"): pmi = np.log((C / tot) / (pr * pc))
    ppmi = np.nan_to_num(np.maximum(pmi, 0)); E = TruncatedSVD(n_components=min(dim, V - 1), random_state=seed).fit_transform(ppmi)
    E /= np.linalg.norm(E, axis=1, keepdims=True) + 1e-12
    lab = KMeans(n_clusters=k, random_state=seed, n_init=10).fit_predict(E)
    first = [w[0] for w in vocab]; last = [w[-1] for w in vocab]; ln = [min(len(w), 8) for w in vocab]
    met = dict(vocab=V, nmi_class_first_letter=round(nmi(first, lab), 3), nmi_class_last_letter=round(nmi(last, lab), 3), nmi_class_length=round(nmi(ln, lab), 3))
    # do words of a class sit at particular places in the line?
    pos = collections.defaultdict(collections.Counter)
    for s in seqs:
        for i, w in enumerate(s):
            if w in vi: pos[lab[vi[w]]]["first" if i == 0 else "last" if i == len(s) - 1 else "inside"] += 1
    met["nmi_class_line_position"] = round(_mi_table(pos), 4)
    cl = [dict(cls=c, n=int((lab == c).sum()), examples=" ".join([vocab[i] for i in np.where(lab == c)[0][:14]]),
               line_first=round(pos[c]["first"] / max(1, sum(pos[c].values())), 3), line_last=round(pos[c]["last"] / max(1, sum(pos[c].values())), 3)) for c in range(k)]
    return met, dict(classes=pd.DataFrame(cl))


def _mi_table(tab):
    rows = list(tab); cols = sorted({c for r in rows for c in tab[r]}); M = np.array([[tab[r][c] for c in cols] for r in rows], float); N = M.sum()
    p = M / N; pr = p.sum(1, keepdims=True); pc = p.sum(0, keepdims=True)
    with np.errstate(divide="ignore", invalid="ignore"): t = np.where(p > 0, p * np.log2(p / (pr * pc)), 0)
    return float(t.sum())


# ---------------------------------------------------------------- 3. unsupervised morphology (Morfessor)
def morphology(pages, seed=0):
    import morfessor
    cnt = collections.Counter(w for p in pages for w in page_words(p))
    io = morfessor.MorfessorIO(); m = morfessor.BaselineModel()
    m.load_data([(c, tuple(tokenize_glyphs(w))) for w, c in cnt.items()], count_modifier=lambda x: int(round(1 + math.log(x))))
    random.seed(seed); m.train_batch()
    segs = {w: m.viterbi_segment(tuple(tokenize_glyphs(w)))[0] for w in cnt}
    morphs = collections.Counter(); stems = collections.defaultdict(set); n_morph = []
    for w, s in segs.items():
        n_morph.append(len(s)); ms = ["".join(x) for x in s]
        for x in ms: morphs[x] += cnt[w]
        if len(ms) >= 2: stems["".join(ms[:-1])].add(ms[-1])
    fam = [len(v) for v in stems.values()]
    met = dict(types=len(cnt), morph_types=len(morphs), morphs_per_word=round(float(np.mean(n_morph)), 2),
               paradigm_mean_endings=round(float(np.mean(fam)), 2) if fam else 0, share_stems_with_3plus_endings=round(float(np.mean([f >= 3 for f in fam])), 3) if fam else 0)
    top = pd.DataFrame([dict(morph=m_, tokens=c) for m_, c in morphs.most_common(40)])
    return met, dict(morphs=top)


# ---------------------------------------------------------------- 4. change points along the book
def change_points(pages, n_top=200, pen=None):
    import ruptures as rpt
    cnt = collections.Counter(w for p in pages for w in page_words(p)); top = [w for w, c in cnt.most_common(n_top)]; ti = {w: i for i, w in enumerate(top)}
    X = np.zeros((len(pages), n_top))
    for k, p in enumerate(pages):
        for w in page_words(p):
            if w in ti: X[k, ti[w]] += 1
    X = X / (X.sum(1, keepdims=True) + 1e-9); X = np.sqrt(X)  # Hellinger
    algo = rpt.Pelt(model="l2", min_size=2).fit(X); bk = algo.predict(pen=pen or 0.8)
    bk = [b for b in bk if b < len(pages)]
    sec_bounds = {k for k in range(1, len(pages)) if pages[k]["section"] != pages[k - 1]["section"]}
    lang_bounds = {k for k in range(1, len(pages)) if pages[k]["lang"] != pages[k - 1]["lang"]}
    hand_bounds = {k for k in range(1, len(pages)) if pages[k]["hand"] != pages[k - 1]["hand"]}
    near = lambda b, S: any(abs(b - s) <= 1 for s in S)
    met = dict(breaks=len(bk), share_at_section_boundary=round(np.mean([near(b, sec_bounds) for b in bk]), 3) if bk else 0,
               share_at_language_boundary=round(np.mean([near(b, lang_bounds) for b in bk]), 3) if bk else 0,
               share_at_hand_boundary=round(np.mean([near(b, hand_bounds) for b in bk]), 3) if bk else 0,
               section_boundaries_found=round(np.mean([near(s, set(bk)) for s in sec_bounds]), 3) if sec_bounds else 0)
    tab = pd.DataFrame([dict(break_before=pages[b]["page"], after=pages[b - 1]["page"], section=f"{pages[b-1]['section']}>{pages[b]['section']}",
                             lang=f"{pages[b-1]['lang']}>{pages[b]['lang']}", hand=f"{pages[b-1]['hand']}>{pages[b]['hand']}") for b in bk])
    return met, dict(breaks=tab)


# ---------------------------------------------------------------- 5. seriation (correspondence analysis)
def seriation(pages, n_top=400):
    cnt = collections.Counter(w for p in pages for w in page_words(p)); top = [w for w, c in cnt.most_common(n_top)]; ti = {w: i for i, w in enumerate(top)}
    X = np.zeros((len(pages), n_top))
    for k, p in enumerate(pages):
        for w in page_words(p):
            if w in ti: X[k, ti[w]] += 1
    P = X / X.sum(); r = P.sum(1, keepdims=True); c = P.sum(0, keepdims=True)
    S = (P - r @ c) / np.sqrt(r @ c + 1e-12); U, s, Vt = np.linalg.svd(S, full_matrices=False)
    f1 = U[:, 0] * s[0] / np.sqrt(r[:, 0]); f2 = U[:, 1] * s[1] / np.sqrt(r[:, 0])
    inertia = s ** 2 / (s ** 2).sum()
    isB = np.array([p["lang"] == "B" for p in pages], float)
    qn = np.array([ord(str(p["quire"])[0]) if p["quire"] else 0 for p in pages], float)
    from scipy.stats import spearmanr
    met = dict(inertia_axis1=round(float(inertia[0]), 3), inertia_axis2=round(float(inertia[1]), 3),
               axis1_vs_langB=round(float(abs(spearmanr(f1, isB)[0])), 3), axis1_vs_quire=round(float(abs(spearmanr(f1, qn)[0])), 3),
               axis2_vs_langB=round(float(abs(spearmanr(f2, isB)[0])), 3))
    tab = pd.DataFrame(dict(page=[p["page"] for p in pages], section=[p["section"] for p in pages], lang=[p["lang"] for p in pages], hand=[p["hand"] for p in pages],
                            quire=[p["quire"] for p in pages], axis1=np.round(f1, 4), axis2=np.round(f2, 4)))
    return met, dict(coordinates=tab)


# ---------------------------------------------------------------- 6. compression distance clustering
def ncd(pages, k=7):
    docs = [" ".join(page_words(p)).encode() for p in pages]; C = [len(zlib.compress(d, 9)) for d in docs]; n = len(docs)
    D = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            cij = len(zlib.compress(docs[i] + b" " + docs[j], 9)); D[i, j] = D[j, i] = (cij - min(C[i], C[j])) / max(C[i], C[j])
    lab = AgglomerativeClustering(n_clusters=k, metric="precomputed", linkage="average").fit_predict(D)
    met = {f"nmi_{f}": round(nmi([p[f] for p in pages], lab), 3) for f in META}
    met["bifolio_z"] = _pair_z(1 - D, pages, lambda a, b: a["bifolio"] == b["bifolio"] and a["bifolio"] is not None)
    return met, {}


# ---------------------------------------------------------------- 7. co-occurrence communities
def communities(pages, min_count=10, seed=0):
    import networkx as nx
    from networkx.algorithms.community import louvain_communities
    cnt = collections.Counter(w for p in pages for w in page_words(p)); vocab = {w for w, c in cnt.items() if c >= min_count}
    co = collections.Counter(); N = 0
    for p in pages:
        for l in p["lines"]:
            ws = [w for w in set(l) if w in vocab]; N += 1
            for a in range(len(ws)):
                for b in range(a + 1, len(ws)): co[tuple(sorted((ws[a], ws[b])))] += 1
    lc = collections.Counter(w for p in pages for l in p["lines"] for w in set(l) if w in vocab)
    G = nx.Graph()
    for (a, b), c in co.items():
        pmi = math.log(c * N / (lc[a] * lc[b]))
        if c >= 3 and pmi > 0.5: G.add_edge(a, b, weight=pmi)
    if G.number_of_nodes() < 5: return dict(nodes=G.number_of_nodes()), {}
    comms = louvain_communities(G, weight="weight", seed=seed)
    # section concentration of each community: share of its tokens in its commonest section
    sec_of = collections.defaultdict(collections.Counter)
    for p in pages:
        for w in page_words(p): sec_of[w][p["section"]] += 1
    conc = []
    for cm in comms:
        tot = collections.Counter()
        for w in cm: tot.update(sec_of[w])
        conc.append(max(tot.values()) / sum(tot.values()))
    from networkx.algorithms.community.quality import modularity
    met = dict(nodes=G.number_of_nodes(), edges=G.number_of_edges(), communities=len(comms), modularity=round(modularity(G, comms, weight="weight"), 3),
               mean_section_concentration=round(float(np.mean(conc)), 3))
    tab = pd.DataFrame([dict(community=i, size=len(c), concentration=round(conc[i], 3), examples=" ".join(sorted(c, key=lambda w: -cnt[w])[:12])) for i, c in enumerate(comms)]).sort_values("size", ascending=False)
    return met, dict(communities=tab)


# ---------------------------------------------------------------- 8. vocabulary richness and the unseen
def richness(pages):
    cnt = collections.Counter(w for p in pages for w in page_words(p)); f = collections.Counter(cnt.values()); n = sum(cnt.values()); S = len(cnt)
    f1, f2 = f[1], f[2]; chao1 = S + (f1 * f1) / (2 * f2) if f2 else S
    met = dict(tokens=n, types=S, singletons=f1, doubletons=f2, chao1=round(chao1), unseen_share=round(1 - S / chao1, 3), good_turing_missing_mass=round(f1 / n, 3))
    # types expected in 14 more leaves (about 14/102 more text), Good-Toulmin
    t = 14 / 102; gt = sum(((-1) ** (i + 1)) * (t ** i) * f[i] for i in range(1, 30)); met["new_types_in_lost_leaves"] = round(gt)
    return met, {}


# ---------------------------------------------------------------- 9. layout layer: what the line position changes
def layout(pages):
    rows = []
    for p in pages:
        for l in p["lines"]:
            if len(l) < 3: continue
            for i, w in enumerate(l):
                g = tokenize_glyphs(w)
                rows.append(dict(pos="first" if i == 0 else "last" if i == len(l) - 1 else "inside", first=g[0], last=g[-1], length=min(len(g), 9)))
    df = pd.DataFrame(rows); met = {}
    for f in ("first", "last", "length"):
        tab = collections.defaultdict(collections.Counter)
        for a, b in zip(df["pos"], df[f]): tab[a][b] += 1
        met[f"mi_line_position_{f}"] = round(_mi_table(tab), 4)
    met["mean_len_first"] = round(float(df[df.pos == "first"].length.mean()), 2); met["mean_len_inside"] = round(float(df[df.pos == "inside"].length.mean()), 2)
    met["mean_len_last"] = round(float(df[df.pos == "last"].length.mean()), 2)
    return met, {}


# ---------------------------------------------------------------- 10. text and picture: vocabulary against pigment
def pigment(pages, paint_csv, seed=0, n_perm=300):
    pm = pd.read_csv(paint_csv); pm["page"] = "f" + pm["folio"].astype(str)
    pm = pm.set_index("page"); keep = [p for p in pages if p["page"] in pm.index]
    if len(keep) < 30: return dict(pages=len(keep)), {}
    cnt = collections.Counter(w for p in keep for w in page_words(p)); top = [w for w, c in cnt.most_common(300) if c >= 20]
    F = np.array([[collections.Counter(page_words(p))[w] / max(1, len(page_words(p))) for w in top] for p in keep])
    out = {}; tabs = []
    from scipy.stats import spearmanr
    rng = np.random.RandomState(seed); secs = np.array([p["section"] for p in keep])
    for col in ("green", "blue", "red"):
        y = pm.loc[[p["page"] for p in keep], col].values
        r = np.array([spearmanr(F[:, j], y)[0] if F[:, j].std() > 0 else 0 for j in range(len(top))]); r = np.nan_to_num(r)
        obs = np.sort(np.abs(r))[-10:].mean()
        null = []
        for _ in range(n_perm):  # permute pigment within section, so section vocabulary cannot drive it
            yp = y.copy()
            for s in set(secs): ix = np.where(secs == s)[0]; yp[ix] = rng.permutation(yp[ix])
            rp = np.nan_to_num(np.array([spearmanr(F[:, j], yp)[0] if F[:, j].std() > 0 else 0 for j in range(len(top))]))
            null.append(np.sort(np.abs(rp))[-10:].mean())
        out[f"{col}_top10_abs_rho"] = round(float(obs), 3); out[f"{col}_z_within_section"] = round(float((obs - np.mean(null)) / (np.std(null) + 1e-12)), 2)
        order = np.argsort(-np.abs(r))[:8]; tabs += [dict(pigment=col, word=top[j], rho=round(float(r[j]), 3)) for j in order]
    out["pages"] = len(keep)
    return out, dict(words=pd.DataFrame(tabs))


METHODS = dict(topics=topics, word_classes=word_classes, morphology=morphology, change_points=change_points, seriation=seriation,
               ncd=ncd, communities=communities, richness=richness, layout=layout)


# ---------------------------------------------------------------- 11. stylometry (Burrows' Delta)
def stylometry(pages, n_words=150, min_len=80, k=5):
    ps = [p for p in pages if len(page_words(p)) >= min_len]
    cnt = collections.Counter(w for p in ps for w in page_words(p)); top = [w for w, c in cnt.most_common(n_words)]
    F = np.array([[collections.Counter(page_words(p))[w] / len(page_words(p)) for w in top] for p in ps])
    Z = (F - F.mean(0)) / (F.std(0) + 1e-12); D = np.abs(Z[:, None, :] - Z[None, :, :]).mean(2)
    np.fill_diagonal(D, np.inf); nn = D.argmin(1)
    met = {f"nn_same_{f}": round(float(np.mean([ps[i][f] == ps[j][f] for i, j in enumerate(nn)])), 3) for f in META}
    np.fill_diagonal(D, 0); lab = AgglomerativeClustering(n_clusters=k, metric="precomputed", linkage="average").fit_predict(D)
    met.update({f"nmi_{f}": round(nmi([p[f] for p in ps], lab), 3) for f in META}); met["pages"] = len(ps)
    return met, {}


# ---------------------------------------------------------------- 12. anomalous pages (isolation forest)
def anomalies(pages, seed=0):
    from sklearn.ensemble import IsolationForest
    units = sorted({g for p in pages for w in page_words(p) for g in tokenize_glyphs(w)})
    ui = {u: i for i, u in enumerate(units)}; X = []
    for p in pages:
        ws = page_words(p); g = np.zeros(len(units))
        for w in ws:
            for u in tokenize_glyphs(w): g[ui[u]] += 1
        g /= g.sum() + 1e-9
        X.append(list(g) + [np.mean([len(tokenize_glyphs(w)) for w in ws]), len(set(ws)) / len(ws)])
    X = np.array(X); s = -IsolationForest(random_state=seed, n_estimators=300).fit(X).score_samples(X)
    order = np.argsort(-s)[:12]
    tab = pd.DataFrame([dict(page=pages[i]["page"], section=pages[i]["section"], lang=pages[i]["lang"], hand=pages[i]["hand"], words=len(page_words(pages[i])), score=round(float(s[i]), 3)) for i in order])
    return dict(top_score=round(float(s[order[0]]), 3), median_score=round(float(np.median(s)), 3)), dict(top=tab)


# ---------------------------------------------------------------- 13. paradigm structure (Kober-style stem x ending grid)
ENDINGS = ["aiin", "iin", "ain", "in", "dy", "ey", "y", "ol", "al", "or", "ar", "l", "r", "s", "m", "n", "o", "d"]


def paradigms(pages, min_stem_tokens=20):
    cnt = collections.Counter(w for p in pages for w in page_words(p)); grid = collections.defaultdict(collections.Counter)
    for w, c in cnt.items():
        for e in ENDINGS:
            if w.endswith(e) and len(w) > len(e): grid[w[:-len(e)]][e] += c; break
    stems = {s: g for s, g in grid.items() if sum(g.values()) >= min_stem_tokens}
    if not stems: return dict(stems=0), {}
    tot = collections.Counter(); [tot.update(g) for g in stems.values()]
    N = sum(tot.values()); He = -sum(v / N * math.log2(v / N) for v in tot.values())
    mi = 0.0
    for s, g in stems.items():
        ns = sum(g.values())
        for e, v in g.items(): mi += v / N * math.log2((v / N) / ((ns / N) * (tot[e] / N)))
    # classes: cluster stems by their ending profile; how many distinct profiles are there?
    E = sorted(tot); M = np.array([[stems[s][e] / sum(stems[s].values()) for e in E] for s in stems])
    km = KMeans(n_clusters=min(6, len(M)), random_state=0, n_init=10).fit(M); inertia_ratio = km.inertia_ / (((M - M.mean(0)) ** 2).sum() + 1e-12)
    met = dict(stems=len(stems), ending_entropy=round(He, 3), stem_ending_mi=round(mi, 3), dependence=round(mi / He, 3), class_fit_residual=round(float(inertia_ratio), 3))
    rows = []
    for c in range(km.n_clusters):
        idx = np.where(km.labels_ == c)[0]; prof = M[idx].mean(0)
        rows.append(dict(cls=c, stems=len(idx), examples=" ".join(list(stems)[i] for i in idx[:10]), profile=" ".join(f"{E[j]}:{prof[j]:.2f}" for j in np.argsort(-prof)[:5])))
    return met, dict(classes=pd.DataFrame(rows))


# ---------------------------------------------------------------- 14. A <-> B alignment of distributional spaces
def ab_alignment(pages, min_count=6, dim=30, seed=0):
    def space(ps):
        seqs = [l for p in ps for l in p["lines"]]; cnt = collections.Counter(w for s in seqs for w in s)
        vocab = [w for w, c in cnt.most_common() if c >= min_count]; vi = {w: i for i, w in enumerate(vocab)}; V = len(vocab)
        C = np.zeros((V, 2 * V))
        for s in seqs:
            for i, w in enumerate(s):
                if w not in vi: continue
                if i > 0 and s[i - 1] in vi: C[vi[w], vi[s[i - 1]]] += 1
                if i + 1 < len(s) and s[i + 1] in vi: C[vi[w], V + vi[s[i + 1]]] += 1
        tot = C.sum() + 1e-9; pr = C.sum(1, keepdims=True) / tot; pc = C.sum(0, keepdims=True) / tot
        with np.errstate(divide="ignore", invalid="ignore"): pm = np.nan_to_num(np.maximum(np.log((C / tot) / (pr * pc)), 0))
        Em = TruncatedSVD(n_components=dim, random_state=seed).fit_transform(pm); Em /= np.linalg.norm(Em, axis=1, keepdims=True) + 1e-12
        return vocab, Em
    A = [p for p in pages if p["lang"] == "A"]; B = [p for p in pages if p["lang"] == "B"]
    va, Ea = space(A); vb, Eb = space(B); ia = {w: i for i, w in enumerate(va)}; ib = {w: i for i, w in enumerate(vb)}
    shared = [w for w in va if w in ib]
    if len(shared) < 20: return dict(shared=len(shared)), {}
    rng = random.Random(seed); rng.shuffle(shared); train, test = shared[: int(len(shared) * 0.7)], shared[int(len(shared) * 0.7):]
    Xa = Ea[[ia[w] for w in train]]; Xb = Eb[[ib[w] for w in train]]; U, _, Vt = np.linalg.svd(Xa.T @ Xb); R = U @ Vt
    P = Ea @ R; S = P @ Eb.T
    hit = np.mean([vb[int(S[ia[w]].argmax())] == w for w in test])  # held-out shared words map to themselves?
    def ed(a, b):
        a, b = tokenize_glyphs(a), tokenize_glyphs(b); d = list(range(len(b) + 1))
        for i in range(1, len(a) + 1):
            prev, d[0] = d[0], i
            for j in range(1, len(b) + 1): prev, d[j] = d[j], min(d[j] + 1, d[j - 1] + 1, prev + (a[i - 1] != b[j - 1]))
        return d[len(b)]
    only = [w for w in va if w not in ib]; pairs = [(w, vb[int(S[ia[w]].argmax())]) for w in only]
    near = np.mean([ed(a, b) <= 1 for a, b in pairs]) if pairs else 0
    rnd = np.mean([ed(a, rng.choice(vb)) <= 1 for a, b in pairs]) if pairs else 0
    met = dict(shared=len(shared), heldout_self_alignment=round(float(hit), 3), a_only_words=len(only), a_only_mapped_to_edit1=round(float(near), 3), random_edit1=round(float(rnd), 3))
    return met, dict(pairs=pd.DataFrame([dict(a_word=a, b_nearest=b) for a, b in pairs[:60]]))


METHODS.update(stylometry=stylometry, anomalies=anomalies, paradigms=paradigms, ab_alignment=ab_alignment)