#!/usr/bin/env python3
"""
center_norm.py -- pairwise overlap vs complete cross-sector coherence (Observation C).

For controlled records with cross Gram matrix G (m x m), the complete cross-sector
coherence is the Schur-multiplier norm of G on trace class (the package's D16 identity
alpha_Z = ||S_G||_cb,S1). This script checks the two families used in the note, WITHOUT
an SDP solver, by sandwiching the norm:
  lower bound  ||G o X||_1 / ||X||_1 for explicit test inputs X (here X = all-ones / m),
  upper bound  max_a ||x_a|| * max_b ||y_b|| for an explicit factorisation G_ab = <x_a, y_b>
               (the Grothendieck-Haagerup characterisation of Schur-multiplier norms).
Families: Fourier records G_ab = (r/sqrt m) exp(2 pi i ab/m)  ->  delta = r/sqrt m, alpha = r
          center-aligned records G = r I                        ->  delta = r,        alpha = r
Also prints the entrywise bound delta <= alpha <= min(1, sqrt(m) delta).
No files are written.
"""
import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ.setdefault(_v, "1")          # small dense blocks: extra BLAS threads only add contention
import numpy as np


def trace_norm(A):
    return float(np.linalg.svd(A, compute_uv=False).sum())


def fourier(m, r):
    a = np.arange(m)
    return (r / np.sqrt(m)) * np.exp(2j * np.pi * np.outer(a, a) / m)


def bounds(G):
    m = G.shape[0]
    X = np.ones((m, m)) / m                           # rank one, ||X||_1 = 1
    lower = trace_norm(G * X) / trace_norm(X)
    # factorisation G = U diag(s) V^H  ->  x_a = sqrt(s) U_a^*, y_b = sqrt(s) V_b^*  (always valid)
    U, s, Vh = np.linalg.svd(G)
    xa = U * np.sqrt(s); yb = Vh.conj().T * np.sqrt(s)
    up_svd = np.linalg.norm(xa, axis=1).max() * np.linalg.norm(yb, axis=1).max()
    # trivial factorisation x_a = G_a (row), y_b = e_b
    up_row = np.linalg.norm(G, axis=1).max()
    return lower, min(up_svd, up_row)


def main():
    print("Observation C: entrywise overlap delta vs cross-sector coherence alpha (sandwiched)")
    print(f"{'family':>10} {'m':>4} {'r':>5} {'delta':>9} {'alpha LB':>9} {'alpha UB':>9} {'sqrt(m) delta':>13}")
    ok = True
    for fam in ("fourier", "aligned"):
        for m in (1, 4, 16, 48):
            for r in (0.5, 0.8):
                G = fourier(m, r) if fam == "fourier" else r * np.eye(m)
                d = np.abs(G).max()
                lb, ub = bounds(G)
                print(f"{fam:>10} {m:>4d} {r:>5.2f} {d:>9.5f} {lb:>9.5f} {ub:>9.5f} {min(1, np.sqrt(m) * d):>13.5f}")
                ok &= abs(lb - r) < 1e-12 and abs(ub - r) < 1e-12 and d - 1e-12 <= lb and ub <= min(1, np.sqrt(m) * d) + 1e-12
    print("both families: alpha = r exactly (lower = upper bound); Fourier saturates alpha = sqrt(m) delta"
          if ok else "CHECK FAILED")


if __name__ == "__main__":
    main()