#!/usr/bin/env python3
"""voynich-evaluate.py — recompute headline numbers of the paper
"How the Voynich Manuscript was written" (voynich.pdf on this site) from the
Zandbergen–Landini transliteration, and check them against the values
recorded by the paper's archived analysis (the Figure 3 legend; the
per-page results behind Section 3.1). Nothing is echoed from stored
results: every number below is rebuilt from the raw transliteration bytes
at each run.

Input (downloaded separately, verified fail-closed)
---------------------------------------------------
ZL3b-n.txt — the Zandbergen–Landini (ZL) transliteration of the Voynich
manuscript, IVTFF Eva- 2.0 format, version 3b of 13/05/2025, 411,671 bytes.
The file is not distributed with this bundle: download it from voynich.nu
and place it beside this script, or point --data at a copy. The paper's
analysis reads exactly this file, and the script checks the copy it is
given against the recorded fingerprint
    sha256 bf5b6d4ac1e3a51b1847a9c388318d609020441ccd56984c901c32b09beccafc
before anything is computed, refusing any other bytes (exit 2) — so a
passing run always refers to the same public input of record.

What is recomputed and checked
------------------------------
1. Corpus slice: 227 text pages, 38,718 kept word tokens, 8,089 word types,
   under the paper's tokenization rules (below).
2. Dialect anchors: 112 pages / 10,749 tokens (dialect A: Scribe 1 pages
   Currier-labelled A) and 43 pages / 10,508 tokens (dialect B: Scribe 2
   pages Currier-labelled B), astronomical/cosmological/zodiac pages set
   aside, per the paper's role assignment.
3. The per-page composition map over all 227 pages: 116 pure dialect A,
   66 pure dialect B, 43 blends (25 decisive + 18 leaning), 2 ambiguous —
   the counts of the paper's Figure 3 legend.
4. Folio 26v, the page the paper walks through the classifier:
   log10 BF(A:B) = -46.378, decisively pure B — the value in the archived
   per-page results behind Section 3.1 (the paper's text prints it as
   "more than 10^20").
5. Numerics cross-check on folio 26v: the fast log-space evidence used for
   the sweep is recomputed with exact big-integer arithmetic and must agree
   to 1e-9 in log10.
6. The label register: 864 primary label tokens, observed o-initial rate
   0.59375, against the exact page-matched null expectation 0.31977
   (the paper's 59.4% vs 32.0%).

Tokenization (the paper's pin, applied to the raw IVTFF text)
-------------------------------------------------------------
Comment lines dropped; latin-1 encoding. A page is each <fN> header line,
with page variables $Q,$P,$I,$L,$H recorded ('-' when absent). All locus
lines of a page are pooled; the locus type letter (P paragraph, L label, ...)
is kept per locus. In-text cleanup, in order: inline comments <!...> removed;
paragraph marks <%> <$> removed; gaps <-> <~> become word breaks; alternate
readings [x:y] -> x (the ZL primary); ligature braces {...} keep content;
rare-glyph codes @NNN; become a placeholder; alignment filler ! removed; any
other <...> tag removed. Word breaks are '.', ',' and whitespace. Tokens
containing an unreadable glyph '?' or a rare-glyph placeholder are excluded
from all counts.

Models (as in the paper)
------------------------
Dialect profiles are frozen Dirichlet(1)-smoothed frequency vectors over the
full k = 8,089 word-type vocabulary, pooled over the anchor pages:
p_v = (c_v + 1)/(T + k). For a page with word counts u:
    pure evidence   Z_A = prod_v p_A,v^{u_v}   (same for B)
    blend evidence  Z_bl = integral_0^1 prod_v (w p_A,v + (1-w) p_B,v)^{u_v} dw
(uniform Beta(1,1) prior on the mixing weight w; the integral is a polynomial
in w and is evaluated exactly). Classification bars: a page is dialect A when
log10 BF(A:B) >= 2, dialect B when <= -2, ambiguous otherwise; it is a blend
when log10 BF(blend : best pure) > 0, decisively so at >= 2, and pure
otherwise. The sweep evaluates the blend integral by a log-space
Bernstein-lattice recursion (all coefficients nonnegative, so log-sum-exp is
stable); check 5 recomputes folio 26v with exact rational arithmetic — the
smallest margin any page has to a classification bar is 0.024 in log10,
about ten orders of magnitude above the observed agreement of the two routes.

Label register: primary labels are the tokens of single-token L loci
(864 of them). The null pairs each labelled page with its own non-label
tokens and asks how often a token drawn from that pool starts with 'o'; the
reported null is the exact expectation of that page-matched draw,
    E = (1/n) * sum_pages n_page * (o-initial fraction of the page's pool),
computed with exact rationals (drawing with or without replacement leaves
the expectation unchanged). Pages whose pool is empty fall back to the pool
of same-illustration same-hand pages (no page needs the fallback on this
input).

Modes and exit codes
--------------------
default      verify the sha256 pin, rebuild everything, print each
             recomputed value beside its pinned expectation with PASS/FAIL.
             Exit 0 only if every check passes; exit 1 on any mismatch.
--data FILE  read the transliteration from FILE instead of ZL3b-n.txt
             beside the script; the fingerprint check applies unchanged.
--mutate     control run: verify the sha256 pin, then perturb the decoded
             text in memory before anything is computed — (a) every locus
             line of folio 26v is deleted, (b) the leading 'o' of each
             label locus on page f88r becomes 'a' — and run the same
             computation. The corpus, anchor, composition, f26v and label
             checks must then genuinely fail: the expected outcome is FAIL
             lines and exit 1. If the perturbed input somehow still passes,
             the control itself is broken and the script exits 4.
-h, --help   print usage and exit 0.
exit 2       input file missing (the message says where to download it) or
             its sha256/length does not match the pin.
exit 3       usage error: unknown flag, or --data without a file.

Python 3.8+ standard library only; no network, no other files. The default
run finishes in a few seconds on ordinary hardware.
"""
import hashlib
import os
import re
import sys
import time
from collections import Counter
from fractions import Fraction
from math import exp, inf, lgamma, log, log10

