#!/usr/bin/env python3
"""
PIST Fiedler-Aware Chiral Boundary Detection

Extends PIST spectral analysis (SpectralN.lean) with Fiedler vector
sign-pattern analysis for chiral boundary classification.

References:
    - formal/SilverSight/PIST/SpectralN.lean (shift-deflation, Fiedler)
    - formal/SilverSight/PIST/CartanConnection.lean (D=1792, crossing weights)
    - formal/CoreFormalism/BraidStateN.lean (chiral state enum)
"""

import numpy as np
from typing import Tuple, Optional

# ── PIST constants (from I₂) ──────────────────────────────────────────

SIGMA_Q16 = 9984   # 39/256 in Q16_16 units
TAU_Q16   = 9362   # 1/7 ≈ 9362/65536
D         = 1792   # lcm(7, 256)
SCALE     = 65536

CHIRAL_LABELS = ["achiral_stable", "left_handed", "right_handed", "chiral_scarred"]


def build_laplacian_8x8(cross_coupling: float = 1e-6) -> np.ndarray:
    """Build graph Laplacian from the Sidon crossing matrix.
    
    The crossing matrix C has:
      C[i,i] = σ = 39/256      (self-weight, diagonal)
      C[i,j] = τ = 1/7         (paired strands, same block)
      C[i,j] = ε               (cross-block, small coupling)
    
    The small cross-block coupling ε breaks the 4-block degeneracy
    so the Fiedler vector is well-defined.
    
    Args:
        cross_coupling: tiny cross-block weight to regularize (default 1e-6)
    
    Adjacency A = off-diagonal entries.
    Degree D[i,i] = sum_j A[i,j].
    Laplacian L = D - A.
    """
    C = np.zeros((8, 8), dtype=np.float64)
    for i in range(8):
        C[i, i] = 39 / 256
        for j in range(8):
            if i != j:
                if i // 2 == j // 2:
                    C[i, j] = 1 / 7
                else:
                    C[i, j] = cross_coupling

    A = C.copy()
    np.fill_diagonal(A, 0.0)
    D = np.diag(A.sum(axis=1))
    L = D - A
    return L


def power_iteration(mat: np.ndarray, max_iter: int = 100,
                    tol: float = 1e-8) -> Tuple[float, np.ndarray]:
    """Dominant eigenvalue and eigenvector via power iteration.
    
    Args:
        mat: n×n symmetric matrix
        max_iter: maximum iterations
        tol: convergence tolerance (residual)
    
    Returns:
        (eigenvalue, eigenvector)
    """
    n = mat.shape[0]
    v = np.arange(1.0, n + 1.0)
    
    for _ in range(max_iter):
        mv = mat @ v
        eig = np.dot(v, mv) / np.dot(v, v)
        norm = np.linalg.norm(mv)
        if norm < 1e-15:
            break
        v_new = mv / norm
        resid = np.linalg.norm(mv - eig * v) / n
        v = v_new
        if resid < tol:
            break
    
    mv = mat @ v
    eig = np.dot(v, mv) / np.dot(v, v)
    return eig, v


def fiedler_vector(L: np.ndarray) -> Tuple[float, np.ndarray]:
    """Compute Fiedler value (2nd smallest eigenvalue) and vector.
    
    Uses full eigendecomposition. For n=8 this is trivially small.
    For larger n, use shift-deflation power iteration (SpectralN.lean).
    
    Args:
        L: n×n Laplacian matrix
    
    Returns:
        (fiedler_value, fiedler_vector)
    """
    eig_vals, eig_vecs = np.linalg.eigh(L)
    # Fiedler = second smallest eigenvalue
    fiedler_val = eig_vals[1]
    fiedler_vec = eig_vecs[:, 1]
    return fiedler_val, fiedler_vec


