#!/usr/bin/env python3
"""
PROVENANCE: NUMERICAL

Computation 215 -- the native sector cannot build an apparatus: Derrick's scaling on P3's
functional, and the regulated lattice at PST's own discreteness, three orders below pinning
==============================================================================
WHAT THIS DOES.  sec:open-localisation ends on one apparatus question: PST does not derive
why T has clusters of six or more, and "the dynamical half, why a bound state forms and
persists, is inherited with the projected field theory".  The backlog carries it as "an
apparatus is a dynamical object the theory does not model".  This file asks the sharper
question underneath: could the theory's OWN field, the order parameter psi of P3 with its
Landau-Ginzburg functional (eq. LG-postulate, radially reduced: F = int [ c |grad psi|^2
- eps psi^2 + psi^4/4 ]), hold a localised static object at all?  If it could, an
apparatus might be a structure of the native sector and "inherited" would be a choice.
It cannot, and the reason is a theorem plus a number.

THE THEOREM (Derrick 1964).  For a static configuration psi(x) in d spatial dimensions
with energy E = E_grad + E_pot, E_pot = int (V(psi) - V_vac) >= 0, the rescaling
psi(lambda x) gives E(lambda) = lambda^(2-d) E_grad + lambda^(-d) E_pot, which for d = 3
is strictly decreasing in lambda: every localised configuration lowers its energy by
shrinking, there is no stationary point, and the only static finite-energy solution is
the vacuum.  Part 1 evaluates E(lambda) on three profiles (a dip, a bubble of the other
vacuum, a ring): monotone, as it must be.  V >= 0 is load-bearing (tamper ii), and the
theorem is dimensional: in d = 1 the kink is static and stable (tamper iii).

THE NUMBER.  Derrick's shrinking stops at the lattice spacing, so on the regulated theory
(comp 200: the scalar sector with spacing l_*, the grain) a localised object could survive
only by PINNING, which happens when the on-site potential dominates the coupling, i.e.
when the discreteness kappa = (m l)^2 is large.  PST's value: with the grain
l_* = hbar c / M_* and the Higgs mass, kappa = (m_h / M_*)^2 = 0.0095; with the coherence
length d_0 = sqrt 2 hbar c / m_h in place of the grain, kappa = 2 exactly.  Part 2 runs the
gradient flow of the lattice energy at kappa = 0.0095 from four localised initial data (a
bubble of radius 6 and one of radius 12, a Gaussian dip, a single flipped site): every one
relaxes to the vacuum.  Part 3 finds where pinning sets in by bisection: a single flipped
site persists only for kappa > 47.0, a four-grain bubble only for kappa > 14.7.  PST's
kappa is 2e-4 of the first threshold and 6e-4 of the second with the grain, and 1/24 and
1/7 of them with d_0.  Whichever of its two lengths is the spacing, the native sector is
in the continuum regime where Derrick rules.

THE LOOPHOLE, NAMED.  Derrick is about static states.  The 3D phi^4 theory has
oscillons: localised oscillations that keep half their core energy for 10^3 to 10^4
periods (Part 4: at t = 1000/m the core still holds 50 to 70 percent).  In PST's units a
period is 1/m_h = 5e-27 s and the core is a few 1/m_h = 1e-17 m: a breather of the order
parameter, carrying no content, lasting 1e-23 s.  Not an apparatus.

WHAT THIS CHANGES.  "Inherited" is forced, not chosen: the native sector cannot hold a
device, so an apparatus is a bound state of the inherited matter sector (the cutoff
Standard Model's own, comp 200), realised in T as a cluster, which comp 156's identity
(eq. metric: compact in M is tight in T) then reads as a device.  What stays is exactly
what the book says: why T(C) has those clusters is the configuration's content.

Tampers: (i) at twice the threshold the flipped site persists, so pinning is real and
the number is what excludes it; (ii) with the potential's sign flipped E(lambda) has a
stationary point; (iii) the 1D kink survives the same flow that erases every 3D object.

Runs in about two minutes; numpy only.
"""
import math
import numpy as np

