#!/usr/bin/env python3
"""
chaos_game.py — Deterministic Chaos Game Engine

Deterministic chaos game using IFS contraction on 8×8 state matrix.
4 basins: q_void (rows 0-1), q_orbit (rows 2-3), q_braid (rows 4-5), q_observer (rows 6-7).
"""

import hashlib
import json
import math
from typing import Dict, List, Optional

from sidon_address import (
    SIDON_ADDRESSES,
    _ADDRESS_TO_STRAND,
    address_to_strand,
    compute_full_address,
    structural_hash,
    verify_sidon_property,
)
from spectral_profile import compute_spectral_profile

EPSILON = 1e-14
N = 8
IFS_ALPHA = 0.75  # Strong contraction for fast convergence

BASIN_ROWS = {
    "q_void": (0, 1),
    "q_orbit": (2, 3),
    "q_braid": (4, 5),
    "q_observer": (6, 7),
}

LCG_A = 1664525
LCG_C = 1013904223
LCG_M = 2**32

DEFAULT_CONVERGENCE_THRESHOLD = 0.95
DEFAULT_MAX_STEPS = 10000
DEFAULT_CONVERGENCE_WINDOW = 10


class LCG:
    def __init__(self, seed: int):
        self.state = seed & 0xFFFFFFFF
    def next(self) -> int:
        self.state = (LCG_A * self.state + LCG_C) % LCG_M
        return self.state
    def next_float(self) -> float:
        return self.next() / LCG_M


