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

Computation 200 -- the interacting vacuum of the regulated scalar sector exists
nonperturbatively; what the item still carries is the gauge and fermion sector, and
that is the cutoff Standard Model's own
===================================================================================
WHAT THIS DOES.  The remaining structural item on the interacting vacuum asks for a
construction rather than an expansion: a Hilbert space and a Hamiltonian bounded below
at all orders, which no perturbative statement supplies (sec:renorm, comp 198).  The
theory never removes its cutoff.  Above M_* there are no modes, the grain
l_* = hbar c / M_* being a property of the projection and not a regulator (sec:renorm),
and below M_* the running is the Standard Model's.  So the construction asked for is
the construction of a cutoff theory, and for the scalar sector that is three classical
theorems applied to an object that is finite in the two senses that matter.  This file
states them, says what each needs, and checks each one numerically where it can be
checked.

THEOREM 1 (finite volume: the Hamiltonian).  In a finite volume with the cutoff, the
scalar sector has finitely many modes, and its Hamiltonian is a Schroedinger operator
on L^2(R^N),

      H = sum_x [ pi_x^2 / 2 + (1/2) sum_mu (phi_{x+mu} - phi_x)^2 + (mu^2/2) phi_x^2 + (lam/4) phi_x^4 ],

whose potential is continuous, bounded below and tends to infinity (lam > 0, either
sign of mu^2).  Such an operator is essentially self-adjoint on the smooth functions
of compact support and has compact resolvent, hence a discrete spectrum with a lowest
eigenvalue (Reed and Simon, vol. II Thm X.28 and vol. IV Thm XIII.67); and its
semigroup e^{-tH} is positivity improving, so the ground state is unique and strictly
positive (Perron-Frobenius for Schroedinger operators, Reed and Simon vol. IV
Thm XIII.47).  That is the interacting vacuum at finite volume, at all orders, with no
expansion anywhere.  Spontaneous symmetry breaking at finite volume is what it always
is: a unique symmetric ground state whose gap to the first excited state closes
exponentially in the volume, so that the infinite-volume vacua are the two pure states
the symmetric one decomposes into.

