#!/usr/bin/env python3
"""
spectral_profile.py — 8D Spectral Profile from Byte-Level Co-occurrence Statistics

⚠️  DISCONNECTED FROM LEAN PIPELINE — This module computes a byte-level
spectral profile that is NOT used by the formal SilverSight system.
The Lean PIST pipeline uses 8×8 braid adjacency matrices from
rrc_pist_predictions_250_v1.json, not byte co-occurrence profiles.

Retained for experimental/Python-side use only.
"""

import math
from typing import List


def compute_spectral_profile(equation: str) -> List[float]:
    """Compute 8D spectral profile from equation structure.
    
    Uses byte-level co-occurrence statistics (not semantic parsing).
    Produces a profile that the chaos game converges to.
    
    The 8 dimensions:
    0 — normalized_byte_entropy: Shannon entropy of byte distribution
    1 — diagonal_strength: Self-transition ratio
    2 — spectral_gap: Dominant / subdominant singular value ratio  
    3 — length_complexity: Length-sensitive complexity score
    4 — asymmetry: Non-symmetry of co-occurrence matrix
    5 — symbol_diversity: Unique byte ratio
    6 — run_structure: Mean run length of repeated bytes
    7 — edge_activity: Activity at byte boundaries (256 wrap)
    """
    if not equation:
        return [0.125] * 8

    data = equation.encode('utf-8')
    n = len(data)
    if n == 0:
        return [0.125] * 8

    # ── Byte frequency ──────────────────────────────────────────────
    freq = [0] * 256
    for b in data:
        freq[b] += 1

    # ── Dimension 0: byte entropy (Shannon, normalized) ────────────
    entropy = 0.0
    for count in freq:
        if count > 0:
            p = count / n
            entropy -= p * math.log2(p)
    max_entropy = math.log2(min(n, 256))
    dim0 = entropy / max_entropy if max_entropy > 0 else 0

    # ── Co-occurrence matrix ────────────────────────────────────────
    C = [[0] * 256 for _ in range(256)]
    for i in range(n):
        C[data[i]][data[(i + 1) % n]] += 1

    total_trans = n

    # ── Dimension 1: diagonal strength ─────────────────────────────
    diag_sum = sum(C[i][i] for i in range(256))
    dim1 = diag_sum / total_trans if total_trans > 0 else 0

    # ── Dimension 2: spectral gap via power iteration ──────────────
    # Build a small representative matrix: collapse 256→8 by byte class
    # Classify bytes into 8 buckets and compute 8x8 transition matrix
    buckets = [0] * 256
    for i in range(256):
        if i < 32:       buckets[i] = 0   # control
        elif i < 48:     buckets[i] = 1   # punctuation/special
        elif i < 58:     buckets[i] = 2   # digits
        elif i < 65:     buckets[i] = 3   # more punctuation
        elif i < 91:     buckets[i] = 4   # uppercase
        elif i < 97:     buckets[i] = 5   # more punctuation
        elif i < 123:    buckets[i] = 6   # lowercase
        elif i < 128:    buckets[i] = 7   # extended ascii
        else:            buckets[i] = i % 8  # unicode spread

    M8 = [[0.0] * 8 for _ in range(8)]
    for i in range(256):
        for j in range(256):
            if C[i][j] > 0:
                bi, bj = buckets[i], buckets[j]
                M8[bi][bj] += C[i][j]
    
    # Normalize
    for i in range(8):
        row_sum = sum(M8[i])
        if row_sum > 0:
            for j in range(8):
                M8[i][j] /= row_sum

    # Power iteration for top 2 eigenvalues
    def power_iter(M, iters=30):
        n = len(M)
        v = [math.sin(i * 1.324717957) + 0.01 for i in range(n)]
        # Normalize
        norm = math.sqrt(sum(x*x for x in v))
        v = [x/norm for x in v]
        
        for _ in range(iters):
            new_v = [0.0] * n
            for i in range(n):
                s = 0.0
                for j in range(n):
                    s += M[i][j] * v[j]
                new_v[i] = s
            norm = math.sqrt(sum(x*x for x in new_v))
            if norm < 1e-15:
                break
            v = [x/norm for x in new_v]
        
        # Rayleigh quotient
        Av = [sum(M[i][j] * v[j] for j in range(n)) for i in range(n)]
        val = sum(v[i] * Av[i] for i in range(n))
        return val, v

    val1, v1 = power_iter(M8)
    
    # Deflate for second eigenvalue
    for i in range(8):
        for j in range(8):
            M8[i][j] -= val1 * v1[i] * v1[j]
    
    val2, _ = power_iter(M8)
    val2 = abs(val2)
    
    if val2 > 1e-10:
        gap = val1 / val2
        dim2 = min(1.0, (gap - 1.0) / 10.0)  # normalize: gap of 11 → 1.0
    else:
        dim2 = 1.0

    # ── Dimension 3: length complexity ──────────────────────────────
    # Short equations → low, long equations → high, with diminishing returns
    dim3 = min(1.0, n / 50.0)

    # ── Dimension 4: asymmetry ──────────────────────────────────────
    # Frobenius norm of (C - C^T)
    diff_sq = 0.0
    total_sq = 0.0
    for i in range(256):
        for j in range(256):
            d = C[i][j] - C[j][i]
            s = C[i][j] + C[j][i]
            diff_sq += d * d
            total_sq += s * s
    if total_sq > 0:
        dim4 = min(1.0, math.sqrt(diff_sq) / math.sqrt(total_sq))
    else:
        dim4 = 0.0

    # ── Dimension 5: symbol diversity ──────────────────────────────
    unique_bytes = sum(1 for f in freq if f > 0)
    dim5 = unique_bytes / min(n, 256) if n > 0 else 0

    # ── Dimension 6: run structure ─────────────────────────────────
    # Mean run length of identical consecutive bytes
    if n > 0:
        runs = []
        current_run = 1
        for i in range(1, n):
            if data[i] == data[i-1]:
                current_run += 1
            else:
                runs.append(current_run)
                current_run = 1
        runs.append(current_run)
        mean_run = sum(runs) / len(runs)
        dim6 = min(1.0, mean_run / (n ** 0.5)) if n > 0 else 0
    else:
        dim6 = 0.0

    # ── Dimension 7: edge_activity ─────────────────────────────────
    # Transitions that cross byte-class boundaries
    cross_count = 0
    for i in range(n):
        if buckets[data[i]] != buckets[data[(i+1) % n]]:
            cross_count += 1
    dim7 = cross_count / n if n > 0 else 0

    profile = [dim0, dim1, dim2, dim3, dim4, dim5, dim6, dim7]
    
    # Normalize
    s = sum(profile)
    if s > 0:
        profile = [p / s for p in profile]
    
    return profile


def profile_to_basin_hint(profile: List[float]) -> str:
    """Get a rough basin hint from the spectral profile."""
    d0, d1, d2, d3, d4, d5, d6, d7 = profile

    void_score = d0 + d6       # entropy + runs (trivial patterns)
    orbit_score = d2 + d7      # spectral gap + edge activity
    braid_score = d1 + d4      # diagonal + asymmetry
    observer_score = d3 + d5   # length + diversity

    scores = {
        "q_void": void_score,
        "q_orbit": orbit_score,
        "q_braid": braid_score,
        "q_observer": observer_score,
    }
    return max(scores, key=scores.get)