M_H, M_STAR = 125.1, 1285.0            # GeV
KAPPA_GRAIN = (M_H / M_STAR) ** 2      # (m_h l_*)^2 with l_* = hbar c / M_*
KAPPA_D0 = 2.0                         # (m_h d_0)^2 with d_0 = sqrt 2 hbar c / m_h (comp 201)


# ------------------------------------------------------------ lattice flows
def energy(phi, kappa):
    d = phi.ndim
    g = sum(((np.roll(phi, -1, ax) - phi) ** 2).sum() for ax in range(d)) * 0.5
    return float(g + (kappa / 8) * ((phi ** 2 - 1) ** 2).sum())


def grad(phi, kappa):
    d = phi.ndim
    lap = sum(np.roll(phi, 1, ax) + np.roll(phi, -1, ax) - 2 * phi for ax in range(d))
    return -lap + (kappa / 2) * phi * (phi ** 2 - 1)


def flow(phi, kappa, steps):
    dt = 0.9 / (2 * phi.ndim + kappa)
    for _ in range(steps):
        phi = phi - dt * grad(phi, kappa)
    return phi


def radial_grid(L):
    x = np.indices((L, L, L)) - L // 2
    return np.sqrt((x ** 2).sum(0))


def flip_persists(kappa, L=12, steps=4000):
    phi = np.ones((L, L, L)); phi[L // 2, L // 2, L // 2] = -1.0
    return flow(phi, kappa, steps)[L // 2, L // 2, L // 2] < 0


def bubble_persists(kappa, R=4, L=20, steps=5000):
    phi = np.where(radial_grid(L) < R, -1.0, 1.0)
    return flow(phi, kappa, steps)[L // 2, L // 2, L // 2] < 0


def bisect(pred, lo, hi, iters=12):
    for _ in range(iters):
        mid = math.sqrt(lo * hi)
        if pred(mid):
            hi = mid
        else:
            lo = mid
    return lo, hi


# ---------------------------------------------------------- Derrick scaling
def derrick(prof, r, dr, sign=+1.0):
    dp = np.gradient(prof, dr)
    Eg = 0.5 * (dp ** 2 * 4 * np.pi * r ** 2 * dr).sum()
    Ep = sign * 0.125 * ((prof ** 2 - 1) ** 2 * 4 * np.pi * r ** 2 * dr).sum()
    lam = np.array([0.5, 0.7, 1.0, 1.4, 2.0, 3.0, 5.0])
    return Eg, Ep, lam, Eg / lam + Ep / lam ** 3          # psi(lambda x), d = 3


# --------------------------------------------------------------- oscillon
def oscillon(A, w, T=1000.0, R=300.0, dr=0.05, dt=0.02):
    n = int(R / dr); r = (np.arange(n) + 0.5) * dr
    phi = 1 - A * np.exp(-r ** 2 / (2 * w ** 2)); phi_old = phi.copy()
    damp = np.where(r > 0.8 * R, 0.05 * ((r - 0.8 * R) / (0.2 * R)) ** 2, 0.0)
    core = r < 10; out = []; E0 = None; signs = 0; last = 0.0
    for step in range(int(T / dt) + 1):
        pp = np.concatenate([[phi[0]], phi, [phi[-1]]])
        lap = (pp[2:] - 2 * pp[1:-1] + pp[:-2]) / dr ** 2 + (pp[2:] - pp[:-2]) / (2 * dr * r)
        acc = lap - 0.5 * phi * (phi ** 2 - 1) - damp * (phi - phi_old) / dt
        phi, phi_old = 2 * phi - phi_old + dt ** 2 * acc, phi
        c0 = phi[0] - 1
        if c0 * last < 0:
            signs += 1
        last = c0
        if step % int(50 / dt) == 0:
            v = (phi - phi_old) / dt; dp = np.gradient(phi, dr)
            e = (0.5 * v ** 2 + 0.5 * dp ** 2 + 0.125 * (phi ** 2 - 1) ** 2) * 4 * np.pi * r ** 2 * dr
            Ec = e[core].sum(); E0 = E0 or Ec; out.append((step * dt, Ec / E0))
    return out, signs, float(np.abs(phi[core] - 1).max())


def main():
    print("=" * 100)
    print("Computation 215: the native sector cannot build an apparatus")
    print("=" * 100)
    checks = {}
    print(f"  PST's discreteness: kappa = (m_h l_*)^2 = {KAPPA_GRAIN:.4f} with the grain; (m_h d_0)^2 = {KAPPA_D0:.1f} with the coherence length")
    print()
    # ---- Part 1: Derrick
    print("PART 1.  Derrick's scaling on three localised profiles (continuum, d = 3)")
    r = np.linspace(1e-3, 40, 40001); dr = r[1] - r[0]
    profs = {"dip": 1 - 1.5 * np.exp(-r ** 2 / 8), "bubble": np.tanh((r - 6) / 1.4), "ring": 1 - np.exp(-(r - 5) ** 2 / 2)}
    okD = True
    for name, p in profs.items():
        Eg, Ep, lam, E = derrick(p, r, dr)
        mono = bool(np.all(np.diff(E) < 0)); okD &= mono and Ep > 0
        print(f"  {name:7s}: E_grad = {Eg:8.3f}, E_pot = {Ep:8.3f}; E(lambda) for lambda = 0.5 .. 5: {np.round(E, 2)}; strictly decreasing: {mono}")
    checks['derrick'] = okD
    print()
    # ---- Part 2: flows at PST's kappa
    print("PART 2.  gradient flow of the lattice energy at kappa = 0.0095 from four localised initial data (32^3)")
    L = 32; rr = radial_grid(L); okF = True
    inits = {"bubble R=6": np.where(rr < 6, -1.0, 1.0), "bubble R=12": np.where(rr < 12, -1.0, 1.0),
             "gaussian dip": 1 - 1.5 * np.exp(-rr ** 2 / 32), "single-site flip": None}
    for name, init in inits.items():
        phi = init.copy() if init is not None else np.ones((L, L, L))
        if init is None:
            phi[L // 2, L // 2, L // 2] = -1.0
        E0 = energy(phi, KAPPA_GRAIN); path = []
        for _ in range(6):
            phi = flow(phi, KAPPA_GRAIN, 500); path.append(float(np.abs(phi - 1).max()))
        E1 = energy(phi, KAPPA_GRAIN); gone = path[-1] < 2e-2 and E1 < 1e-2 * E0
        okF &= gone
        print(f"  {name:17s}: E 0 -> 3000 steps: {E0:9.3f} -> {E1:.2e}; max|phi - 1| along the flow {[f'{x:.3f}' for x in path]}; vacuum reached: {gone}")
    checks['flows'] = okF
    print()
    # ---- Part 3: pinning thresholds
    print("PART 3.  where pinning sets in (bisection in kappa)")
    lo1, hi1 = bisect(flip_persists, 1.0, 1000.0)
    lo2, hi2 = bisect(bubble_persists, 0.5, 400.0)
    print(f"  a single flipped site persists for kappa > {lo1:.2f} (bracket [{lo1:.2f}, {hi1:.2f}]); PST: {KAPPA_GRAIN / lo1:.1e} of it with the grain, {KAPPA_D0 / lo1:.3f} with d_0")
    print(f"  a four-grain bubble persists for kappa > {lo2:.2f} (bracket [{lo2:.2f}, {hi2:.2f}]); PST: {KAPPA_GRAIN / lo2:.1e} of it with the grain, {KAPPA_D0 / lo2:.3f} with d_0")
    checks['pinning'] = 40 < lo1 < 60 and 10 < lo2 < 20 and KAPPA_D0 < lo2 and KAPPA_GRAIN < 1e-3 * lo1
    print()
    # ---- Part 4: the oscillon
    print("PART 4.  the time-dependent loophole: a radial oscillon (m = 1 units, t to 1000)")
    okO = True
    hbar_s = 6.582e-25 / M_H
    for A, w in ((1.5, 3.0), (1.0, 4.5)):
        out, signs, amp = oscillon(A, w)
        fr = {t: f for t, f in out}
        okO &= 0.3 < fr[1000.0] < 0.9 and signs > 50 and amp > 0.1
        print(f"  A = {A}, w = {w}: core energy fraction at t = 100, 300, 600, 1000: {fr[100.0]:.2f}, {fr[300.0]:.2f}, {fr[600.0]:.2f}, {fr[1000.0]:.2f}; the core crosses the vacuum {signs} times; amplitude at the end {amp:.2f}")
    print(f"  in PST's units: a period 1/m_h = {hbar_s:.1e} s, a core of a few 1/m_h = {3 * 197.3e-18 / M_H:.1e} m, a life of 10^3 to 10^4 periods = 1e-23 s: a breather of the order parameter, not an apparatus")
    checks['oscillon'] = okO
    print()
    # ---- Tampers
    print("TAMPERS")
    t1 = flip_persists(2 * hi1)
    print(f"  (i)   at kappa = {2 * hi1:.1f}, twice the threshold, the flipped site persists: {t1}")
    Eg, Ep, lam, E = derrick(profs['dip'], r, dr, sign=-1.0)
    t2 = not bool(np.all(np.diff(E) < 0))
    print(f"  (ii)  Derrick with the potential's sign flipped: E(lambda) = {np.round(E, 2)}, no longer monotone: {t2}")
    x = np.arange(200) - 100; kink = np.tanh(0.5 * x * math.sqrt(KAPPA_D0) / 2)
    k1 = flow(kink.copy(), KAPPA_D0, 3000); t3 = k1[100 + 20] > 0.5 and k1[100 - 20] < -0.5 and abs(k1[100]) < 0.05
    print(f"  (iii) a 1D kink at kappa = 2 survives 3000 flow steps (phi(-20), phi(0), phi(+20) = {k1[80]:+.3f}, {k1[100]:+.3f}, {k1[120]:+.3f}): Derrick is dimensional: {t3}")
    checks['tampers'] = t1 and t2 and t3
    print()
    ok = all(checks.values())
    print("SUMMARY")
    print(f"  [{'PASS' if checks['derrick'] else 'FAIL'}] Derrick: every localised configuration of P3's functional lowers its energy by shrinking; no static lump in d = 3")
    print(f"  [{'PASS' if checks['flows'] else 'FAIL'}] on the regulated lattice at PST's discreteness every localised initial datum relaxes to the vacuum")
    print(f"  [{'PASS' if checks['pinning'] else 'FAIL'}] pinning needs kappa above 47 (a site) or 15 (a four-grain bubble); PST is below both with either of its lengths")
    print(f"  [{'PASS' if checks['oscillon'] else 'FAIL'}] the one long-lived localised object is an oscillon: a breather at 1/m_h, not a device")
    print(f"  [{'PASS' if checks['tampers'] else 'FAIL'}] tampers rejected")
    print()
    print("  CONCLUSION: the order parameter cannot hold a localised static object, by Derrick's theorem in the")
    print("  continuum and by three orders of discreteness short of pinning on the regulated lattice, so an")
    print("  apparatus is a bound state of the inherited matter sector, realised in T as a cluster that comp 156's")
    print("  identity reads as a device.  'Inherited' is forced, not chosen; what the configuration supplies is")
    print("  which clusters, and that stays T(C)-contingent.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()