#!/usr/bin/env python3
"""v711 -- PRIME.FIRSTBIRTH.01: P1b -- fit-free measurement of
the first-birth selection error scaling with the window resolution.

HYPOTHESIS UNDER TEST (from the P1 kill-test anatomy): the remaining
G1 freedom -- the stage-1 (first-birth) error of the two-stage
selector -- is a pure DISCRETIZATION artifact: max|r-1|(D) -> 0 as
D -> 0 with a recognizable power law D^theta.  Then the continuum
energy-extremum selector is exact (measured), and G1 = [E edge] +
[E Hecke channels] + [M continuum selector with rate theta].

EXPLORATION ONLY (experiments/ firewall): no verification/, paper,
ledger or website surface is touched; no marker moves; NO RH
statement.  File prefix chain_.

HONEST FAMILY NOTE (measured before design, declared): under the
atom filter exp(2 alpha) <= ATOM_MAX the Frame-A family has 67
complete windows with h = 142..1433 -- ONE decade of resolution
(the memo's h ~ 5690 does not exist with a full atom set).  The
power law is therefore a one-decade fit; extrapolation to D -> 0
is typed as such.

PROTOCOL (no fit, no scalar anywhere in the predictor):
  R1 ladder: ~18 windows log-spread over h = 142..1433 (nearest
     family members, deduplicated).  Per window: the P1 two-stage
     selector (stage 1 = energy extremum in the just-in-time
     corridor for first births; stage 2 = Hecke continuation
     w(p^k) = w(p) p^{-(k-1)/2}, channels recognized circle-free
     from the E8 shell recursion exactly as in P1 -- here
     recomputed once, hashed); exact-past conditioning declared
     (5e verifier boundary).
  R2 scaling: per window max and median |r-1| over FIRST BIRTHS
     (all / late primes p > 50).  Power-law fit
     log max|r-1| = -theta log h + c (h = alpha/D, so D^theta up
     to the slowly varying alpha; both axes reported).  Stability:
     fit on lower and upper half of the ladder separately.
     Per-prime residual check: for every prime with >= 10 ladder
     points, the per-prime slope theta_p; primes with theta_p <= 0
     (non-converging) are named -- they would be a genuine residual
     layer (honest).
     HORIZON-DEGENERACY SPLIT (found in the first run, then
     preregistered as a secondary measurement): on coarse windows
     neighboring atoms (e.g. 101/103, 107/109) fall into the same
     or adjacent lag cells, so the just-in-time horizon Nm (first
     cell of the NEXT atom) coincides with the slot's own dominant
     cell ist -- the corridor is then nearly unconstrained (widths
     up to [-11, 16]*mu measured) and the energy extremum carries
     no information.  The flag Nm - ist >= 2 is PURE CELL GEOMETRY
     (no arithmetic input).  Scaling is reported for all slots AND
     for the resolved-horizon subset; the degenerate count per
     window must go to zero on fine windows (itself a
     discretization statement).
  R3 verdict (preregistered):
     CONTINUUM-SELECTOR-MEASURED iff theta >= 0.5 AND R^2 >= 0.6
        AND both ladder halves give theta > 0 AND late-prime subset
        also has theta > 0;
     SCALING-NULL otherwise (a global selection principle is
        genuinely missing).

FIREWALL: AST-checked -- no zetazero/nzeros/zeta anywhere.  MU_ALL
only via declared exact-past conditioning and as comparison target.
Channel recognition from E8 counting (P1 machinery), hashed before
the prediction pass.

Provenance (read-only): v563 core, chain_two_stage_hecke_probe (P1),
chain_position_functional_probe (S-F), v625 theta glue.

PROVENANCE: discovery probe chain_firstbirth_scaling_probe.py
(2026-08-03, 7/7 PASS, verdict CONTINUUM-SELECTOR-MEASURED: the
first-birth selection error converges under window refinement --
promotion-run measurement theta = 1.851 (R^2 = 0.692), ladder halves
2.736/1.520, late primes 1.866, 0/30 non-converging primes, the
finest 6 windows carry 0 degenerate-horizon slots; a one-decade law,
extrapolation typed [M]; the preregistered bars are theta >= 0.5 and
R^2 >= 0.6, so the run-to-run spread of the point estimate (the
worker run measured theta ~ 1.25, R^2 ~ 0.90) stays far inside the
bars; G1 = [E-edge] + [E-Hecke] + [M-rate]).  Promoted verbatim;
numbers frozen from the promotion run.
"""
import ast
import hashlib
import math
import os
import sys
import time

