CCFU Script 4 — verify_gl7_orbit.py

"""
verify_gl7_orbit.py — GL(7) orbit certificates
=======================================================

Numerical verifier for the four C2-compatible split patterns.
Compared with the original uploaded script, this version uses scipy.optimize
least_squares on the 35 independent 3-form coefficients, which avoids the long
stall observed in pair 2 with L-BFGS-B.

Run: python3 verify_gl7_orbit.py
"""

from itertools import permutations, combinations
import time

import numpy as np
from scipy.optimize import least_squares

canonical = {
    (0, 1, 2): +1, (0, 3, 4): -1, (0, 5, 6): -1,
    (1, 3, 5): -1, (2, 4, 5): +1, (1, 4, 6): +1, (2, 3, 6): +1,
}

BASIS_3 = [(i, j, k) for i in range(7) for j in range(i + 1, 7) for k in range(j + 1, 7)]
GL_POSITIONS = [(a, b) for a in range(7) for b in range(7)]


def permutation_sign(seq):
    sign = 1
    lst = list(seq)
    for a in range(len(lst)):
        for b in range(a + 1, len(lst)):
            if lst[a] > lst[b]:
                sign *= -1
    return sign


def build_omega_tensor(pattern_triples_signed):
    Omega = np.zeros((7, 7, 7))
    for triple, s in pattern_triples_signed.items():
        for p in permutations(triple):
            Omega[p[0], p[1], p[2]] = s * permutation_sign(p)
    return Omega


def assign_signs_to_pattern(triple_list):
    """Return the first sign assignment with full orbit-rank 35."""
    triple_list = sorted(triple_list)
    n = len(triple_list)
    for signs_int in range(2 ** n):
        signs = [(1 if (signs_int >> k) & 1 else -1) for k in range(n)]
        signed = {triple_list[k]: signs[k] for k in range(n)}
        Omega = build_omega_tensor(signed)
        M = np.zeros((35, 49))
        for ri, (i, j, k) in enumerate(BASIS_3):
            for ci, (a, b) in enumerate(GL_POSITIONS):
                val = 0.0
                if b == i:
                    val += Omega[a, j, k]
                if b == j:
                    val += Omega[i, a, k]
                if b == k:
                    val += Omega[i, j, a]
                M[ri, ci] = val
        if np.linalg.matrix_rank(M, tol=1e-9) == 35:
            return signed
    return None


test_patterns_lists = [
    [(0, 1, 2), (0, 1, 4), (0, 2, 5), (0, 5, 6), (1, 4, 6), (2, 3, 5), (3, 4, 6)],
    [(0, 1, 2), (0, 1, 4), (0, 2, 5), (0, 5, 6), (1, 4, 6), (2, 3, 6), (3, 4, 5)],
    [(0, 1, 2), (0, 1, 4), (0, 2, 6), (0, 5, 6), (1, 3, 5), (2, 4, 5), (3, 4, 6)],
    [(0, 1, 2), (0, 1, 4), (0, 3, 5), (0, 5, 6), (1, 2, 6), (2, 4, 5), (3, 4, 6)],
]


def apply_h(h, O):
    return np.einsum('ai,bj,ck,abc->ijk', h, h, h, O)


def normalize_h(h):
    d = np.linalg.det(h)
    if abs(d) < 1e-12:
        return None
    return h * (abs(d) ** (-1.0 / 7.0))


def best_c_and_residual(h_normed, O1, O2):
    applied = apply_h(h_normed, O1)
    c = float(np.sum(applied * O2)) / float(np.sum(O2 * O2))
    diff = applied - c * O2
    return c, float(np.sum(diff ** 2))


def tensor_basis_vector(O):
    return np.array([O[i, j, k] for i, j, k in BASIS_3], dtype=float)


def residual_vector(params, O1, O2_basis):
    h = params[:49].reshape(7, 7)
    c = params[49]
    h_normed = normalize_h(h)
    if h_normed is None:
        return np.ones(len(BASIS_3), dtype=float) * 1e4
    applied = apply_h(h_normed, O1)
    applied_basis = tensor_basis_vector(applied)
    return applied_basis - c * O2_basis