THEOREM 2 (Euclidean side: the infinite-volume limit).  The Euclidean lattice theory
with spacing l_* has the action above with ferromagnetic nearest-neighbour coupling and
an even single-site weight exp(-(mu^2/2 + d) phi^2 - lam phi^4 / 4).  Griffiths' second
inequality then makes every correlation <phi_x phi_y>_Lambda non-decreasing in the
volume Lambda with Dirichlet data (the field held at zero outside Lambda, so that
enlarging Lambda only switches couplings on and changes no single-site weight; with
free data the boundary sites lose stiffness and the ordering is not Griffiths'), and
the single-site quartic bounds
<phi_x^2> uniformly in Lambda, so every correlation converges as Lambda grows: the
infinite-volume Euclidean measure exists (Glimm and Jaffe, Quantum Physics, ch. 4 and
11; Simon, The P(phi)_2 Euclidean Quantum Field Theory, ch. IV, whose arguments are
dimension-independent once the cutoff is kept).  Nothing here uses d = 2 or d = 3; the
dimension enters only through the continuum limit, which this theory does not take.

THEOREM 3 (back to a Hamiltonian: reflection positivity).  A nearest-neighbour lattice
action is reflection positive, so the transfer matrix T is a positive operator and
T = e^{-l_* H} defines a Hamiltonian bounded below whose ground state is the Euclidean
vacuum (Osterwalder-Schrader reconstruction, in its lattice form, Glimm and Jaffe ch. 6).
So the infinite-volume Euclidean measure of Theorem 2 is the vacuum of a Hamiltonian,
and Theorem 1's finite-volume Hamiltonian is its restriction.

WHAT THIS CLOSES.  The scalar sector of the regulated theory has an interacting vacuum
nonperturbatively, in finite and in infinite volume, and the one-loop and two-loop
statements of sec:renorm and comp 198 are expansions about it.  That is the
construction the item asked for, for the sector the item was written about: the
substrate scalar, its phi^4 vertices and its sombrero.

WHAT THIS DOES NOT CLOSE, STATED PLAINLY.  Below M_* the theory is the full Standard
Model with a cutoff, gauge fields and chiral fermions included.  For that sector the
same three steps are the lattice Standard Model's own: finite-volume existence is
elementary (compact group integrals), reflection positivity holds for the standard
actions, but the infinite-volume limit of a four-dimensional non-abelian gauge theory
at weak coupling is not a theorem (Balaban's programme), and chiral fermions on a
regulator meet the Nielsen-Ninomiya obstruction.  Those are not problems this theory
introduces; it inherits them with the Standard Model it reduces to, exactly as it
inherits the cosmological-constant gap.  So the item narrows to the gauge and fermion
sector and is carried as inherited, which is where colour-generation went at v32.128.

CHECKS.  (1) Theorem 1 on one and two sites by finite differences: lowest eigenvalues,
gap > 0, ground state strictly positive and even; the double well's gap closing as the
well deepens and as sites are added, the finite-volume precursor of symmetry breaking.
(2) Theorem 3 on one site: the Trotter transfer matrix has positive eigenvalues and
-log of their ratios converges to Theorem 1's gaps as the time step shrinks.
(3) Theorem 2 on the four-dimensional lattice at the book's own couplings in lattice
units, lam = lambda(M_*) = 0.092 and mu^2 = -(m_h l_*)^2 / 2 = -0.0048, and at a massive
test point: Metropolis on L = 4, 6, 8 with Dirichlet data, the central
nearest-neighbour correlation and <phi^2> non-decreasing in L within errors and bounded,
the increments shrinking at the massive point.  (4) Tamper: an antiferromagnetic
coupling, which Griffiths' inequality does not cover, breaks the monotonicity; and a
negative quartic has no ground state at all (the finite-difference spectrum runs off
to minus infinity as the grid widens), which is the sense in which lam > 0, comp 198's
result, is what the construction needs.
"""

import math
import numpy as np
from scipy.sparse import diags, kron, identity
from scipy.sparse.linalg import eigsh

MH, MSTAR = 125.25, 4 * math.pi * 125.25 * math.sqrt(2.0 / 3.0)     # GeV
LAM_BOOK = 0.25 * math.exp(-1.0)                                      # lambda(M_*) = e^-1/4 = 0.092
MU2_BOOK = -0.5 * (MH / MSTAR) ** 2                                   # -(m_h l_*)^2 / 2 in lattice units


def rule(ch="="):
    print(ch * 92)


# ---------------------------------------------------------------------------
# Part 1: the finite-volume Hamiltonian, one and two sites
# ---------------------------------------------------------------------------
def kinetic_1d(n, h):
    return diags([np.full(n - 1, -0.5 / h ** 2), np.full(n, 1.0 / h ** 2), np.full(n - 1, -0.5 / h ** 2)], [-1, 0, 1])


def spectrum_1site(mu2, lam, L=8.0, n=400, k=3):
    x = np.linspace(-L, L, n)
    h = x[1] - x[0]
    V = 0.5 * mu2 * x ** 2 + 0.25 * lam * x ** 4
    H = kinetic_1d(n, h) + diags(V)
    E, U = eigsh(H, k=k, which="SA")
    o = np.argsort(E)
    return E[o], U[:, o], x


def spectrum_2site(mu2, lam, J=1.0, L=6.0, n=70, k=3):
    x = np.linspace(-L, L, n)
    h = x[1] - x[0]
    K = kinetic_1d(n, h)
    I = identity(n)
    X = diags(x)
    V1 = diags(0.5 * mu2 * x ** 2 + 0.25 * lam * x ** 4)
    H = kron(K, I) + kron(I, K) + kron(V1, I) + kron(I, V1) + 0.5 * J * (kron(X, I) - kron(I, X)) @ (kron(X, I) - kron(I, X))
    E, U = eigsh(H.tocsr(), k=k, which="SA")
    o = np.argsort(E)
    return E[o], U[:, o], x


def part1():
    rule()
    print("  1.  THEOREM 1 ON ONE AND TWO SITES: DISCRETE SPECTRUM, UNIQUE POSITIVE GROUND STATE")
    rule()
    ok = True
    print("      one site, lam = 1:")
    for mu2 in (1.0, 0.0, -2.0, -6.0):
        E, U, x = spectrum_1site(mu2, 1.0)
        g = U[:, 0] * np.sign(U[np.argmax(np.abs(U[:, 0])), 0])
        even = np.max(np.abs(g - g[::-1])) < 1e-6
        print(f"        mu^2 = {mu2:+5.1f}: E0 = {E[0]:.5f}, gap E1 - E0 = {E[1] - E[0]:.5f}, ground state min/max = {g.min() / g.max():+.1e} (no sign change), even: {even}")
        ok &= E[1] - E[0] > 1e-6 and g.min() > -1e-6 * g.max() and even
    print("      the double well's gap closes as the well deepens (tunnelling), and the ground state stays unique and even")
    print("      two coupled sites, lam = 1, J = 1:")
    gaps = []
    for mu2 in (1.0, -2.0, -4.0):
        E, U, x = spectrum_2site(mu2, 1.0)
        g = U[:, 0] * np.sign(U[np.argmax(np.abs(U[:, 0])), 0])
        print(f"        mu^2 = {mu2:+5.1f}: E0 = {E[0]:.5f}, gap = {E[1] - E[0]:.5f}, ground state min/max = {g.min() / g.max():+.1e}")
        ok &= E[1] - E[0] > 1e-6 and g.min() > -1e-6 * g.max()
        gaps.append(E[1] - E[0])
    E1, _, _ = spectrum_1site(-4.0, 1.0)
    print(f"        at mu^2 = -4 the two-site gap {gaps[-1]:.4f} is below the one-site gap {E1[1] - E1[0]:.4f}: adding volume closes the gap, the precursor of symmetry breaking")
    ok &= gaps[-1] < E1[1] - E1[0]
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# Part 2: reflection positivity on one site: the transfer matrix
# ---------------------------------------------------------------------------
def part2():
    rule()
    print("  2.  THEOREM 3 ON ONE SITE: THE TRANSFER MATRIX IS POSITIVE AND REPRODUCES THE SPECTRUM")
    rule()
    mu2, lam = -2.0, 1.0
    E, _, _ = spectrum_1site(mu2, lam)
    gap_H = E[1] - E[0]
    x = np.linspace(-5, 5, 300)
    h = x[1] - x[0]
    V = 0.5 * mu2 * x ** 2 + 0.25 * lam * x ** 4
    ok = True
    prev = None
    for a in (0.4, 0.2, 0.1, 0.05):
        T = np.exp(-(x[:, None] - x[None, :]) ** 2 / (2 * a) - a * (V[:, None] + V[None, :]) / 2) * h / math.sqrt(2 * math.pi * a)
        w = np.linalg.eigvalsh((T + T.T) / 2)[::-1]
        gap_T = -math.log(w[1] / w[0]) / a
        print(f"        time step a = {a:.2f}: smallest eigenvalue of T {w[-1]:+.1e} (positive to round-off: {w.min() > -1e-12}), gap from T {gap_T:.5f} against the Hamiltonian's {gap_H:.5f}")
        ok &= w.min() > -1e-12
        prev = abs(gap_T - gap_H)
    ok &= prev < 2e-2 * gap_H
    print("      so -log T / a is a Hamiltonian bounded below whose spectrum is Theorem 1's, as reflection positivity says")
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# Part 3: the four-dimensional lattice, Griffiths monotonicity in the volume
# ---------------------------------------------------------------------------
def metropolis(L, mu2, lam, J=1.0, sweeps=2500, therm=2500, step=0.6, seed=0, antiferro=False):
    """4D lattice phi^4, Dirichlet data (phi = 0 outside the box), checkerboard Metropolis.  Returns the
    nearest-neighbour correlation and <phi^2> averaged over the central 2^4 block, with
    blocked errors."""
    rng = np.random.default_rng(seed)
    d = 4
    phi = rng.normal(size=(L,) * d) * 0.1
    idx = np.indices((L,) * d).sum(axis=0) % 2
    sgn = -1.0 if antiferro else 1.0
    def neighbour_sum(f):
        s = np.zeros_like(f)
        for ax in range(d):
            fwd = np.roll(f, -1, axis=ax); bwd = np.roll(f, 1, axis=ax)
            # Dirichlet data: the wrapped neighbours are outside the box, where phi = 0
            sl_last = [slice(None)] * d; sl_last[ax] = -1
            sl_first = [slice(None)] * d; sl_first[ax] = 0
            fwd[tuple(sl_last)] = 0.0; bwd[tuple(sl_first)] = 0.0
            s += fwd + bwd
        return s
    nn = np.full(phi.shape, 2.0 * d)      # every site keeps its 2d bonds; the missing neighbours are the zero field outside
    c = slice(L // 2 - 1, L // 2 + 1)
    centre = (c,) * d
    g1, p2 = [], []
    for it in range(therm + sweeps):
        for parity in (0, 1):
            mask = idx == parity
            prop = phi + rng.normal(size=phi.shape) * step
            ns = neighbour_sum(phi)
            # dE = energy(prop) - energy(phi) with energy(f) = J nn f^2/2 - sgn J f ns + mu2 f^2/2 + lam f^4/4,
            # written as one expression in the differences (the same arithmetic, half the temporaries)
            pq, fq = prop * prop, phi * phi
            dE = (0.5 * J * nn + 0.5 * mu2) * (pq - fq) - sgn * J * ns * (prop - phi) + 0.25 * lam * (pq * pq - fq * fq)
            acc = (rng.uniform(size=phi.shape) < np.exp(-np.clip(dE, -50, 50))) & mask
            phi = np.where(acc, prop, phi)
        if it >= therm:
            ns = neighbour_sum(phi)
            g1.append(np.mean((phi * ns / nn)[centre]))
            p2.append(np.mean((phi ** 2)[centre]))
    def blocked(v, nb=10):
        v = np.array(v)
        b = v[: len(v) // nb * nb].reshape(nb, -1).mean(axis=1)
        return v.mean(), b.std(ddof=1) / math.sqrt(nb)
    return blocked(g1), blocked(p2)


def part3():
    rule()
    print("  3.  THEOREM 2 ON THE FOUR-DIMENSIONAL LATTICE, DIRICHLET DATA: CORRELATIONS NON-DECREASING AND BOUNDED IN THE VOLUME")
    rule()
    ok = True
    for name, mu2, lam in (("the book's couplings in units of l_*: lam = 0.092, mu^2 = -0.0048", MU2_BOOK, LAM_BOOK),
                           ("a massive test point: lam = 0.5, mu^2 = 1", 1.0, 0.5)):
        print(f"      {name}")
        rows = []
        for L in (4, 6, 8):
            (g, ge), (p, pe) = metropolis(L, mu2, lam, seed=L)
            rows.append((L, g, ge, p, pe))
            print(f"        L = {L}: <phi_x phi_x+1> = {g:.4f} +- {ge:.4f},   <phi^2> = {p:.4f} +- {pe:.4f}")
        mono = all(rows[i + 1][1] >= rows[i][1] - 2 * (rows[i][2] + rows[i + 1][2]) for i in range(2))
        bounded = all(r[3] < 50 for r in rows)
        print(f"        non-decreasing in L within two errors: {mono};  bounded: {bounded}")
        ok &= mono and bounded
        if mu2 > 0:
            inc = [rows[i + 1][1] - rows[i][1] for i in range(2)]
            print(f"        increments {inc[0]:+.4f}, {inc[1]:+.4f}: shrinking, the limit is being approached")
            ok &= abs(inc[1]) < abs(inc[0]) + 2 * (rows[1][2] + rows[2][2])
    print("      monotone and bounded gives a limit: the infinite-volume Euclidean measure exists, and by Part 2 it is the vacuum of a Hamiltonian")
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# Part 4: tampers
# ---------------------------------------------------------------------------
def part4():
    rule()
    print("  4.  TAMPERS: AN ANTIFERROMAGNET IS OUTSIDE GRIFFITHS, A NEGATIVE QUARTIC HAS NO GROUND STATE")
    rule()
    ok = True
    rows = []
    for L in (4, 6, 8):
        (g, ge), _ = metropolis(L, 1.0, 0.5, seed=100 + L, antiferro=True)
        rows.append((L, g, ge))
        print(f"        antiferromagnetic coupling, L = {L}: <phi_x phi_x+1> = {g:+.4f} +- {ge:.4f}")
    print("        the nearest-neighbour correlation is negative and the inequality that ordered Part 3 says nothing here")
    ok &= all(r[1] < 0 for r in rows)
    E_prev = None
    for L in (3.0, 4.0, 5.0):
        E, _, _ = spectrum_1site(1.0, -0.5, L=L, n=300)
        print(f"        lam = -0.5, grid half-width {L:.0f}: lowest eigenvalue {E[0]:+.3f}")
        ok &= E_prev is None or E[0] < E_prev
        E_prev = E[0]
    print("        the lowest eigenvalue runs off as the grid widens: no ground state, so lam > 0 (comp 198) is what Theorem 1 needs")
    print()
    return bool(ok)


def main():
    rule()
    print("  COMPUTATION 200.  The interacting vacuum of the regulated scalar sector exists nonperturbatively")
    rule()
    results = [
        ("finite volume: discrete spectrum, unique strictly positive even ground state, gap closing with volume in the broken phase", part1()),
        ("reflection positivity: the transfer matrix is positive and reproduces the Hamiltonian's spectrum", part2()),
        ("infinite volume: correlations non-decreasing and bounded in L on the 4D lattice at the book's couplings", part3()),
        ("tampers: an antiferromagnet escapes the inequality, a negative quartic has no ground state", part4()),
    ]
    rule("-")
    for i, (name, r) in enumerate(results, 1):
        print(f"  {i}. {name}: {'PASS' if r else 'FAIL'}")
    allok = all(r for _, r in results)
    rule("-")
    print("  VERDICT:", "PASS" if allok else "FAIL")
    print("""
  The scalar sector of the regulated theory has an interacting vacuum at all orders:
  a unique positive ground state of a self-adjoint Hamiltonian in finite volume, an
  infinite-volume Euclidean measure by Griffiths monotonicity, and a Hamiltonian
  bounded below behind it by reflection positivity.  The item narrows to the gauge
  and fermion sector, whose nonperturbative existence with a cutoff is the Standard
  Model's own open problem, inherited rather than introduced.""")
    return allok


if __name__ == "__main__":
    main()