PIN_SHA = 'bf5b6d4ac1e3a51b1847a9c388318d609020441ccd56984c901c32b09beccafc'
PIN_BYTES = 411671
SRC_NAME = 'ZL3b-n.txt'

# ---- pinned expectations (the paper's recorded values) ----------------------
EXPECT = {
    'corpus': (227, 38718, 8089),          # pages, kept tokens, word types
    'anchor_A': (112, 10749),              # pages, tokens
    'anchor_B': (43, 10508),
    'composition': (116, 66, 43),          # pure A, pure B, blends
    'blend_split': (25, 18),               # decisive, leaning
    'ambiguous': 2,
    'f26v_BF': -46.378,                    # log10 BF(A:B), 3 decimals
    'labels_n': 864,
    'labels_obs': 0.59375,                 # o-initial rate, 4 decimals
    'labels_null': 0.31977,                # exact matched expectation, 4 dec.
}
TOL_BF = 5e-4      # 3 decimals
TOL_RATE = 5e-5    # 4 decimals
TOL_XCHK = 1e-9    # float vs exact arithmetic, log10

USAGE = """\
usage: voynich-evaluate.py [--data FILE] [--mutate] [-h | --help]

Recompute the composition-map and label-register numbers of "How the
Voynich Manuscript was written" from the ZL transliteration and check them
against the values recorded by the paper's archived analysis.

  --data FILE  transliteration file to read (default: ZL3b-n.txt beside
               this script); its sha256 and length must match the pin
  --mutate     control run on a perturbed in-memory copy of the input; the
               checks must then fail and the exit code is nonzero
  -h, --help   print this usage and exit 0

The default run finishes in a few seconds on ordinary hardware.
Exit codes: 0 all checks pass; 1 a check failed; 2 input missing or not
matching the pinned fingerprint; 3 usage error; 4 broken mutate control.
"""