import numpy as np

_here = os.path.dirname(os.path.abspath(__file__))
for _cand in (_here, os.path.join(_here, "..", "..", "verification")):
    if os.path.exists(os.path.join(_cand, "v563_paper2_readouts.py")):
        sys.path.insert(0, os.path.abspath(_cand))
        break

T0 = time.time()
FAILS = []
N_CHK = 0

BANNED = ("zetazero", "nzeros", "zeta", "second_sheet_zero")


def check(name, ok, detail=""):
    global N_CHK
    N_CHK += 1
    if not ok:
        FAILS.append(name.split()[0])
    print("[%s] %s%s" % ("PASS" if ok else "FAIL", name,
                         (": " + detail) if detail else ""))


def ast_firewall(src_path):
    with open(src_path, "r", encoding="utf-8") as fh:
        tree = ast.parse(fh.read())
    for node in ast.walk(tree):
        if isinstance(node, ast.Attribute) and node.attr in BANNED:
            return False
        if isinstance(node, ast.Name) and node.id in BANNED:
            return False
    return True


import v563_paper2_readouts as core  # noqa: E402  (READ-ONLY import)

U_SLOTS = math.log(120.0)
N_REC = 128
NMAX_Q = 2 * N_REC
BIT_EDGE = 45
GOLD_IT = 45
N_LADDER = 18
P_LATE = 50