def init_state_matrix(seed: int) -> List[List[float]]:
    """Initialize 8×8 state matrix deterministically from seed."""
    lcg = LCG(seed)
    A = [[0.0] * N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            if i == j:
                A[i][j] = SIDON_ADDRESSES[i] / 128.0 + 0.5
            elif abs(i - j) == 1:
                A[i][j] = -0.1 + 0.04 * (lcg.next_float() - 0.5)
            else:
                A[i][j] = 0.05 * (lcg.next_float() - 0.5)
    return A


def mat_copy(A): return [row[:] for row in A]
def mat_diff_norm(A, B):
    return math.sqrt(sum((A[i][j] - B[i][j])**2 for i in range(N) for j in range(N)))
def mat_norm(A):
    return math.sqrt(sum(A[i][j]**2 for i in range(N) for j in range(N)))


def strand_to_basin(strand: int) -> str:
    if strand < 2:   return "q_void"
    elif strand < 4: return "q_orbit"
    elif strand < 6: return "q_braid"
    else:            return "q_observer"


def get_basin_rows(strand: int):
    return BASIN_ROWS[strand_to_basin(strand)]


def ifs_contract(A, strand, step, eq_hash):
    """Apply IFS contraction toward a strand's quadrant.
    
    Pure IFS contraction: A <- (1-alpha)*A + alpha*T where T is the
    target matrix with strong energy in the target strand's basin and
    suppressed energy elsewhere. No post-step modifications.
    """
    alpha = IFS_ALPHA
    r0, r1 = get_basin_rows(strand)
    lcg = LCG((eq_hash + step * 104729 + strand * 7919) & 0xFFFFFFFF)
    
    for i in range(N):
        for j in range(N):
            in_basin = (r0 <= i <= r1)
            on_diag = (i == j)
            is_strand = (i == strand)
            
            if is_strand and on_diag:
                target = 10.0  # maximum energy at strand diagonal
            elif is_strand:
                target = 3.0 + 0.5 * lcg.next_float()
            elif in_basin and on_diag:
                target = 4.0 + 0.5 * lcg.next_float()
            elif in_basin:
                target = 1.5 + 0.3 * lcg.next_float()
            elif on_diag:
                target = 0.02 + 0.01 * lcg.next_float()
            else:
                target = 0.005 * lcg.next_float()
            
            A[i][j] = (1 - alpha) * A[i][j] + alpha * target


def quadrant_energy(A):
    """Compute Frobenius energy in each basin (2-row block)."""
    energy = {}
    for basin, (r0, r1) in BASIN_ROWS.items():
        e = sum(A[i][j]**2 for i in range(r0, r1 + 1) for j in range(N))
        energy[basin] = math.sqrt(e)
    energy["total"] = sum(v for k, v in energy.items())
    return energy


def energy_ratio(A):
    """Ratio of dominant basin energy to total energy."""
    qe = quadrant_energy(A)
    total = qe["total"]
    if total < EPSILON:
        return 0.0
    basin_energies = {k: v for k, v in qe.items() if k != "total"}
    return max(basin_energies.values()) / total


def dominant_basin(A):
    qe = quadrant_energy(A)
    del qe["total"]
    return max(qe, key=qe.get)


def detect_quarantine(equation):
    if not equation or not equation.strip():
        return "empty_equation"
    eq = equation.strip()
    normalized = eq.replace(" ", "").replace("\t", "")
    contradictions = {"0=1", "1=0", "false=true", "true=false",
                      "False=True", "True=False", "⊥=⊤", "⊤=⊥"}
    if normalized in contradictions:
        return "explicit_contradiction"
    if eq in {"0 = 1", "1 = 0", "False = True", "True = False", "⊥ = ⊤", "⊤ = ⊥"}:
        return "explicit_contradiction"
    return None


def sidon_guided_chaos_game(
    target_address: list,
    max_steps: int = DEFAULT_MAX_STEPS,
    convergence_threshold: float = DEFAULT_CONVERGENCE_THRESHOLD,
    convergence_window: int = DEFAULT_CONVERGENCE_WINDOW,
    equation: str = "",
):
    """Deterministic chaos game guided by Sidon address."""
    # Quarantine check
    quarantine = detect_quarantine(equation) if equation else None
    if quarantine:
        return {
            "converged": False,
            "basin": "QUARANTINE",
            "steps": -1,
            "energy_ratio": 0.0,
            "address": target_address,
            "hash": hex(structural_hash(equation))[2:18] if equation else "",
            "target_strand": -1,
            "quarantine": quarantine,
        }
    
    primary = target_address[0] if target_address else SIDON_ADDRESSES[0]
    primary_strand = address_to_strand(primary)
    eq_hash = structural_hash(equation) if equation else 42
    A = init_state_matrix(eq_hash & 0xFFFFFFFF)
    
    converged = False
    basin_history = []
    trajectory = []
    
    for step in range(max_steps):
        # Adaptive: emphasize target strand more over time
        progress = min(step / max(max_steps // 3, 1), 1.0)
        target_prob = 0.5 + 0.45 * progress
        
        lcg = LCG((eq_hash + step * 104729) & 0xFFFFFFFF)
        if lcg.next_float() < target_prob:
            chosen = primary_strand
        else:
            others = [s for s in range(N) if s != primary_strand]
            chosen = others[step % len(others)]
        
        ifs_contract(A, chosen, step, eq_hash)
        trajectory.append(chosen)
        
        # Check convergence every 4 steps
        if step % 4 == 0 and step > 0:
            ratio = energy_ratio(A)
            current_basin = dominant_basin(A)
            basin_history.append(current_basin)
            
            if ratio >= convergence_threshold:
                if len(basin_history) >= convergence_window:
                    recent = basin_history[-convergence_window:]
                    if len(set(recent)) == 1:
                        converged = True
                        break
    
    steps = step + 1 if converged else max_steps
    final_ratio = energy_ratio(A)
    final_basin = dominant_basin(A)
    
    if not converged and final_ratio >= convergence_threshold:
        converged = True
    
    qe = quadrant_energy(A)
    profile = energy_to_profile(qe)
    from sidon_address import spectral_to_sidon_address
    achieved = spectral_to_sidon_address(profile, eq_hash)
    
    return {
        "converged": converged,
        "basin": final_basin,
        "steps": steps,
        "energy_ratio": round(final_ratio, 6),
        "address": achieved,
        "hash": hex(eq_hash)[2:18],
        "target_strand": primary_strand,
        "quarantine": None,
        "trajectory": trajectory[:100],
    }


def energy_to_profile(energy):
    """Convert quadrant energies to 8D profile."""
    total = energy.get("total", 1.0)
    if total < EPSILON:
        total = 1.0
    v = energy.get("q_void", 0.0) / total
    o = energy.get("q_orbit", 0.0) / total
    b = energy.get("q_braid", 0.0) / total
    ob = energy.get("q_observer", 0.0) / total
    profile = [v, v*v, o, o*o, b, b*b, ob, ob*ob]
    s = sum(profile)
    if s > 0:
        profile = [p / s for p in profile]
    else:
        profile = [0.125] * 8
    return profile


def generate_receipt(results, schema_version="stage3_v1"):
    import datetime
    converged = sum(1 for r in results if r.get("converged"))
    quarantined = sum(1 for r in results if r.get("quarantine"))
    basin_counts = {"q_void": 0, "q_orbit": 0, "q_braid": 0, "q_observer": 0}
    for r in results:
        b = r.get("basin", "")
        if b in basin_counts:
            basin_counts[b] += 1
    receipt = {
        "schema": f"rrc_chaos_game_search_{schema_version}",
        "sidon_property_verified": verify_sidon_property(),
        "total_searches": len(results),
        "converged": converged,
        "quarantined": quarantined,
        "failed": len(results) - converged - quarantined,
        "basin_distribution": basin_counts,
        "convergence_rate": round(converged / len(results), 4) if results else 0.0,
        "parameters": {
            "matrix_size": N,
            "ifs_alpha": IFS_ALPHA,
            "convergence_threshold": DEFAULT_CONVERGENCE_THRESHOLD,
            "max_steps": DEFAULT_MAX_STEPS,
            "convergence_window": DEFAULT_CONVERGENCE_WINDOW,
        },
        "results": results,
        "computed_at": datetime.datetime.now(datetime.timezone.utc).isoformat(),
    }
    canonical = json.dumps(receipt, sort_keys=True, separators=(",", ":"))
    receipt["receipt_sha256"] = hashlib.sha256(canonical.encode()).hexdigest()
    return receipt


def search_equation(equation: str, **kwargs):
    """Full pipeline: equation → profile → address → chaos game → result."""
    profile = compute_spectral_profile(equation)
    address = compute_full_address(equation)
    result = sidon_guided_chaos_game(
        target_address=address,
        equation=equation,
        **kwargs,
    )
    result["spectral_profile"] = [round(x, 6) for x in profile]
    return result