# =====================================================================
#  berkeDet BENCHMARK — Google Colab edition
#  Four complete Leibniz-determinant evaluations, pure Python, fair.
#
#  HOW TO USE ON COLAB:
#    1. Runtime -> Change runtime type -> CPU  (GPU/TPU give NO benefit:
#       this is pure-Python loop code, no tensor ops)
#    2. Paste this whole cell and run.
#    3. Edit the CONFIG block to push n higher.
#
#  HONEST LIMITS (so you can verify them yourself):
#    - berkeDet materializes ALL n! permutations.
#        n=10:  3.6M perms  ~0.7 GB   fine
#        n=11: 39.9M perms  ~7 GB     fits standard Colab (12.7 GB), tight
#        n=12:  479M perms  ~70+ GB   IMPOSSIBLE to materialize in RAM
#      (the script prints this estimate and refuses politely)
#    - Naive / Narayana with inversion counting: no memory wall,
#      but ~2.5 min each at n=11 and ~35 min each at n=12.
#    - SJT: n=11 ~1 min, n=12 ~15 min. No memory wall (one perm in RAM).
# =====================================================================

import gc, sys, time, random, statistics

# ----------------------------- CONFIG --------------------------------
SEED  = 7
RUNS  = 3          # median of RUNS timings per method per n
N_MIN = 6

MAX_N = {          # raise these yourself; the walls above are real
    "Naive lex + inv-count": 10,   # 11 -> ~2.5 min, 12 -> ~35 min
    "Narayana + inv-count":  10,   # 11 -> ~2.5 min, 12 -> ~35 min
    "SJT + sign flip":       11,   # 12 -> ~15 min
    "berkeDet":              11,   # 12 -> memory-impossible (see math)
}

# Test-your-hypothesis lever: store berkeDet's permutations as
# 'list' (as in the paper) or 'tuple' (smaller, faster to iterate).
BERKEDET_STORAGE = "list"          # try "tuple" and compare!

# ----------------------- shared: 0-based, fair -----------------------
# ALL methods use permutations of {0..n-1} so every method indexes the
# matrix directly with no per-access arithmetic. Identical treatment.

def inv_sign(p):
    """Sign by inversion counting: the classical per-permutation cost."""
    m = 0; n = len(p)
    for i in range(n):
        pi = p[i]
        for j in range(i + 1, n):
            if pi > p[j]:
                m += 1
    return 1 - 2 * (m & 1)

# ------------------------- 1. naive lex + inv ------------------------
def det_naive(M):
    n = len(M); p = list(range(n)); det = 0
    while True:
        s = inv_sign(p); prod = 1
        for i in range(n):
            prod *= M[i][p[i]]
        det += s * prod
        i = n - 2
        while i >= 0 and p[i] >= p[i + 1]:
            i -= 1
        if i < 0:
            return det
        best = -1
        for j in range(i + 1, n):
            if p[j] > p[i] and (best < 0 or p[j] < p[best]):
                best = j
        p[i], p[best] = p[best], p[i]
        p[i + 1:] = sorted(p[i + 1:])          # the explicit SORT
    # unreachable

# ---------------------- 2. Narayana Pandita + inv --------------------
def det_narayana(M):
    n = len(M); p = list(range(n)); det = 0
    while True:
        s = inv_sign(p); prod = 1
        for i in range(n):
            prod *= M[i][p[i]]
        det += s * prod
        i = n - 2
        while i >= 0 and p[i] >= p[i + 1]:
            i -= 1
        if i < 0:
            return det
        j = n - 1
        while p[j] <= p[i]:
            j -= 1
        p[i], p[j] = p[j], p[i]
        p[i + 1:] = p[:i:-1]                   # suffix REVERSAL
    # unreachable

# ------------------------- 3. SJT + O(1) flip ------------------------
def det_sjt(M):
    n = len(M)
    p = list(range(n)); dirs = [-1] * n; sign = 1; det = 0
    while True:
        prod = 1
        for i in range(n):
            prod *= M[i][p[i]]
        det += sign * prod
        mobile = -1; mi = -1
        for i in range(n):
            j = i + dirs[i]
            if 0 <= j < n and p[i] > p[j] and p[i] > mobile:
                mobile = p[i]; mi = i
        if mi == -1:
            return det
        j = mi + dirs[mi]
        p[mi], p[j] = p[j], p[mi]
        dirs[mi], dirs[j] = dirs[j], dirs[mi]
        sign = -sign                            # O(1) sign update
        for i in range(n):
            if p[i] > mobile:
                dirs[i] = -dirs[i]