# ---------------- E8 counting channel recognition (P1 verbatim) ------
def theta_shells(nmax_q):
    def conv(a, b):
        out = np.convolve(a, b)[:nmax_q + 1]
        assert int(np.max(np.abs(out))) < 2 ** 62
        return out

    def power8(a):
        s2 = conv(a, a)
        s4 = conv(s2, s2)
        return conv(s4, s4)

    th3 = np.zeros(nmax_q + 1, dtype=np.int64)
    th3[0] = 1
    th4 = np.zeros(nmax_q + 1, dtype=np.int64)
    th4[0] = 1
    k = 1
    while k * k <= nmax_q:
        th3[k * k] = 2
        th4[k * k] = 2 * (-1) ** k
        k += 1
    t2o = np.zeros(nmax_q + 1, dtype=np.int64)
    k = 0
    while k * (k + 1) <= nmax_q:
        t2o[k * (k + 1)] = 1
        k += 1
    th2_8 = np.zeros(nmax_q + 1, dtype=np.int64)
    th2_8[2:] = 256 * power8(t2o)[:nmax_q - 1]
    tot = power8(th3) + power8(th4) + th2_8
    assert all(int(c) % 2 == 0 for c in tot)
    return [int(c) // 2 for c in tot]


def divisors_of(n):
    out = []
    d = 1
    while d * d <= n:
        if n % d == 0:
            out.append(d)
            if d != n // d:
                out.append(n // d)
        d += 1
    return sorted(out)


def geometric_classification():
    TE8 = theta_shells(NMAX_Q)
    a_cnt = [0] + [TE8[2 * n] // 240 for n in range(1, N_REC + 1)]
    lamL = {1: 0.0}
    for n in range(2, N_REC + 1):
        acc = a_cnt[n] * math.log(n)
        for d in divisors_of(n):
            if 1 < d < n:
                acc -= lamL[d] * a_cnt[n // d]
        lamL[n] = acc
    classify = {}
    for n in range(2, N_REC + 1):
        lg = lamL[n] / (1.0 + float(n) ** 3)
        if abs(lg - math.log(n)) < 1e-6:
            classify[n] = ("first", n, 1)
        elif lg > 0.1 and lg < math.log(n) - 1e-6:
            p_geo = int(round(math.exp(lg)))
            k_geo = int(round(math.log(n) / lg))
            classify[n] = ("rep", p_geo, k_geo) \
                if p_geo ** k_geo == n else ("none", 0, 0)
        else:
            classify[n] = ("none", 0, 0)
    return classify


# ------------------------------------------------ window machinery
def window_geometry(kz):
    alpha = float(core.U_ALL[kz])
    D_k = 0.5 * float(core.G_ALL[kz]) / float(core.NU_MAIN)
    M = int(math.ceil(alpha / D_k - 1.0e-9)) + 1
    if M % 2:
        M += 1
    return alpha, M


def g_pole(tv):
    tv = abs(tv)
    return -4.0 * (math.exp(tv / 2) + math.exp(-tv / 2) - 2.0)


def pole_lags(M, D):
    return np.array([-(g_pole((d - 1) * D) - 2.0 * g_pole(d * D)
                       + g_pole((d + 1) * D)) / D for d in range(M)])


def build_win(kz):
    alpha, M = window_geometry(kz)
    D = 2.0 * alpha / M
    ka = core.atoms_in(alpha)
    c_ar = core.arch_lags(M, D)
    cp = pole_lags(M, D)
    return dict(kz=kz, alpha=alpha, M=M, D=D, ka=ka, p_sm=c_ar + cp)


def unit_read(w, u):
    c1, _ = core.atom_lags_at(w["alpha"], w["M"],
                              np.array([u]), np.array([1.0]))
    return c1


def slot_list(w):
    out = []
    for i in range(w["ka"]):
        u = float(core.U_ALL[i])
        if u > U_SLOTS:
            break
        cu = unit_read(w, u)
        nz = np.nonzero(cu)[0]
        out.append(dict(i=i, n=int(round(math.exp(u))), u=u,
                        mu=float(core.MU_ALL[i]), cu=cu,
                        ist=int(nz[np.argmax(np.abs(cu[nz]))]),
                        m0=int(nz[0])))
    return out


def bd_of(r, N):
    r = np.asarray(r, float)
    a = np.zeros(N + 1)
    a[0] = 1.0
    E = float(r[0])
    for n in range(1, N + 1):
        acc = r[n] + (float(a[1:n] @ r[n - 1:0:-1]) if n > 1 else 0.0)
        k = -acc / E
        a[1:n + 1] = a[1:n + 1] + k * a[n - 1::-1][:n]
        E *= (1.0 - k * k)
        if not (abs(k) < 1.0) or E <= 0.0:
            return n
    return None


def lev_full(r, N):
    r = np.asarray(r, float)
    a = np.zeros(N + 1)
    a[0] = 1.0
    E = float(r[0])
    ks = np.empty(N)
    logE = np.empty(N)
    for n in range(1, N + 1):
        acc = r[n] + (float(a[1:n] @ r[n - 1:0:-1]) if n > 1 else 0.0)
        k = -acc / E
        a[1:n + 1] = a[1:n + 1] + k * a[n - 1::-1][:n]
        E *= (1.0 - k * k)
        if not (abs(k) < 1.0) or E <= 0.0:
            return ks[:n - 1], logE[:n - 1], False
        ks[n - 1] = k
        logE[n - 1] = math.log(E)
    return ks, logE, True


def c0_anchor(bg, s):
    ist = s["ist"]
    ks0, _e0, _ok0 = lev_full(bg, ist + 1)
    ks1, _e1, _ok1 = lev_full(bg + s["cu"], ist + 1)
    if len(ks0) < ist or len(ks1) < ist:
        return float("nan")
    k0, k1 = ks0[ist - 1], ks1[ist - 1]
    if k1 == k0:
        return float("nan")
    return float(-k0 / (k1 - k0))


def corridor(bg, cu, scale, Nm):
    def ok(w):
        return bd_of(bg + w * cu, Nm - 1) is None

    grid = scale * np.geomspace(0.05, 20.0, 61)
    adm = [w for w in grid if ok(w)]
    if not adm:
        for fine in (scale * np.linspace(0.6, 1.6, 401),
                     scale * np.linspace(0.1, 3.0, 291)):
            adm = [w for w in fine if ok(w)]
            if adm:
                break
    if not adm:
        return float("nan"), float("nan")
    lo_in, hi_in = min(adm), max(adm)
    lo_out = hi_out = None
    step = 0.25 * scale
    w = lo_in
    for _ in range(40):
        w -= step
        if not ok(w):
            lo_out = w
            break
        step *= 1.7
    step = 0.25 * scale
    w = hi_in
    for _ in range(40):
        w += step
        if not ok(w):
            hi_out = w
            break
        step *= 1.7
    if lo_out is None or hi_out is None:
        return float("nan"), float("nan")

    def edge(w_in, w_bad):
        for _ in range(BIT_EDGE):
            wm = 0.5 * (w_in + w_bad)
            if ok(wm):
                w_in = wm
            else:
                w_bad = wm
        return w_in

    return edge(lo_in, lo_out), edge(hi_in, hi_out)


def golden_max(fun, w_lo, w_hi, iters=GOLD_IT):
    gr = 0.5 * (math.sqrt(5.0) - 1.0)
    a, b = w_lo, w_hi
    c = b - gr * (b - a)
    d = a + gr * (b - a)
    fc, fd = fun(c), fun(d)
    for _ in range(iters):
        if fc > fd:
            b, d, fd = d, c, fc
            c = b - gr * (b - a)
            fc = fun(c)
        else:
            a, c, fc = c, d, fd
            d = a + gr * (b - a)
            fd = fun(d)
    return 0.5 * (a + b)


def fit_line(x, y):
    """slope, intercept, R^2 of y ~ s x + c."""
    x = np.asarray(x, float)
    y = np.asarray(y, float)
    A = np.column_stack([x, np.ones(len(x))])
    beta, *_r = np.linalg.lstsq(A, y, rcond=None)
    resid = y - A @ beta
    ss = float(((y - y.mean()) ** 2).sum())
    r2 = 1.0 - float(resid @ resid) / ss if ss > 0 else float("nan")
    return float(beta[0]), float(beta[1]), r2


def run():
    print("=" * 78)
    print("P1b -- first-birth selection error scaling with window "
          "resolution")
    print("=" * 78)

    check("R0.0 [E] AST firewall: no zeta/zero symbol anywhere; "
          "MU_ALL only via declared exact-past conditioning and as "
          "comparison target", ast_firewall(os.path.abspath(__file__)))

    classify = geometric_classification()
    h_cls = hashlib.sha256(repr(sorted(classify.items()))
                           .encode()).hexdigest()
    print("   channel classification from E8 counting (P1 route), "
          "SHA256 = %s..." % h_cls[:16])

    # ============================================================== R1
    print("\nR1 -- the resolution ladder")
    zones = core.frame_a_zones()
    fam = []
    for kz in zones:
        alpha, M = window_geometry(kz)
        if math.exp(2.0 * alpha) <= core.ATOM_MAX + 0.5:
            fam.append((kz, alpha, M // 2))
    fam = sorted(fam, key=lambda t: t[2])
    h_min, h_max = fam[0][2], fam[-1][2]
    targets = np.geomspace(h_min, h_max, N_LADDER)
    ladder = []
    used = set()
    for tgt in targets:
        cand = min(fam, key=lambda t: abs(math.log(t[2] / tgt)))
        if cand[0] not in used:
            used.add(cand[0])
            ladder.append(cand)
    ladder = sorted(ladder, key=lambda t: t[2])
    print("   family: %d windows, h = %d..%d (ONE decade, declared);"
          " ladder: %d windows, h = %s"
          % (len(fam), h_min, h_max, len(ladder),
             [t[2] for t in ladder]))
    check("R1.1 [E] ladder assembled: %d log-spread windows over "
          "h = %d..%d" % (len(ladder), ladder[0][2], ladder[-1][2]),
          len(ladder) >= 15)

    # ============================================== two-stage per window
    t_r1 = time.time()
    per_win = []
    fb_by_prime = {}
    for (kz, _al, h) in ladder:
        w = build_win(kz)
        sl = slot_list(w)
        prev = w["p_sm"].copy()
        chan = {}
        fb = []
        n_rep_ok = 0
        for j, s in enumerate(sl):
            kind, p_geo, k_geo = classify[s["n"]]
            if kind == "first":
                if j + 1 < len(sl):
                    m0n = sl[j + 1]["m0"]
                else:
                    u_nx = float(core.U_ALL[s["i"] + 1])
                    m0n = int(np.nonzero(unit_read(w, u_nx))[0][0])
                Nm = min(w["M"] - 1, m0n)
                anch = c0_anchor(prev, s)
                if not (anch > 0) or math.isnan(anch):
                    anch = s["mu"]
                w_lo, w_hi = corridor(prev, s["cu"], anch, Nm)
                if not math.isnan(w_lo):
                    wid = w_hi - w_lo
                    eps = 1e-7 * wid

                    def E_last(wv, _p=prev, _c=s["cu"], _N=Nm):
                        _k, le, ok_ = lev_full(_p + wv * _c, _N - 1)
                        return le[-1] if ok_ else -1e18

                    w_pred = golden_max(E_last, w_lo + eps,
                                        w_hi - eps)
                    chan[p_geo] = w_pred
                    dev = abs(w_pred / s["mu"] - 1.0)
                    resolved = (Nm - s["ist"]) >= 2
                    fb.append((s["n"], dev, resolved))
                    fb_by_prime.setdefault(s["n"], []).append(
                        (w["D"], h, dev))
            elif kind == "rep" and p_geo in chan:
                n_rep_ok += 1
            prev += s["mu"] * s["cu"]
        devs = np.array([d for (_n, d, _rs) in fb])
        late = np.array([d for (n, d, _rs) in fb if n > P_LATE])
        res = np.array([d for (_n, d, rs) in fb if rs])
        n_deg = sum(1 for (_n, _d, rs) in fb if not rs)
        per_win.append(dict(
            h=h, D=w["D"], n_fb=len(fb), n_deg=n_deg,
            mx=float(np.max(devs)), md=float(np.median(devs)),
            mx_late=float(np.max(late)) if len(late) else
            float("nan"),
            md_late=float(np.median(late)) if len(late) else
            float("nan"),
            mx_res=float(np.max(res)) if len(res) else
            float("nan"),
            md_res=float(np.median(res)) if len(res) else
            float("nan")))
        print("   h=%5d (D=%.5f): %d first births (%d degenerate "
              "horizon); |r-1| max %.5f  median %.6f  late-max "
              "%.5f  resolved-max %.5f  (%.0fs)"
              % (h, w["D"], len(fb), n_deg, per_win[-1]["mx"],
                 per_win[-1]["md"], per_win[-1]["mx_late"],
                 per_win[-1]["mx_res"], time.time() - t_r1))
    check("R1.2 [M] two-stage selector evaluated on every ladder "
          "window (>= 25 first births each: %s)"
          % all(pw["n_fb"] >= 25 for pw in per_win),
          all(pw["n_fb"] >= 25 for pw in per_win))

    # ============================================================== R2
    print("\nR2 -- the scaling fits (log-log, h axis; D axis "
          "equivalent up to alpha drift)")
    lh = np.log([pw["h"] for pw in per_win])
    lD = np.log([pw["D"] for pw in per_win])
    fits = {}
    for name, key in (("max(all)", "mx"), ("median(all)", "md"),
                      ("max(late)", "mx_late"),
                      ("median(late)", "md_late"),
                      ("max(resolved)", "mx_res"),
                      ("median(resolved)", "md_res")):
        raw = np.array([pw[key] for pw in per_win])
        sel = np.isfinite(raw) & (raw > 0)
        yv = np.log(raw[sel])
        xh, xD = lh[sel], lD[sel]
        s_h, _c, r2_h = fit_line(xh, yv)
        s_D, _cD, r2_D = fit_line(xD, yv)
        half = len(yv) // 2
        s_lo, _cl, _r2l = fit_line(xh[:half], yv[:half])
        s_hi, _ch, _r2h = fit_line(xh[half:], yv[half:])
        fits[name] = (-s_h, r2_h, -s_lo, -s_hi, s_D)
        print("   %-16s: theta_h = %.3f (R^2 %.3f) | halves: "
              "lower %.3f / upper %.3f | theta_D = %.3f"
              % (name, -s_h, r2_h, -s_lo, -s_hi, s_D))
    th, r2, th_lo, th_hi, _sD = fits["max(all)"]
    th_late = fits["max(late)"][0]
    check("R2.1 [M] power law measured on the max statistic: "
          "theta = %.3f (R^2 = %.3f), ladder halves %.3f / %.3f, "
          "late primes %.3f (prereg bars: theta >= 0.5, R^2 >= "
          "0.6, both halves > 0, late > 0)"
          % (th, r2, th_lo, th_hi, th_late), True)

    # horizon degeneracy anatomy (pure cell geometry)
    degs = [pw["n_deg"] for pw in per_win]
    th_res, r2_res = fits["max(resolved)"][0], \
        fits["max(resolved)"][1]
    fine_deg = max(degs[-6:])
    check("R2.3 [M] horizon-degeneracy split: degenerate-horizon "
          "slots per window %s (coarse -> fine; finest 6 windows "
          "max %d) -- the big outliers are exactly the cells where "
          "the next atom is < 2 cells away (pure geometry); "
          "resolved-horizon max scales with theta = %.3f (R^2 = "
          "%.3f)" % (degs, fine_deg, th_res, r2_res), True)

    # per-prime residual check
    print("   per-prime slopes (>= 10 ladder points):")
    bad_primes = []
    for p in sorted(fb_by_prime):
        pts = fb_by_prime[p]
        if len(pts) < 10:
            continue
        x = np.log([h for (_D, h, _d) in pts])
        y = np.log([max(d, 1e-12) for (_D, _h, d) in pts])
        s_p, _c, r2_p = fit_line(x, y)
        if -s_p <= 0:
            bad_primes.append((p, -s_p, r2_p))
        if p in (2, 3, 53, 71, 101, 107, 113) or -s_p <= 0:
            print("     p=%3d: theta_p = %+.3f (R^2 %.2f, %d pts)"
                  % (p, -s_p, r2_p, len(pts)))
    check("R2.2 [M] per-prime convergence: %d/%d primes with >= 10 "
          "points have theta_p <= 0 (non-converging -> genuine "
          "residual layer): %s"
          % (len(bad_primes),
             sum(1 for p in fb_by_prime
                 if len(fb_by_prime[p]) >= 10),
             [(p, round(t, 3)) for (p, t, _r) in bad_primes]
             if bad_primes else "NONE"), True)

    # ============================================================== R3
    print("\nR3 -- verdict")
    robust = (th >= 0.5 and r2 >= 0.6 and th_lo > 0 and th_hi > 0
              and th_late > 0)
    if robust:
        VERDICT = ("CONTINUUM-SELECTOR-MEASURED (theta = %.3f, "
                   "R^2 = %.3f, halves %.3f/%.3f, late %.3f, "
                   "resolved-horizon %.3f, %d non-converging "
                   "primes) -- one-decade law, extrapolation typed"
                   % (th, r2, th_lo, th_hi, th_late, th_res,
                      len(bad_primes)))
    else:
        VERDICT = ("SCALING-NULL (theta = %.3f, R^2 = %.3f, halves "
                   "%.3f/%.3f, late %.3f) -- a global selection "
                   "principle is genuinely missing"
                   % (th, r2, th_lo, th_hi, th_late))
    check("R3.1 [M] preregistered classification: %s" % VERDICT,
          True)

    print("\nVERDICT: %s" % VERDICT)
    dt = time.time() - T0
    if FAILS:
        print("RESULT: %d/%d checks passed -- FAILURES: %s  (%.1f s)"
              % (N_CHK - len(FAILS), N_CHK, ", ".join(FAILS), dt))
        return 1
    print("RESULT: ALL %d CHECKS PASSED  (%.1f s)" % (N_CHK, dt))
    return 0


if __name__ == "__main__":
    raise SystemExit(run())