"""The scattering side: the theory's expected response, worked out from the equations (October 14, 2026).

A calculation, not a simulation: no oscillators are run and no saved data are read. It uses only the model's list of
natural rates (the same formula as nudge_test.natural_rates) and the large-group dispersion relation for the model's
two-group equations,

    1 = (K_mode / 2) F(lambda),    F(lambda) = (1/N) sum over pulled oscillators of 1 / (lambda + i omega_j),

with K_mode = K + eps = 2.0 for A and B moving together and K - eps = 1.8 for A and B moving oppositely. A root with
positive real part is a pattern of order that grows; its imaginary part is the pattern's turning rate. A large set of
quantiles (400,000) stands in for the large-group limit; the model's own 4,000 are shown beside it.

Run from the project's top folder: python analysis/cancer-research/scattering-side/theory_numbers.py
It writes theory_numbers.txt beside itself.
"""
from pathlib import Path

import numpy as np

HERE = Path(__file__).resolve().parent
N, FAST, K, EPS, RAISE = 4000, 20.0, 1.9, 0.1, 0.005
LARGE, MU = 400_000, 3e-4
LEVELS = (1.0, 0.8, 0.7, 0.65, 0.6, 0.55, 0.5, 0.45, 0.4)
ROOT_LEVELS = (0.6, 0.55, 0.52, 0.51, 0.5, 0.45, 0.4)


def natural_rates(n):
    """The model's natural rates: evenly spaced quantiles of a Lorentzian of half-width 1 (nudge_test.natural_rates)."""
    return np.tan(np.pi * (np.arange(1, n + 1) - 0.5) / n - np.pi / 2)


def pulled_and_edge(n):
    """The pulled rates (|omega| <= 20) and b, the largest natural rate in the slower half (|omega| <= the median)."""
    w = natural_rates(n)
    wp = w[np.abs(w) <= FAST]
    a = np.abs(wp)
    return wp, a[a <= np.median(a)].max(), int((a <= np.median(a)).sum())


def faster_half_moved(wp, cf, b):
    """The faster half moved toward the edge: |omega| becomes b + cf (|omega| - b), sign kept; the slower half as is."""
    s = np.abs(wp)
    out = wp.copy()
    f = s > b
    out[f] = np.sign(wp[f]) * (b + cf * (s[f] - b))
    return out


def F(lams, ws, n, chunk=16):
    lams = np.atleast_1d(np.asarray(lams, complex))
    out = np.empty(len(lams), complex)
    for i in range(0, len(lams), chunk):
        out[i:i + chunk] = (1.0 / (lams[i:i + chunk, None] + 1j * ws[None, :])).sum(1) / n
    return out


def growing_root(k_mode, ws, n, b):
    """The fastest-growing turning root near the edge, by Newton's method from a few starting points (None if none)."""
    best = None
    for g0, d in ((0.003, 0.003), (0.01, 0.005), (0.03, 0.01)):
        lam = complex(g0, -(b + d))
        for _ in range(60):
            q = 1.0 / (lam + 1j * ws)
            val = 1 - k_mode / 2 * q.sum() / n
            der = k_mode / 2 * (q * q).sum() / n
            step = val / der
            lam -= step
            if not np.isfinite(lam) or abs(step) < 1e-13:
                break
        ok = np.isfinite(lam) and abs(1 - k_mode / 2 * F(lam, ws, n)[0]) < 1e-9
        if ok and lam.real > 0 and abs(lam.imag) > 0.5 and (best is None or lam.real > best.real):
            best = lam
    return best