def parse_args(argv):
    """Strict flag parsing: unknown arguments are refused (exit 3)."""
    mutate = False
    data = None
    i = 0
    while i < len(argv):
        a = argv[i]
        if a in ('-h', '--help'):
            print(USAGE, end='')
            sys.exit(0)
        elif a == '--mutate':
            mutate = True
        elif a == '--data':
            i += 1
            if i == len(argv):
                sys.stderr.write('error: --data needs a file argument\n\n'
                                 + USAGE)
                sys.exit(3)
            data = argv[i]
        elif a.startswith('--data='):
            data = a[len('--data='):]
        else:
            sys.stderr.write(f'error: unknown argument {a!r}\n\n' + USAGE)
            sys.exit(3)
        i += 1
    return mutate, data

# ---- IVTFF parsing (the paper's tokenization pin) ---------------------------
HDR = re.compile(r'^<(f[^.>]+)>\s*<!(.*)>\s*$')
LOC = re.compile(r'^<(f[^.>]+)\.([0-9a-zA-Z]+),([@+*=&~/])([A-Za-z])'
                 r'([^>]*)>\s?(.*)$')
VAR = re.compile(r'\$(\w+)=(\S+)')


def tokens_of(text):
    """Cleanup + word-break rules for one locus; returns kept tokens."""
    t = re.sub(r'<![^>]*>', '', text)              # inline comments
    t = t.replace('<%>', '').replace('<$>', '')    # paragraph marks
    t = t.replace('<->', '.').replace('<~>', '.')  # gaps -> word break
    t = re.sub(r'<[^>]*>', '', t)                  # any other tag
    t = re.sub(r'\[([^:\]]*):[^\]]*\]', r'\1', t)  # alternates -> primary
    t = t.replace('{', '').replace('}', '')        # ligature braces
    t = re.sub(r'@[0-9]+;', '\x01', t)             # rare-glyph placeholder
    t = t.replace('!', '')                         # alignment filler
    raw = [w for w in re.split(r'[., \t]+', t) if w]
    return [w for w in raw if '?' not in w and '\x01' not in w]


def parse(text):
    """Returns (page_order, pages); pages[pg] = {'vars', 'loci'} with
    loci entries (type_letter, kept_tokens)."""
    order, pages = [], {}
    cur = None
    for line in text.split('\n'):
        if line.startswith('#'):
            continue
        m = HDR.match(line)
        if m:
            cur = m.group(1)
            if cur in pages:
                raise ValueError('duplicate page ' + cur)
            v = dict(VAR.findall(m.group(2)))
            pages[cur] = {'vars': {k: v.get(k, '-')
                                   for k in ('Q', 'P', 'I', 'L', 'H')},
                          'loci': []}
            order.append(cur)
            continue
        m = LOC.match(line)
        if m:
            if m.group(1) != cur:
                raise ValueError('locus outside its page: ' + repr(line))
            pages[cur]['loci'].append((m.group(4), tokens_of(m.group(6))))
        elif line.strip():
            raise ValueError('unparsed line: ' + repr(line))
    return order, pages


# ---- evidence engine --------------------------------------------------------
def blend_log10(u_items, la, lb):
    """log10 of the Beta(1,1) blend evidence, log-space lattice.

    u_items: (vocab index, count) pairs; la/lb: natural-log profiles.
    lc[j] tracks the log coefficient of w^j (1-w)^(n-j) in the product
    polynomial prod_v (w p_A,v + (1-w) p_B,v)^{u_v}; all coefficients are
    nonnegative so log-sum-exp accumulation is stable."""
    lc = [0.0]
    n = 0
    for v, uv in u_items:
        av, bv = la[v], lb[v]
        if uv == 1:
            ch = (bv, av)
        else:
            ch = tuple(lgamma(uv + 1) - lgamma(i + 1) - lgamma(uv - i + 1)
                       + i * av + (uv - i) * bv for i in range(uv + 1))
        new = [-inf] * (n + uv + 1)
        for j, s in enumerate(lc):
            if s == -inf:
                continue
            for i, t in enumerate(ch):
                x, y = new[j + i], s + t
                if x < y:
                    x, y = y, x
                if y > -inf:
                    d = y - x
                    x = x + log(1.0 + exp(d)) if d > -45.0 else x
                new[j + i] = x
        lc = new
        n += uv
    # integrate against Beta(1,1): sum_j c_j * j! (n-j)! / (n+1)!
    terms = [lc[j] + lgamma(1 + j) + lgamma(1 + n - j) for j in range(n + 1)]
    mx = max(terms)
    tot = mx + log(sum(exp(t - mx) for t in terms))
    return (tot - lgamma(2 + n)) / log(10.0)


