#!/usr/bin/env python3
"""
walsh_dirac_lipnorm.py
=======================
Tests whether the Dirac-commutator Lip-norm

    L_Dirac(T) := ||[D, T]||_op,    D = sum_a chi_a  (Boolean-CAR Dirac)

is bounded above and below by an EXPONENTIAL Walsh-weight Lip-norm
e^{c|S|}.  If yes, the finding of walsh_round_trip2.py
(L_round needs exponential Lip-norm) is automatic in PST's
spectral-triple framework, without any artificial weight choice.

Tests for T = chi_S at varying weights k = |S|, at D in {4, 6, 8, 10, 12}.

Construction:
  D     = sum_a chi_a  acting on H_D = C^{2^D}  (Jordan-Wigner)
  [D, chi_S] = sum_a [chi_a, chi_S]
  [chi_a, chi_S] = 0 if a not in S
  [chi_a, chi_S] = 2 chi_a chi_S if a in S  (anticommute by 1 Pauli sign change)

So ||[D, chi_S]||_op <= 2|S|  trivially.

Below we measure ||[D, chi_S]||_op as a function of |S|.  The question
is whether it grows EXPONENTIALLY (matching e^{c|S|} with c > log 2)
or LINEARLY (matching the naive bound).
"""
import math
import numpy as np
import numpy.linalg as la
from itertools import combinations

sz = np.array([[1, 0], [0, -1]], dtype=complex)
sx = np.array([[0, 1], [1, 0]], dtype=complex)
sy = np.array([[0, -1j], [1j, 0]], dtype=complex)
I2 = np.eye(2, dtype=complex)
sp = 0.5 * (sx - 1j * sy)


def kron_chain(ops):
    out = ops[0]
    for op in ops[1:]:
        out = np.kron(out, op)
    return out


def clifford_chi_a(D, a):
    """chi_a = c_a + c_a^dag via JW:  sigma_z^{<a} sigma_x^{(a)} I^{>a}."""
    chain = [sz] * a + [sx] + [I2] * (D - 1 - a)
    return kron_chain(chain)


def walsh_chi_S(D, S):
    """chi_S as 2^D x 2^D diagonal matrix."""
    factors = [sz if a in S else I2 for a in range(D)]
    return kron_chain(factors)


def discrete_dirac(D):
    """D_disc = sum_a chi_a as a 2^D x 2^D matrix."""
    N = 1 << D
    out = np.zeros((N, N), dtype=complex)
    for a in range(D):
        out = out + clifford_chi_a(D, a)
    return out


def commutator(A, B):
    return A @ B - B @ A


def op_norm(M):
    """Spectral norm (largest singular value)."""
    return float(la.norm(M, ord=2))


# ---------- Sweep ----------
print("=" * 80)
print("  walsh_dirac_lipnorm.py  --  Dirac-commutator Lip-norm vs Walsh weight")
print("=" * 80)
print()
print("  L_Dirac(chi_S) := ||[sum_a chi_a, chi_S]||_op")
print()
print(f"  {'D':>3} {'|S|':>4} {'||[D, chi_S]||':>15} {'/ 2|S|':>10} {'/ e^{|S| log 2}':>16} {'/ e^{1.5|S|}':>14}")
print(f"  {'-'*3} {'-'*4} {'-'*15} {'-'*10} {'-'*16} {'-'*14}")
for D in (4, 6, 8):  # D = 10, 12 dropped for Pyodide-friendly runtime (2^12 = 4096 matrices push op_norm cost over the budget)
    D_op = discrete_dirac(D)
    for k in range(1, D + 1):
        # Sample one S of weight k -- use S = {0, 1, ..., k-1}
        S = tuple(range(k))
        chi_S_op = walsh_chi_S(D, S)
        comm = commutator(D_op, chi_S_op)
        n = op_norm(comm)
        # Ratios vs linear (2|S|), exponential e^{log 2 |S|} = 2^|S|, and e^{1.5|S|}
        r_lin = n / (2 * k) if k > 0 else float('nan')
        r_e1  = n / math.exp(math.log(2) * k)
        r_e15 = n / math.exp(1.5 * k)
        print(f"  {D:>3} {k:>4} {n:>15.4f} {r_lin:>10.4f} {r_e1:>16.4f} {r_e15:>14.4f}")
    print()

print("=" * 80)
print("  Interpretation")
print("=" * 80)
print()
print("  If ||[D, chi_S]||_op is linear in |S| (the naive bound 2|S|), then L_Dirac")
print("  is POLYNOMIAL in Walsh weight.  Hence L_Dirac is NOT exponential and does")
print("  NOT match the e^{c|S|} Lip-norm that walsh_round_trip2.py found necessary.")
print()
print("  In that case, the natural single-Dirac-commutator Lip-norm is")
print("  INSUFFICIENT to close L_round.  Either:")
print("    (a) iterated commutators ||[D, [D, ... [D, T] ...]]||  are needed, or")
print("    (b) the Lip-norm is sup over a richer family, e.g.")
print("         L_iter(T) := sup_k ||[D^k, T]||^{1/k},")
print("        which gives the spectral radius (asymptotic growth rate).")
print()
print("  If ||[D, chi_S]||_op grows exponentially with |S|, L_Dirac IS exponential")
print("  and matches the finding directly.")