# ----------------------------- 4. berkeDet ---------------------------
def meryemSign(n):
    b = [1]
    for e in range(1, n):
        c = b[:]
        for d in range(e):
            f = [-x for x in c] if d % 2 == 0 else c[:]
            b.extend(f)
    return b

def meryemPer(n, storage="list"):
    """All n! permutations of {0..n-1}, lexicographic.
    Block replication + index mapping. No sort/reverse/scan anywhere."""
    wrap = tuple if storage == "tuple" else list
    a = [wrap([0])]
    for b in range(2, n + 1):
        c = list(range(b))
        d = []
        for f in c:
            g = [e for e in c if e != f]       # remaining, increasing
            for h in a:
                d.append(wrap([f] + [g[i] for i in h]))
        a = d
    return a

def det_berkeDet(M):
    n = len(M)
    signs = meryemSign(n)
    perms = meryemPer(n, BERKEDET_STORAGE)
    det = 0
    for s, p in zip(signs, perms):
        prod = 1
        for i in range(n):
            prod *= M[i][p[i]]
        det += s * prod
    return det

# ------------------------ memory guard for berkeDet ------------------
def berkedet_mem_estimate_gb(n):
    """Rough RAM to materialize n! permutations (+ signs + prev layer)."""
    import math
    f = math.factorial(n)
    per = (56 + 8 * n) if BERKEDET_STORAGE == "list" else (40 + 8 * n)
    total = f * (per + 8)                      # items + outer-list slots
    total += math.factorial(n - 1) * (per + 8) # previous layer, peak
    total += f * 8                             # signs list slots
    return total / 1e9

def available_ram_gb():
    try:
        import psutil
        return psutil.virtual_memory().available / 1e9
    except Exception:
        return None

# ------------------------------ runner -------------------------------
random.seed(SEED)
methods = [
    ("Naive lex + inv-count", det_naive),
    ("Narayana + inv-count", det_narayana),
    ("SJT + sign flip",      det_sjt),
    ("berkeDet",             det_berkeDet),
]

N_MAX_GLOBAL = max(MAX_N.values())
print(f"berkeDet storage mode: {BERKEDET_STORAGE!r}   (try 'tuple' too!)")
print(f"{'n':>3} | " + " | ".join(f"{name:>22}" for name, _ in methods))
print("-" * (6 + 25 * len(methods)))

results = {name: {} for name, _ in methods}
for n in range(N_MIN, N_MAX_GLOBAL + 1):
    M = [[random.randint(-9, 9) for _ in range(n)] for _ in range(n)]
    ref = None
    row = []
    for name, fn in methods:
        if n > MAX_N[name]:
            row.append(None); continue
        if name == "berkeDet":
            need = berkedet_mem_estimate_gb(n)
            avail = available_ram_gb()
            if avail is not None and need > 0.75 * avail:
                print(f"  [berkeDet n={n}: needs ~{need:.1f} GB, "
                      f"available ~{avail:.1f} GB -> SKIPPED. "
                      f"This is the materialization wall.]")
                row.append(None); continue
        times = []
        for _ in range(RUNS):
            gc.collect()
            t0 = time.perf_counter()
            d = fn(M)
            t1 = time.perf_counter()
            times.append(t1 - t0)
        if ref is None:
            ref = d
        assert d == ref, f"MISMATCH {name} n={n}: {d} != {ref}"
        med = statistics.median(times)
        results[name][n] = med
        row.append(med)
    print(f"{n:>3} | " + " | ".join(
        (f"{t:>20.4f}s" if t is not None else f"{'-':>21}") for t in row))

print("\nAll executed methods agree on every determinant (asserted).")

# ------------------------------ speedups ------------------------------
for n in sorted(results["berkeDet"]):
    line = [f"n={n}: berkeDet speedup"]
    bt = results["berkeDet"][n]
    for name in ("Naive lex + inv-count", "Narayana + inv-count",
                 "SJT + sign flip"):
        if n in results[name]:
            line.append(f"vs {name.split(' +')[0]}: "
                        f"{results[name][n] / bt:.2f}x")
    print("   ".join(line))

# ------------------------------- plot ---------------------------------
try:
    import matplotlib.pyplot as plt
    plt.figure(figsize=(8, 5))
    for name, _ in methods:
        xs = sorted(results[name]); ys = [results[name][x] for x in xs]
        if xs:
            plt.plot(xs, ys, marker="o", label=name)
    plt.yscale("log"); plt.xlabel("n"); plt.ylabel("seconds (log)")
    plt.title("Leibniz determinant: full-enumeration methods (pure Python)")
    plt.grid(True, which="both", alpha=0.3); plt.legend()
    plt.show()
except Exception as e:
    print(f"(plot skipped: {e})")