def blend_log10_exact(u_items, numsA, numsB, D):
    """The same blend evidence with exact big-integer arithmetic (both
    profiles share the denominator D = T + k of their own pool; numsA/numsB
    are the integer numerators c_v + 1). Used as the cross-check route."""
    S = [1]                       # S[j] = c_j * DA^j * DB^(n-j), integers
    DA, DB = D
    n = 0
    for v, uv in u_items:
        a, b = numsA[v], numsB[v]
        ch = [0] * (uv + 1)
        binom = 1
        for i in range(uv + 1):
            ch[i] = binom * a ** i * b ** (uv - i)
            binom = binom * (uv - i) // (i + 1)
        new = [0] * (n + uv + 1)
        for j, s in enumerate(S):
            if s:
                for i, t in enumerate(ch):
                    new[j + i] += s * t
        S = new
        n += uv
    fact = [1] * (n + 2)
    for j in range(1, n + 2):
        fact[j] = fact[j - 1] * j
    powA = [1] * (n + 1)
    powB = [1] * (n + 1)
    for j in range(1, n + 1):
        powA[j] = powA[j - 1] * DA
        powB[j] = powB[j - 1] * DB
    tot = 0
    for j in range(n + 1):
        tot += S[j] * powA[n - j] * powB[j] * fact[j] * fact[n - j]
    den = powA[n] * powB[n] * fact[n + 1]
    z = Fraction(tot, den)
    return log10(z.numerator) - log10(z.denominator)


# ---- checks -----------------------------------------------------------------
class Report:
    def __init__(self):
        self.fails = 0

    def check(self, name, ok, detail):
        tag = 'PASS' if ok else 'FAIL'
        if not ok:
            self.fails += 1
        print(f'  [{tag}] {name}: {detail}', flush=True)