def classify_chiral_boundary(fiedler_vec: np.ndarray) -> str:
    """Classify chiral boundary state from Fiedler vector sign pattern.
    
    The Fiedler vector has one component per strand (8 total, 4 pairs).
    Sign pattern across paired strands determines chirality:
    
    Pattern                    | Chirality
    ---------------------------------------------------------
    All + (or all -)           | achiral_stable (no boundary)
    Mixed per pair (+, -)      | left_handed (mass bias)
    Mixed per pair (-, +)      | right_handed (vector bias)
    Both pairs strongly mixed  | chiral_scarred (topological defect)
    
    Args:
        fiedler_vec: 8-component Fiedler eigenvector
    
    Returns:
        chiral label string
    """
    sign = np.sign(fiedler_vec)
    
    # Count sign flips within each pair
    intra_flips = 0
    for k in range(4):
        if sign[2*k] != sign[2*k+1]:
            intra_flips += 1
    
    # Count sign flips between adjacent pairs
    inter_flips = 0
    for k in range(3):
        if sign[2*k+1] != sign[2*k+2]:
            inter_flips += 1
    
    # Compute pair-wise net sign: bias within each pair
    bias = []
    for k in range(4):
        pair_sign = fiedler_vec[2*k] + fiedler_vec[2*k+1]
        bias.append(pair_sign)
    
    net_bias = sum(bias)
    
    # Classification rules
    if intra_flips == 0 and inter_flips == 0:
        return "achiral_stable"      # all same sign
    elif intra_flips > 0 and net_bias < 0:
        return "left_handed"         # mass bias (negative)
    elif intra_flips > 0 and net_bias > 0:
        return "right_handed"        # vector bias (positive)
    else:
        return "chiral_scarred"      # mixed topological defect


def compute_chiral_boundary_profile(C_matrix: Optional[np.ndarray] = None) -> dict:
    """Full chiral boundary analysis of the Sidon crossing matrix.
    
    Returns:
        dict with keys: fiedler_value, fiedler_vector, chiral_label,
                        intra_pair_flips, inter_pair_flips, spectral_gap
    """
    if C_matrix is not None:
        L = build_laplacian_from_matrix(C_matrix)
    else:
        L = build_laplacian_8x8()
    
    # Fiedler analysis
    f_val, f_vec = fiedler_vector(L)
    chiral_label = classify_chiral_boundary(f_vec)
    
    # Spectral gap
    lambda_max, _ = power_iteration(L)
    spectral_gap = lambda_max - f_val
    
    return {
        "fiedler_value": float(f_val),
        "fiedler_vector": f_vec.tolist(),
        "chiral_label": chiral_label,
        "intra_pair_flips": sum(1 for k in range(4) if np.sign(f_vec[2*k]) != np.sign(f_vec[2*k+1])),
        "inter_pair_flips": sum(1 for k in range(3) if np.sign(f_vec[2*k+1]) != np.sign(f_vec[2*k+2])),
        "spectral_gap": float(spectral_gap),
        "dominant_eigenvalue": float(lambda_max),
    }


def build_laplacian_from_matrix(mat: np.ndarray) -> np.ndarray:
    """Build Laplacian from arbitrary 8x8 matrix.
    
    Args:
        mat: 8×8 adjacency/intensity matrix (Int or float)
    
    Returns:
        8×8 Laplacian
    """
    A = np.abs(mat).astype(np.float64)
    np.fill_diagonal(A, 0.0)
    D = np.diag(A.sum(axis=1))
    return D - A


# ── Demo ──────────────────────────────────────────────────────────────

def demo():
    print("=" * 60)
    print("PIST Fiedler-Aware Chiral Boundary Detection")
    print("=" * 60)
    
    L = build_laplacian_8x8()
    print(f"\nLaplacian L:\n{L}")
    
    lambda_max, v1 = power_iteration(L)
    print(f"\nλ_max (dominant): {lambda_max:.6f}")
    
    f_val, f_vec = fiedler_vector(L)
    print(f"Fiedler value (λ₂): {f_val:.6f}")
    print(f"Spectral gap: {lambda_max - f_val:.6f}")
    print(f"Fiedler vector: {np.array2string(f_vec, precision=6, suppress_small=True)}")
    print(f"Sign pattern: {np.array2string(np.sign(f_vec), precision=0, suppress_small=True)}")
    
    profile = compute_chiral_boundary_profile()
    print(f"\nChiral classification: {profile['chiral_label']}")
    print(f"Intra-pair sign flips: {profile['intra_pair_flips']}")
    print(f"Inter-pair sign flips: {profile['inter_pair_flips']}")
    
    # Test on perturbed matrices
    print("\n--- Perturbation analysis ---")
    
    # Left-handed perturbation: add negative bias to pair (0,1)
    L_pert = L.copy()
    L_pert[0, 0] += 0.5  # increase degree for strand 0
    fv, _ = fiedler_vector(L_pert)
    print(f"Left-bias perturbation: Fiedler={fv:.6f}, chiral={classify_chiral_boundary(_)}")


if __name__ == "__main__":
    demo()