def find_gl7_map(O1, O2, n_attempts=60, max_nfev=2000):
    """Find h and scalar c with h*O1 ~ c*O2.

    least_squares is much more reliable here than finite-difference L-BFGS-B
    on the squared loss; the original method stalls on the second test pair.
    """
    norm_O2 = float(np.sum(O2 ** 2))
    O2_basis = tensor_basis_vector(O2)
    best = (float('inf'), None, None)

    for attempt in range(n_attempts):
        rng = np.random.default_rng(attempt)
        if attempt < 15:
            h_init = np.eye(7) + 0.3 * rng.standard_normal((7, 7))
        elif attempt < 30:
            h_init = np.eye(7) + 0.8 * rng.standard_normal((7, 7))
        elif attempt < 50:
            h_init = rng.standard_normal((7, 7))
        else:
            perm = rng.permutation(7)
            h_init = np.eye(7)[perm] + 0.3 * rng.standard_normal((7, 7))

        h_normed_init = normalize_h(h_init)
        if h_normed_init is None:
            continue

        c_init, _ = best_c_and_residual(h_normed_init, O1, O2)
        x0 = np.concatenate([h_normed_init.flatten(), [c_init]])

        try:
            result = least_squares(
                residual_vector,
                x0,
                args=(O1, O2_basis),
                method='trf',
                max_nfev=max_nfev,
                ftol=1e-12,
                xtol=1e-12,
                gtol=1e-12,
            )
            h_final = result.x[:49].reshape(7, 7)
            h_normed = normalize_h(h_final)
            if h_normed is None:
                continue
            c, res = best_c_and_residual(h_normed, O1, O2)
            rel = res / norm_O2
            if rel < best[0]:
                best = (rel, h_normed, c)
            if rel < 1e-14:
                break
        except Exception as exc:
            print(f"    attempt {attempt}: optimizer exception: {exc}")
            continue

    return best


def hitchin_lambda(O):
    """Compute Hitchin bilinear form b_Omega and determinant."""
    idx = list(range(7))
    b = np.zeros((7, 7))
    for a in range(7):
        for bb in range(7):
            total = 0.0
            for c2a in combinations(idx, 2):
                rem = [x for x in idx if x not in c2a]
                for c2b in combinations(rem, 2):
                    c3 = tuple(x for x in rem if x not in c2b)
                    if len(c3) != 3:
                        continue
                    va = O[a, c2a[0], c2a[1]]
                    vb = O[bb, c2b[0], c2b[1]]
                    vo = O[c3[0], c3[1], c3[2]]
                    if va == 0 or vb == 0 or vo == 0:
                        continue
                    perm = list(c2a) + list(c2b) + list(c3)
                    if len(set(perm)) != 7:
                        continue
                    total += permutation_sign(perm) * va * vb * vo
            b[a][bb] = total / 6.0
    eigs = np.linalg.eigvalsh(b)
    pos = int(np.sum(eigs > 0.01))
    neg = int(np.sum(eigs < -0.01))
    det_b = np.linalg.det(b)
    return (pos, neg), det_b


def main():
    O1 = build_omega_tensor(canonical)
    print(f"Canonical Ω built. ||Ω₁||² = {np.sum(O1 ** 2)}")

    hsig_can, _ = hitchin_lambda(O1)
    orbit_can = 'SPLIT' if hsig_can in [(3, 4), (4, 3)] else 'COMPACT'
    print(f"Canonical: sig(b_Ω)={hsig_can}, orbit={orbit_can}")
    print()

    results = []
    for idx, pat_list in enumerate(test_patterns_lists):
        print(f"\nPair {idx + 1}: Canonical ↔ {pat_list}")
        pattern_signed = assign_signs_to_pattern(pat_list)
        if pattern_signed is None:
            print("  No valid sign assignment. Skipping.")
            continue
        O2 = build_omega_tensor(pattern_signed)
        start = time.time()
        rel_residual, h, c = find_gl7_map(O1, O2, n_attempts=60)
        elapsed = time.time() - start
        det_h = np.linalg.det(h) if h is not None else 0.0
        hsig2, _ = hitchin_lambda(O2)
        orbit = 'SPLIT' if hsig2 in [(3, 4), (4, 3)] else 'COMPACT'
        status = '✓' if rel_residual < 1e-8 else '✗'
        print(f"  {status} residual={rel_residual:.2e}, c={c:.3f}, det(h)={det_h:.3f}, time={elapsed:.1f}s")
        print(f"    sig(b_Ω)={hsig2} → {orbit}")
        results.append((idx + 1, rel_residual < 1e-8, rel_residual, c))

    print(f"\nSUMMARY: {sum(1 for _, s, _, _ in results if s)}/{len(results)} verified.")
    return 0 if all(s for _, s, _, _ in results) else 1


if __name__ == '__main__':
    raise SystemExit(main())