def main(mutate, data):
    t0 = time.time()
    print('voynich-evaluate: recomputing the composition map and label '
          'register from the pinned ZL transliteration'
          + (' [MUTATE CONTROL: perturbed input, checks must fail]'
             if mutate else ''), flush=True)

    here = os.path.dirname(os.path.abspath(__file__))
    src = data if data else os.path.join(here, SRC_NAME)
    if not os.path.exists(src):
        print(f'REFUSED: {src} not found — download ZL3b-n.txt (IVTFF '
              'Eva- 2.0, v3b of 13/05/2025) from voynich.nu and place it '
              'beside this script, or pass --data FILE', flush=True)
        return 2
    blob = open(src, 'rb').read()
    sha = hashlib.sha256(blob).hexdigest()
    if sha != PIN_SHA or len(blob) != PIN_BYTES:
        print(f'REFUSED: {src} sha256 {sha} ({len(blob)} bytes) does '
              f'not match the pin {PIN_SHA[:16]}... ({PIN_BYTES} bytes)',
              flush=True)
        return 2
    print(f'  input {os.path.basename(src)}: {len(blob)} bytes, sha256 '
          f'verified ({sha[:16]}...)', flush=True)

    text = blob.decode('latin-1')
    if mutate:
        lines = text.split('\n')
        n_dropped = len([l for l in lines if l.startswith('<f26v.')])
        lines = [l for l in lines if not l.startswith('<f26v.')]
        lab = re.compile(r'^(<f88r\.[0-9a-zA-Z]+,[@+*=&~/]L[^>]*>\s*)o')
        n_flipped = 0
        for i, l in enumerate(lines):
            m = lab.match(l)
            if m:
                lines[i] = m.group(1) + 'a' + l[len(m.group(0)):]
                n_flipped += 1
        text = '\n'.join(lines)
        print(f'  mutation applied in memory: {n_dropped} f26v locus lines '
              f'deleted, {n_flipped} f88r label initials o->a', flush=True)

    rep = Report()

    # 1. corpus slice
    order, pages = parse(text)
    toks = {pg: [w for _, kept in pages[pg]['loci'] for w in kept]
            for pg in order}
    vocab = sorted({w for ws in toks.values() for w in ws})
    vidx = {w: i for i, w in enumerate(vocab)}
    K = len(vocab)
    ntok = sum(len(ws) for ws in toks.values())
    ep, et, ev = EXPECT['corpus']
    rep.check('corpus slice',
              (len(order), ntok, K) == EXPECT['corpus'],
              f'{len(order)} pages, {ntok} tokens, {K} types '
              f'(expect {ep}/{et}/{ev})')

    # 2. anchors (role pin: astronomical/cosmological/zodiac pages set aside,
    #    then Scribe 1 x Currier A and Scribe 2 x Currier B)
    def role(v):
        if v['I'] in ('A', 'C', 'Z'):
            return 'C'
        if v['H'] == '1' and v['L'] == 'A':
            return 'A'
        if v['H'] == '2' and v['L'] == 'B':
            return 'B'
        return 'T'

    roles = {pg: role(pages[pg]['vars']) for pg in order}
    Apg = [p for p in order if roles[p] == 'A']
    Bpg = [p for p in order if roles[p] == 'B']
    tA = sum(len(toks[p]) for p in Apg)
    tB = sum(len(toks[p]) for p in Bpg)
    rep.check('dialect-A anchor', (len(Apg), tA) == EXPECT['anchor_A'],
              f'{len(Apg)} pages / {tA} tokens '
              f'(expect {EXPECT["anchor_A"][0]}/{EXPECT["anchor_A"][1]})')
    rep.check('dialect-B anchor', (len(Bpg), tB) == EXPECT['anchor_B'],
              f'{len(Bpg)} pages / {tB} tokens '
              f'(expect {EXPECT["anchor_B"][0]}/{EXPECT["anchor_B"][1]})')

    # 3. frozen profiles and the per-page sweep
    def pooled(pgs):
        C = [0] * K
        for p in pgs:
            for w in toks[p]:
                C[vidx[w]] += 1
        return C

    CA, CB = pooled(Apg), pooled(Bpg)
    denA, denB = sum(CA) + K, sum(CB) + K
    la = [log(c + 1) - log(denA) for c in CA]
    lb = [log(c + 1) - log(denB) for c in CB]
    L10 = log(10.0)

    nA = nB = nBl = nAmb = nDec = nLean = 0
    f26v = None
    for pg in order:
        cnt = Counter(toks[pg])
        u_items = sorted(((vidx[w], c) for w, c in cnt.items()),
                         key=lambda x: -x[1])
        lZA = sum(c * la[v] for v, c in u_items) / L10
        lZB = sum(c * lb[v] for v, c in u_items) / L10
        lBl = blend_log10(u_items, la, lb) if u_items else 0.0
        ab = lZA - lZB
        bl = lBl - max(lZA, lZB)
        if bl > 0:
            nBl += 1
            if bl >= 2:
                nDec += 1
            else:
                nLean += 1
        elif ab >= 2:
            nA += 1
        elif ab <= -2:
            nB += 1
        else:
            nAmb += 1
        if pg == 'f26v':
            f26v = (ab, bl, u_items)
    eA, eB, eBlend = EXPECT['composition']
    rep.check('composition map',
              (nA, nB, nBl, nAmb) == (eA, eB, eBlend, EXPECT['ambiguous']),
              f'{nA} pure A / {nB} pure B / {nBl} blends / {nAmb} ambiguous '
              f'(expect {eA}/{eB}/{eBlend}/{EXPECT["ambiguous"]}; '
              f'total {nA + nB + nBl + nAmb})')
    rep.check('blend split',
              (nDec, nLean) == EXPECT['blend_split'],
              f'{nDec} decisive + {nLean} leaning '
              f'(expect {EXPECT["blend_split"][0]}+{EXPECT["blend_split"][1]})')

    # 4. folio 26v
    if f26v is None:
        rep.check('folio 26v', False, 'page absent from the parsed corpus')
        ab26 = bl26 = 0.0
        u26 = []
    else:
        ab26, bl26, u26 = f26v
        ok = abs(ab26 - EXPECT['f26v_BF']) <= TOL_BF
        rep.check('folio 26v log10 BF(A:B)', ok,
                  f'{ab26:+.6f} (expect {EXPECT["f26v_BF"]:+.3f} '
                  f'to 3 decimals); blend vs best pure {bl26:+.3f} -> '
                  + ('pure B' if ab26 <= -2 and bl26 <= 0 else 'NOT pure B'))

    # 5. exact-arithmetic cross-check on folio 26v
    numsA = None
    if u26:
        numsA = [c + 1 for c in CA]
        numsB = [c + 1 for c in CB]
        exact_pure = (sum(c * (log10(numsA[v]) - log10(denA))
                          for v, c in u26)
                      - sum(c * (log10(numsB[v]) - log10(denB))
                            for v, c in u26))
        exact_bl = blend_log10_exact(u26, numsA, numsB, (denA, denB))
        fast_bl = blend_log10(u26, la, lb)
        d1 = abs(exact_pure - ab26)
        d2 = abs(exact_bl - fast_bl)
        rep.check('f26v numerics cross-check',
                  d1 <= TOL_XCHK and d2 <= TOL_XCHK,
                  f'fast vs exact arithmetic: |dBF|={d1:.2e}, '
                  f'|dlog10 Z_blend|={d2:.2e} (bar {TOL_XCHK:.0e})')

    # 6. label register
    by_page, pool_by_page = {}, {}
    for pg in order:
        for tl, kept in pages[pg]['loci']:
            if tl == 'L':
                if len(kept) == 1:
                    by_page.setdefault(pg, []).append(kept[0])
            else:
                pool_by_page.setdefault(pg, []).extend(kept)
    pools = {}
    for pg in by_page:
        pool = pool_by_page.get(pg, [])
        if not pool:
            ih = (pages[pg]['vars']['I'], pages[pg]['vars']['H'])
            pool = [w for q, ws in pool_by_page.items() for w in ws
                    if (pages[q]['vars']['I'], pages[q]['vars']['H']) == ih]
        pools[pg] = pool
    n_lab = sum(len(v) for v in by_page.values())
    if n_lab == 0:
        rep.check('label register', False, 'no primary label tokens found')
    else:
        n_o = sum(1 for ws in by_page.values() for w in ws if w[0] == 'o')
        obs = n_o / n_lab
        exact = sum(len(by_page[pg])
                    * Fraction(sum(1 for w in pools[pg] if w[0] == 'o'),
                               len(pools[pg]))
                    for pg in by_page) / n_lab
        exact_f = float(exact)
        rep.check('label count', n_lab == EXPECT['labels_n'],
                  f'{n_lab} primary label tokens '
                  f'(expect {EXPECT["labels_n"]})')
        rep.check('label o-initial rate',
                  abs(obs - EXPECT['labels_obs']) <= TOL_RATE,
                  f'observed {obs:.5f} ({n_o}/{n_lab}; expect '
                  f'{EXPECT["labels_obs"]:.5f} to 4 decimals)')
        rep.check('label matched null',
                  abs(exact_f - EXPECT['labels_null']) <= TOL_RATE,
                  f'exact page-matched expectation {exact_f:.8f} '
                  f'(expect {EXPECT["labels_null"]:.5f} to 4 decimals)')

    wall = time.time() - t0
    verdict = 'PASS' if rep.fails == 0 else 'FAIL'
    print(f'OVERALL: {verdict} — {rep.fails} of the checks failed '
          f'({wall:.1f} s)', flush=True)
    if mutate:
        if rep.fails:
            print('mutate control: the perturbed input fails as it must; '
                  'the nonzero exit below is the expected outcome.',
                  flush=True)
            return 1
        print('MUTATE CONTROL BROKEN: the perturbed input still passes '
              'every check.', flush=True)
        return 4
    return 0 if rep.fails == 0 else 1


if __name__ == '__main__':
    sys.exit(main(*parse_args(sys.argv[1:])))