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

Computation 207 -- the spectral dimension of the projected geometry: which way it flows
=======================================================================================
WHAT THIS DOES.  v32.144 withdrew the claim that PST's scale-dependent dimensionality
agrees with Causal Dynamical Triangulations (CDT): CDT's spectral dimension rises from
two at short scales to four, the book's "high near threshold, four far from it" runs
the other way, and the two quantities (a count of directions, a diffusion exponent) were
never shown to be the same thing.  The section then said the spectral dimension of the
projected geometry had not been computed.  This file computes it, exactly, in the finite
models the Mosco chapter already uses, and finds in which regime it rises, in which it
sits at the Mosco value, and in which it falls.

THE MODEL (the book's own objects, nothing added).
  Configuration space P(D) = {0,1}^|D| with the Bernoulli(1/2) measure; the Boolean
  Dirichlet form eq. boolean-dirichlet, whose generator is the hypercube Laplacian: each
  bit flips independently.  The Mosco rescaling (1/|D|) E^B makes the rate per bit 1/|D|,
  so the configuration walks at one flip per unit of Mosco time s.  The realisation map
  pushes the walk to M: flipping bit a displaces the point by xi(a), with the A4
  covariance (d_0^2/4) g^{mu nu}, which is the i.i.d. surrogate of comps 70 and 83 (flat
  bulk, Gaussian xi).  By time s the bits flipped an odd number of times number
  K ~ Binomial(|D|, p), p = (1 - exp(-2s/|D|))/2, and the displacement given K is
  Gaussian with covariance K (d_0^2/4) I_d, whichever bits they were.  The return
  probability at resolution eps (the walker counts as returned if it is within eps of its
  start) is therefore the exact finite mixture

      P_eps(s) = sum_K Binom(K; |D|, p) * chi2_d.cdf( 4 eps^2 / (K d_0^2) ),   K = 0 term 1,

  and the spectral dimension is Carlip's d_s(s) = -2 d log P / d log s.  Two resolutions
  are physical: eps = d_0, the finest resolution at which rho distinguishes spacetime
  points (sec:grain-vs-coherence), and eps = ell_* = d_0/14.3, the grain (canonical
  normalisation, d_0 = 2.2 am, ell_* = 1.5e-19 m).  eps -> 0 is exact return on the
  hypercube, P = ((1 + exp(-2s/|D|))/2)^|D|.  Controls: the continuous-time walk on the
  cubic lattice Z^d at the same total rate, P = e^{-s} I_0(s/d)^d, valency 2d fixed.

RESULTS (d = 4, the dimension of M; d = 3 for the spatial surrogate of comp 70 agrees).
  1. THE MOSCO PLATEAU IS THERE AND WIDENS WITH |D|.  At every resolution d_s sits
     within 0.1 of 4 on a window of Mosco time whose lower end is fixed (s = 2 at the
     grain, s = 5 at d_0) and whose upper end grows with |D|: 25, 55, 144, 464, 1690 for
     |D| = 256 to 65536, the ratio per quadrupling tending to 4, so the plateau's top end
     grows in proportion to |D|, where bits start to be flipped twice and the walk
     saturates.  This is the Mosco limit seen as a diffusion exponent.
  2. AT THE RESOLUTION d_0 THE SPECTRAL DIMENSION DOES NOT RISE.  With eps = d_0 the
     curve climbs from 0 at s -> 0 through 2 near s = 2 to a maximum of 4.46 at s = 9 and
     settles at 4: an overshoot smaller than the cubic lattice's own (Z^4 peaks at 4.87,
     Z^3 at 3.65).  Below the coherence length the projected geometry is, for a probe that
     resolves d_0 and no finer, a discrete structure like any other: the dimension falls
     to zero, not to two.
  3. AT THE GRAIN IT RISES, AND THE RISE IS THE CAUSAL SET MECHANISM.  With eps = ell_*
     = d_0/14.3 the curve climbs through 4 at s = 2.5, reaches 22.9 at s = 13 (a few
     coherence lengths; the rms displacement at Mosco time s is sqrt(s) d_0) and falls
     back to the plateau by s = 30.  The height is a logarithm of the resolution: 4.5,
     6.8, 13.9, 19.9, 25.8, 33.5, 39.3, 49.2, 58.1 for eps/d_0 = 1, 1/2, 1/5, 1/10, 1/20,
     1/50, 1/100, 1/300, 1/1000, about 8.4 per e-fold, unbounded as eps -> 0, where the
     exact-return curve d_s = 2s (bare hypercube, no dimension at all, as comp 4 noted)
     is recovered.  The mechanism is the one Eichhorn, Surya and Versteegen identified on
     causal sets: a walker's first step can go in any of |D| directions, so the return
     probability drops exponentially rather than as a power until the diffusion has
     spread past the resolution, and the log-derivative in that stretch is large.
     Restricting the causal set diffusion to a hypersurface bounds the valency and
     restores reduction; here, resolving at d_0 does the same, because steps inside eps
     count as returns.
  4. THE FIXED-VALENCY CONTROL IS BOUNDED.  Z^4's overshoot is 4.87 at every resolution
     (returns on a lattice are exact), Z^3's 3.65.  The rise in the projected geometry
     below d_0 is therefore the unbounded-valency signature, not a generic discreteness
     artefact.  CDT's two appears in neither model: dimensional reduction to two is a
     property of a sum over geometries, not of a single graph's diffusion.