def main():
    lines = ["The scattering side: the theory's expected response (a calculation from the equations; nothing simulated)", ""]
    wp, b, slower = pulled_and_edge(N)
    a = np.abs(wp)
    lines += [f"The model (N = {N} per group): {len(wp)} pulled oscillators, slower half {slower}, faster half {len(wp) - slower}.",
              f"  b, the largest natural rate in the slower half: {b:.5f}; the smallest in the faster half: {a[a > b].min():.5f}.",
              f"  The faster half's median natural rate: {np.median(a[a > b]):.3f}.", ""]

    for cf in (0.6, 0.4):
        moved = np.abs(faster_half_moved(wp, cf, b))
        lines.append(f"  The faster half moved by cf = {cf}: median {np.median(a[a > b]):.3f} -> {np.median(moved[a > b]):.3f};"
                     f" largest {a.max():.3f} -> {moved.max():.3f}.")
    lines.append("")

    WL, bL, _ = pulled_and_edge(LARGE)
    k_raised = K * (1 + RAISE) + EPS
    c = 2.0 / k_raised
    lo, hi = 1e-4, 0.05                                      # the raised pull's growing central pattern (a real root)
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        lo, hi = (lo, mid) if 1 - k_raised / 2 * F([mid], WL, LARGE)[0].real > 0 else (mid, hi)
    lines += ["Uniform narrowing, the matching estimate:",
              f"  raised pull: K + eps = {k_raised:.4f}, {k_raised / 2:.5f} times the threshold 2 Delta.",
              f"  narrowing every rate by c puts K + eps = 2.0 at 1/c times the threshold: c = {c:.6f}"
              f" (narrowing by {100 * (1 - c):.4f}%).",
              f"  For scale: the raised pull's central rhythm grows at {lo:.5f} per time unit in the large-group theory"
              f" (e-folding {1 / lo:.0f} units).", ""]

    lines += [f"The faster half moved toward the edge by cf (large-group stand-in, {LARGE:,} quantiles; b = {bL:.5f}).",
              "  centre: how crowded the centre is relative to the untouched spread (sets the threshold of the central rhythm);",
              "  J: the central rhythm's response-time factor (untouched about -0.994; a smaller size means a quicker response);",
              "  edge: how crowded the rates are just outside b, relative to the centre;",
              "  needed: the lowest K + eps at which a turning pattern near the edge could grow (the model has K + eps = 2.0).",
              ""]
    for cf in LEVELS:
        ws = faster_half_moved(WL, cf, bL)
        fr = F([2e-3, 4e-3], ws, LARGE).real
        J = (fr[1] - fr[0]) / 2e-3
        centre = F([1e-5], ws, LARGE)[0].real / F([1e-5], WL, LARGE)[0].real
        edge = (1 / (1 + bL ** 2)) / cf
        om = np.linspace(bL + 2e-4, bL + 0.5, 600)
        vals = F(MU - 1j * om, ws, LARGE)
        cross = np.where(np.sign(vals.imag[:-1]) != np.sign(vals.imag[1:]))[0]
        if len(cross):
            i = cross[0]
            needed = f"{2.0 / vals.real[i]:.3f} (at turning rate {om[i]:.4f})"
        else:
            needed = f"at least {2.0 / edge:.2f} (the pattern sits within {MU} of b, not resolved; bound from the edge crowding)"
        lines.append(f"  cf {cf:4.2f}: centre {centre:.4f}  J {J:+.4f}  edge {edge:.3f}  needed {needed}")
    lines.append("")

    for n, label in ((LARGE, "large-group stand-in"), (N, "the model's 4,000")):
        wset, bb, _ = pulled_and_edge(n)
        lines.append(f"Growing turning patterns near the edge ({label}): growth rate per time unit, turning rate, e-folding time")
        for k_mode, mode in ((K + EPS, "A and B together (K + eps = 2.0)"), (K - EPS, "A and B opposite (K - eps = 1.8)")):
            for cf in ROOT_LEVELS:
                r = growing_root(k_mode, faster_half_moved(wset, cf, bb), n, bb)
                lines.append(f"  {mode}, cf {cf:4.2f}: " + ("none growing" if r is None else
                             f"growth {r.real:.5f}, turning {abs(r.imag):.4f}, e-folding {1 / r.real:.0f} units"))
        lines.append("")
    lines.append("For the model's 4,000, growth rates below about 0.001 are near the spacing of its rates at the edge"
                 " (about 0.0008 after moving) and are not reliable.")
    text = "\n".join(lines) + "\n"
    (HERE / "theory_numbers.txt").write_text(text)
    print(text)


if __name__ == "__main__":
    main()