WHAT THIS CHANGES.  The sentence withdrawn at v32.144 can be restated as a computed one:
between the grain and the coherence length the spectral dimension of the projected
geometry rises above four, the direction of the full-causal-set result and against CDT,
by an amount logarithmic in d_0/ell_* (about 23 at the canonical normalisation); above
d_0 it is four over a window that grows with |D|; at the resolution d_0 itself it never
exceeds 4.5.  The book's direction count ("fewer directions as epsilon(C) grows") and
the diffusion exponent flow the same way in the sub-d_0 regime, which is the regime the
Conclusion's withdrawn sentence was about, and remain different quantities.

Tamper: replacing the Binomial number of odd-flipped bits by its mean (a deterministic
K = |D| p) removes the K = 0 term and with it the whole rise; the curve then starts at 4.
The rise lives in the probability of not having moved at all.
"""
import numpy as np
from scipy.special import gammaln, gammainc, ive


def log_binom(D, K, p):
    K = np.asarray(K, float)
    with np.errstate(divide="ignore"):
        return gammaln(D + 1) - gammaln(K + 1) - gammaln(D - K + 1) + K * np.log(p) + (D - K) * np.log1p(-p)


def log_P_push(s, D, d, eps):
    """log of the return probability within eps (units of d_0) of the pushed-forward walk."""
    if eps <= 0:
        return D * np.log((1 + np.exp(-2 * s / D)) / 2)
    p = (1 - np.exp(-2 * s / D)) / 2
    K = np.arange(0, D + 1)
    lw = log_binom(D, K, p)
    ret = np.ones(D + 1)
    ret[1:] = gammainc(d / 2, 2 * eps ** 2 / K[1:])
    m = lw.max()
    return m + np.log(np.sum(np.exp(lw - m) * ret))


def log_P_lattice(s, d):
    """continuous-time walk on Z^d, total rate 1: P = e^{-s} I_0(s/d)^d."""
    return d * np.log(ive(0, s / d))


def ds_curve(logP, S):
    lp = np.array([logP(s) for s in S])
    return -2 * np.gradient(lp, np.log(S))


def plateau(S, ds, d, tol=0.1):
    w = S[np.abs(ds - d) < tol]
    return (w.min(), w.max()) if len(w) else (np.nan, np.nan)


def main():
    S = np.logspace(-2, 5, 400)
    d = 4
    D0 = 4096
    ELL = 1 / 14.3
    print("Computation 207 -- the spectral dimension of the projected geometry")
    print("=" * 72)
    print(f"d = {d}, Mosco time s (one bit flip per unit), lengths in units of d_0; |D| = {D0} unless stated")
    print()
    curves = {}
    print("1. the curves d_s(s) at the three resolutions and for the exact return")
    for name, eps in (("eps = d_0", 1.0), ("eps = ell_* = d_0/14.3", ELL), ("eps = d_0/100", 0.01), ("exact return", 0.0)):
        ds = ds_curve(lambda s: log_P_push(s, D0, d, eps), S)
        curves[name] = ds
        lo, hi = plateau(S, ds, d)
        i = int(np.argmax(ds[S < 1000]))
        print(f"   {name:24s}: max d_s = {ds[S < 1000][i]:6.2f} at s = {S[S < 1000][i]:7.3g};  plateau |d_s - 4| < 0.1 on s in [{lo:.3g}, {hi:.3g}]")
        print("      s   : " + "  ".join(f"{s:7.3g}" for s in (0.1, 1, 2, 3, 5, 10, 13, 20, 30, 100, 1000, 4096, 20000)))
        print("      d_s : " + "  ".join(f"{ds[np.argmin(np.abs(S - s))]:7.2f}" for s in (0.1, 1, 2, 3, 5, 10, 13, 20, 30, 100, 1000, 4096, 20000)))
    ok_d0 = curves["eps = d_0"].max() < 4.6 and curves["eps = d_0"].max() > 4.0
    ok_ell = 20 < curves["eps = ell_* = d_0/14.3"][S < 1000].max() < 26
    ok_exact = abs(curves["exact return"][np.argmin(np.abs(S - 100))] / (2 * 100) - 1) < 0.1
    print(f"   at eps = d_0 the maximum is below 4.6 (no rise beyond a lattice-sized overshoot): {ok_d0}")
    print(f"   at the grain the maximum lies between 20 and 26: {ok_ell}")
    print(f"   the exact-return curve is the bare hypercube's d_s = 2s (checked at s = 100): {ok_exact}")
    print()
    print("2. the Mosco plateau widens with |D| (eps = ell_*)")
    tops = []
    for D in (256, 1024, 4096, 16384, 65536):
        ds = ds_curve(lambda s: log_P_push(s, D, d, ELL), S)
        lo, hi = plateau(S, ds, d)
        tops.append(hi)
        print(f"   |D| = {D:6d}: plateau on s in [{lo:.3g}, {hi:.4g}]   (max d_s {ds[S < 1000].max():.2f})")
    ratios = [tops[i + 1] / tops[i] for i in range(len(tops) - 1)]
    print("   ratio of the top end per quadrupling of |D|: " + ", ".join(f"{r:.2f}" for r in ratios) + "  (tending to 4)")
    ok_plateau = all(r > 1.5 for r in ratios) and ratios[-1] > ratios[0] and tops[-1] > 1000
    print(f"   the top end grows without bound, in proportion to |D| asymptotically: {ok_plateau}")
    print()
    print("3. the height of the rise against the resolution (|D| = 16384)")
    hs = []
    for eps in (1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.003, 0.001):
        ds = ds_curve(lambda s: log_P_push(s, 16384, d, eps), S)
        h = ds[S < 1000].max(); hs.append((eps, h))
        print(f"   eps/d_0 = {eps:6.3f}: max d_s = {h:6.2f} at s = {S[S < 1000][np.argmax(ds[S < 1000])]:5.3g}")
    x = np.log(1 / np.array([e for e, _ in hs[2:]])); y = np.array([h for _, h in hs[2:]])
    slope = np.polyfit(x, y, 1)[0]
    print(f"   growth per e-fold of d_0/eps over eps <= d_0/5: {slope:.2f}; unbounded as eps -> 0")
    ok_log = 6 < slope < 11 and all(hs[i + 1][1] > hs[i][1] for i in range(len(hs) - 1))
    print(f"   the height is monotone in the resolution and grows logarithmically: {ok_log}")
    print()
    print("4. the fixed-valency control")
    lat = {}
    for dd in (3, 4):
        ds = ds_curve(lambda s: log_P_lattice(s, dd), S)
        lat[dd] = ds.max()
        print(f"   Z^{dd}: max d_s = {ds.max():.3f} at s = {S[np.argmax(ds)]:.3g}; at s = 1000, d_s = {ds[np.argmin(np.abs(S - 1000))]:.3f}")
    ok_lat = lat[4] < 5 and lat[3] < 4 and curves["eps = d_0"].max() < lat[4]
    print(f"   bounded overshoots, resolution-independent; the eps = d_0 curve stays below Z^4's: {ok_lat}")
    print()
    print("5. tamper: deterministic K = |D| p (no probability of not having moved)")
    def log_P_det(s):
        K = max(D0 * (1 - np.exp(-2 * s / D0)) / 2, 1e-12)
        return np.log(gammainc(d / 2, 2 * ELL ** 2 / K))
    ds_t = ds_curve(log_P_det, S)
    print(f"   max d_s over s < 1000 = {ds_t[S < 1000].max():.2f}: the rise is gone; the curve starts at the plateau")
    ok_t = ds_t[S < 1000].max() < 4.3
    print()
    ok = ok_d0 and ok_ell and ok_exact and ok_plateau and ok_log and ok_lat and ok_t
    print("SUMMARY")
    print(f"  [{'PASS' if ok_plateau else 'FAIL'}] the Mosco plateau at d_s = 4 is there and its top end grows with |D|")
    print(f"  [{'PASS' if ok_d0 else 'FAIL'}] at the resolution d_0 the spectral dimension never exceeds 4.5: no rise")
    print(f"  [{'PASS' if ok_ell else 'FAIL'}] at the grain ell_* it rises to about 23 before settling at four")
    print(f"  [{'PASS' if ok_log and ok_exact else 'FAIL'}] the rise is logarithmic in d_0/eps and unbounded for exact return")
    print(f"  [{'PASS' if ok_lat else 'FAIL'}] the fixed-valency lattice's overshoot is bounded: the rise is the valency signature")
    print(f"  [{'PASS' if ok_t else 'FAIL'}] tamper: without the not-moved term the rise disappears")
    print()
    print("  CONCLUSION: between the grain and the coherence length the spectral dimension of the projected")
    print("  geometry rises above four, the direction of the full-causal-set result and against CDT, by an")
    print("  amount logarithmic in d_0/eps; above d_0 it is four on a window that grows with |D|; at the")
    print("  resolution d_0 itself it never exceeds 4.5.  Reduction to two appears nowhere.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()
