"""
11D navigation application layer.

Uses the Fractal Correction Engine for error-corrected quantum evolution
within an M-theory-motivated 11-dimensional parameter space.

The 11D structure:
  - Dimensions 0-3: 4D spacetime (Minkowski)
  - Dimensions 4-10: 7 compactified extra dimensions (Calabi-Yau / G2)

Key physics integration:
  - Spectral decomposition: Hamiltonian eigenvalues map to energy bands
    corresponding to each dimension. Per-dimension curvature is computed
    from spectral subspace projections of the density matrix.
  - KK activation: compact dimensions activate when system energy
    approaches their Kaluza-Klein mass threshold (sigmoid gating).
  - Geometric continuity: observed spacetime curvature (dims 0-3)
    extrapolates compact dimension states (dims 4-10) via the FCE's
    curvature transfer function.
  - QFT boundary effects: Casimir energy and Hawking temperature at
    dimensional transitions.
  - Thermodynamic cost: Jarzynski/Crooks free energy accounting.
  - Holographic bounds: Bekenstein-Hawking entropy, MSS chaos bound.

References:
    [1] Horava & Witten, Nucl.Phys.B460:506-524 (1996)
    [2] Alsing & Cafaro, arXiv:2311.18458 (2024)
"""

import logging
import numpy as np
from dataclasses import dataclass, field
from typing import List, Optional, Dict, Any, Tuple

from .config import PHYS
from .engine import FractalCorrectionEngine, EngineConfig, EvolutionResult
from .fidelity import uhlmann_fidelity, validate_density_matrix
from .lindblad import build_lindblad_operators, compute_decoherence_rates
from .frenet_serret import QuantumFrenetSerret
from .holographic import HolographicSystem
from .nonequilibrium import NonEquilibriumDynamics
from .qft_effects import QuantumFieldTheoryEffects

# Optional KK Hamiltonian import (avoids circular dependency)
try:
    from .kk_hamiltonian import KKTowerBuilder
except ImportError:
    KKTowerBuilder = None

logger = logging.getLogger(__name__)


# ======================================================================
# Data structures
# ======================================================================

@dataclass
class DimensionInfo:
    """Physical description of one dimension."""
    index: int
    name: str
    dim_type: str               # 'spacetime' or 'compact'
    compactification_radius: float  # meters (inf for spacetime dims)
    kk_mass_scale: float        # GeV (0 for spacetime dims)
    suppression_factor: float   # Horava-Witten suppression (1.0 for spacetime)


@dataclass
class NavigationState:
    """Current state of the 11D navigation system."""
    coordinates: np.ndarray     # 11D coordinate vector
    momenta: np.ndarray         # 11D momentum vector
    quantum_state: np.ndarray   # Density matrix of the quantum subsystem
    time: float
    dimension_info: List[DimensionInfo]


@dataclass
class DimensionalCurvatureTensor:
    """Per-dimension Frenet-Serret curvature and torsion."""
    curvatures: np.ndarray        # (11,) kappa per dimension
    torsions: np.ndarray          # (11,) tau per dimension
    activation: np.ndarray        # (11,) 0-1 activation level per dim
    subspace_weights: np.ndarray  # (11,) Tr(P_d @ rho) per subspace


@dataclass
class DimensionalTransitionState:
    """Transition state tracking for dimensional activation."""
    active_dimensions: List[int]
    activation_energies: np.ndarray      # (11,) KK threshold per dim
    current_energy: float
    transition_history: List[Tuple[float, int, str]]


@dataclass
class NavigationCost:
    """Thermodynamic cost of navigation via Jarzynski/Crooks."""
    free_energy_change: float
    dissipated_work: float
    entropy_production: float
    landauer_cost: float
    work_values: np.ndarray
    crooks_verification: Optional[Dict[str, float]]


@dataclass
class NavigationConfig:
    """Configuration for the enhanced 11D navigation system."""
    kk_activation_sharpness: float = 10.0
    enable_boundary_effects: bool = True
    navigation_temperature: float = 1.0
    chaos_check_interval: int = 25
    holographic_check_interval: int = 25


@dataclass
class NavigationResult:
    """Result of a navigation trajectory computation."""
    states: List[NavigationState]
    evolution_result: EvolutionResult
    dimension_curvatures: np.ndarray       # (n_steps+1, 11) curvature per dim
    dimension_vulnerabilities: np.ndarray  # Which dims need most correction
    total_fidelity: float
    summary: Dict[str, Any]

    # Full physics integration fields
    curvature_tensors: List[DimensionalCurvatureTensor] = field(
        default_factory=list
    )
    transition_state: Optional[DimensionalTransitionState] = None
    navigation_cost: Optional[NavigationCost] = None
    holographic_metrics: Dict[str, Any] = field(default_factory=dict)
    chaos_metrics: Dict[str, Any] = field(default_factory=dict)
    qft_boundary_effects: Dict[int, Dict[str, float]] = field(
        default_factory=dict
    )


# ======================================================================
# M-theory dimensional scaling (unchanged)
# ======================================================================

class DimensionalScaling:
    """
    M-theory dimensional scaling factors.

    Computes suppression factors for extra dimensions based on:
    - Kaluza-Klein mass tower
    - Horava-Witten warping
    - AdS/CFT holographic weight
    """

    def __init__(self):
        self.planck_mass = PHYS.planck_mass
        self.string_scale = PHYS.string_scale

    def kk_mass(self, dimension: int, mode_number: int) -> float:
        """KK mode mass: m_n = n/R for compactification radius R."""
        if dimension < 4 or mode_number == 0:
            return 0.0
        config_key = dimension + 1
        radius = PHYS.extra_dim_bounds.get(config_key, 1e-21)
        radius_natural = radius * 5.07e15  # meters -> GeV^-1
        return mode_number / radius_natural

    def horava_witten_suppression(self, dimension: int) -> float:
        """Exponential warping suppression from HW domain-wall theory."""
        if dimension < 4:
            return 1.0
        warp_k = self.planck_mass / self.string_scale
        config_key = dimension + 1
        distance = np.sqrt(PHYS.extra_dim_bounds.get(config_key, 1e-21))
        suppression = np.exp(-warp_k * distance * 1e21)
        return max(suppression, 1e-10)

    def ads_cft_weight(self, bulk_dim: int = 11, boundary_dim: int = 4) -> float:
        """Holographic conformal weight."""
        d = boundary_dim
        m2_R2 = (bulk_dim - boundary_dim) ** 2
        delta = (d + np.sqrt(d ** 2 + m2_R2)) / 2
        return delta / (bulk_dim + 1)

    def scaling_factor(self, dimension: int) -> float:
        """Combined scaling factor for a dimension."""
        if dimension < 4:
            return 1.0 + 0.2 * dimension
        hw = self.horava_witten_suppression(dimension)
        ads = self.ads_cft_weight()
        kk_sum = sum(
            1.0 / (1.0 + self.kk_mass(dimension, n))
            for n in range(1, 5)
        )
        return hw * ads * (1.0 + kk_sum)

    def all_factors(self) -> np.ndarray:
        """Scaling factors for all 11 dimensions."""
        return np.array([self.scaling_factor(d) for d in range(11)])


def build_dimension_info() -> List[DimensionInfo]:
    """Construct physical description of all 11 dimensions."""
    scaler = DimensionalScaling()
    dims = []

    spacetime_names = ['time', 'x', 'y', 'z']
    for i in range(4):
        dims.append(DimensionInfo(
            index=i, name=spacetime_names[i], dim_type='spacetime',
            compactification_radius=float('inf'),
            kk_mass_scale=0.0,
            suppression_factor=1.0,
        ))

    for i in range(4, 11):
        config_key = i + 1
        radius = PHYS.extra_dim_bounds.get(config_key, 1e-21)
        dims.append(DimensionInfo(
            index=i, name=f'compact_{i}', dim_type='compact',
            compactification_radius=radius,
            kk_mass_scale=scaler.kk_mass(i, 1),
            suppression_factor=scaler.horava_witten_suppression(i),
        ))

    return dims


# ======================================================================
# Spectral dimension mapper
# ======================================================================

class SpectralDimensionMapper:
    """
    Maps Hamiltonian eigenvalues to 11D dimensions via energy bands.

    The key insight: Hamiltonian eigenvalues define a spectrum whose
    structure naturally encodes dimensional information. By decomposing
    the spectrum into bands and computing Frenet-Serret curvature on
    each spectral subspace, we get genuine per-dimension curvature.

    For systems with fewer eigenvalues than 11 dimensions, the
    extrapolate_compact_from_spacetime method uses geometric continuity
    to infer compact dimension states from observed spacetime curvature.
    """

    def __init__(
        self,
        H: np.ndarray,
        scaler: DimensionalScaling,
        dim_info: List[DimensionInfo],
        kk_builder=None,
    ):
        self.H = H
        self.scaler = scaler
        self.dim_info = dim_info
        self.dim_q = H.shape[0]
        self.kk_builder = kk_builder

        # Eigendecompose H once (cached -- H is fixed)
        self.eigenvalues, self.eigenvectors = np.linalg.eigh(H)
        self.E_min = float(self.eigenvalues[0])
        self.E_max = float(self.eigenvalues[-1])
        self.E_range = self.E_max - self.E_min

        # Compute KK activation energies (normalized to system scale)
        self.activation_energies = self._compute_activation_energies()

        # Compute energy bands and projectors
        if kk_builder is not None:
            self.energy_bands = kk_builder.get_natural_energy_bands()
        else:
            self.energy_bands = self._compute_energy_bands()
        self.projectors = QuantumFrenetSerret.spectral_decompose_hamiltonian(
            H, self.energy_bands
        )

    def _compute_activation_energies(self) -> np.ndarray:
        """
        KK threshold energies normalized to the system's eigenvalue range.

        Spacetime dims (0-3): threshold = -inf (always active).
        Compact dims (4-10): threshold derived from KK mass, mapped
        onto the system's energy scale via logarithmic compression.
        """
        thresholds = np.full(11, -np.inf)

        if self.E_range < 1e-20:
            return thresholds

        # KK masses span ~18 orders of magnitude (1e-13 to 1e+5 GeV)
        # Map them onto the system's eigenvalue range via log-compression
        kk_masses = []
        for d in range(4, 11):
            m = self.dim_info[d].kk_mass_scale
            kk_masses.append(m if m > 0 else 1e-20)

        log_kk = np.log10(np.array(kk_masses))
        log_min, log_max = log_kk.min(), log_kk.max()
        log_range = log_max - log_min if log_max > log_min else 1.0

        for cd in range(7):
            # Normalize KK mass to [0, 1] in log-space
            frac = (log_kk[cd] - log_min) / log_range
            # Map to eigenvalue range: lowest KK → first quartile of spectrum
            thresholds[4 + cd] = self.E_min + frac * self.E_range

        return thresholds

    def _compute_energy_bands(self) -> List[Tuple[float, float]]:
        """
        Assign eigenvalue bands to 11 dimensions.

        Divides the eigenvalue range into 11 bands using the sorted
        activation energies as boundaries. Each band corresponds to
        one dimension.
        """
        # Build band boundaries from activation energies
        boundaries = sorted(
            [self.E_min - 0.01 * max(abs(self.E_min), 1e-10)]
            + [self.activation_energies[d] for d in range(4, 11)]
            + [self.E_max + 0.01 * max(abs(self.E_max), 1e-10)]
        )

        # We need exactly 11 bands.
        # Strategy: spacetime gets 4 bands at the low end,
        # compact gets 7 bands at the high end.
        # Use linspace across eigenvalue range for even distribution.
        band_edges = np.linspace(
            self.E_min - 0.01 * max(abs(self.E_min), 1e-10),
            self.E_max + 0.01 * max(abs(self.E_max), 1e-10),
            12,  # 12 edges → 11 bands
        )

        bands = []
        for i in range(11):
            bands.append((float(band_edges[i]), float(band_edges[i + 1])))

        return bands

    def activation_function(
        self,
        system_energy: float,
        dimension: int,
        sharpness: float = 10.0,
    ) -> float:
        """
        Soft sigmoid activation for dimension accessibility.

        sigma(sharpness * (E_system - E_threshold_d))

        Spacetime dims (0-3): always 1.0.
        Compact dims (4-10): sigmoid from 0 to 1 as energy
        crosses the KK threshold.
        """
        if dimension < 4:
            return 1.0

        threshold = self.activation_energies[dimension]
        if not np.isfinite(threshold):
            return 1.0

        # Normalize energy relative to system scale
        if self.E_range > 1e-20:
            x = sharpness * (system_energy - threshold) / self.E_range
        else:
            x = 0.0

        # Numerically safe sigmoid
        x = np.clip(x, -500.0, 500.0)
        return 1.0 / (1.0 + np.exp(-x))

    def compute_dimensional_curvature(
        self,
        rho: np.ndarray,
        fs: QuantumFrenetSerret,
        system_energy: float,
        sharpness: float = 10.0,
    ) -> DimensionalCurvatureTensor:
        """
        Compute the full per-dimension curvature tensor.

        1. Spectral subspace projection: compute kappa, tau per band
        2. Activation gating: multiply by sigmoid activation
        3. Extrapolation: for empty subspaces, extrapolate from
           spacetime curvature via geometric continuity
        """
        # Per-subspace curvature from spectral decomposition
        curvatures, torsions, weights = fs.compute_subspace_curvatures(
            rho, self.H, self.projectors
        )

        # Activation mask
        activation = np.array([
            self.activation_function(system_energy, d, sharpness)
            for d in range(11)
        ])

        # Extrapolate compact dims from spacetime when subspace is empty
        # Use per-band spacetime curvature if available; otherwise fall
        # back to the total system curvature (handles the common case
        # where each band is rank-1 and gives kappa=0 individually).
        spacetime_kappa_mean = np.mean(curvatures[:4])
        if spacetime_kappa_mean < 1e-15:
            # All spectral bands are rank-1 → use total system curvature
            spacetime_kappa_mean = fs.compute_curvature(rho, self.H)

        if spacetime_kappa_mean > 1e-15:
            for d in range(11):
                if curvatures[d] < 1e-15 and (weights[d] < 1e-12 or d < 4):
                    if d < 4:
                        # Spacetime dim with rank-1 subspace: use total kappa
                        curvatures[d] = spacetime_kappa_mean * (1.0 / 4.0)
                    else:
                        # Geometric continuity: curvature transfers across
                        # dimensions with exponential decay at KK threshold
                        kk_mass = self.dim_info[d].kk_mass_scale
                        if kk_mass > 0 and system_energy > 0:
                            ratio = system_energy / kk_mass
                            curvatures[d] = spacetime_kappa_mean * (
                                1.0 - np.exp(-ratio)
                            )
                        else:
                            curvatures[d] = spacetime_kappa_mean * activation[d]

        # Apply activation gating
        curvatures_gated = curvatures * activation
        torsions_gated = torsions * activation

        return DimensionalCurvatureTensor(
            curvatures=curvatures_gated,
            torsions=torsions_gated,
            activation=activation,
            subspace_weights=weights,
        )

    def extrapolate_compact_from_spacetime(
        self,
        spacetime_curvatures: np.ndarray,
        system_energy: float,
    ) -> np.ndarray:
        """
        Extrapolate compact dimension curvature from spacetime observation.

        Uses geometric continuity: M1-11 is a single connected system.
        The curvature measured in dims 0-3 constrains dims 4-10 via:
          kappa_d = kappa_spacetime * (1 - exp(-E / E_kk[d]))

        At low energy (E << E_kk): kappa_d → 0 (dimension frozen)
        At KK scale (E ~ E_kk): kappa_d → kappa_spacetime * 0.63
        At high energy (E >> E_kk): kappa_d → kappa_spacetime (dimension open)

        Args:
            spacetime_curvatures: shape (4,) curvatures from dims 0-3.
            system_energy: current system energy.

        Returns:
            shape (7,) extrapolated curvatures for dims 4-10.
        """
        kappa_st = np.mean(np.abs(spacetime_curvatures))
        compact = np.zeros(7)

        for cd in range(7):
            dim_idx = cd + 4
            kk_mass = self.dim_info[dim_idx].kk_mass_scale

            if kk_mass > 0 and system_energy > 0:
                ratio = system_energy / kk_mass
                compact[cd] = kappa_st * (1.0 - np.exp(-ratio))
            else:
                compact[cd] = 0.0

        return compact


# ======================================================================
# Main navigation engine
# ======================================================================

class NavigationEngine:
    """
    11D navigation using the FCE for error-corrected trajectory computation.

    Integrates all physics modules:
    - Spectral decomposition for per-dimension curvature
    - KK activation thresholds for compact dimension accessibility
    - Holographic bounds (Bekenstein-Hawking, MSS, Ryu-Takayanagi)
    - Non-equilibrium thermodynamics (Jarzynski, Crooks)
    - QFT boundary effects (Casimir, Hawking temperature)
    - Chaos analysis (per-dimension Lyapunov exponents)

    Usage:
        nav = NavigationEngine(hamiltonian, system_params={...})
        result = nav.navigate(rho_init, dt=1e-6, n_steps=100)
    """

    def __init__(
        self,
        hamiltonian: np.ndarray,
        system_params: Optional[Dict[str, float]] = None,
        engine_config: Optional[EngineConfig] = None,
        nav_config: Optional[NavigationConfig] = None,
        initial_coordinates: Optional[np.ndarray] = None,
        initial_momenta: Optional[np.ndarray] = None,
        kk_builder=None,
    ):
        self.H = hamiltonian
        self.dim_q = hamiltonian.shape[0]
        self.nav_config = nav_config or NavigationConfig()
        self.kk_builder = kk_builder

        # Build FCE engine
        self.engine = FractalCorrectionEngine(
            hamiltonian=hamiltonian,
            system_params=system_params,
            config=engine_config,
        )

        # 11D structure
        self.scaler = DimensionalScaling()
        self.dim_info = build_dimension_info()
        self.scaling_factors = self.scaler.all_factors()

        # KK momentum operators (if KK builder is provided)
        self.momentum_ops = None
        if kk_builder is not None:
            self.momentum_ops = kk_builder.build_momentum_operators()

        # Spectral dimension mapper (replaces scalar projection)
        self.spectral_mapper = SpectralDimensionMapper(
            hamiltonian, self.scaler, self.dim_info,
            kk_builder=kk_builder,
        )

        # Physics modules
        self.holographic = HolographicSystem(dimensions=11)
        self.noneq = NonEquilibriumDynamics(
            temperature=self.nav_config.navigation_temperature
        )
        self.qft = QuantumFieldTheoryEffects(dimensions=4)

        self.coords = (
            initial_coordinates.copy() if initial_coordinates is not None
            else np.zeros(11)
        )
        self.momenta = (
            initial_momenta.copy() if initial_momenta is not None
            else np.zeros(11)
        )

        self.fs = QuantumFrenetSerret()

        # Transition tracking
        self.transition_state = DimensionalTransitionState(
            active_dimensions=list(range(4)),  # spacetime always active
            activation_energies=self.spectral_mapper.activation_energies.copy(),
            current_energy=0.0,
            transition_history=[],
        )

    def _update_coordinates(
        self,
        rho: np.ndarray,
        dt: float,
    ) -> Tuple[np.ndarray, np.ndarray, DimensionalCurvatureTensor]:
        """
        Update 11D coordinates using per-dimension spectral curvature.

        Returns:
            (coordinates, momenta, curvature_tensor)
        """
        coords = self.coords.copy()
        momenta = self.momenta.copy()

        # Energy expectation and speed
        E_mean = np.real(np.trace(self.H @ rho))
        H_centered = self.H - E_mean * np.eye(self.dim_q)
        E2 = np.real(np.trace(H_centered @ H_centered @ rho))
        speed = np.sqrt(max(E2, 0.0))

        # Per-dimension curvature tensor
        dim_curv = self.spectral_mapper.compute_dimensional_curvature(
            rho, self.fs, E_mean,
            sharpness=self.nav_config.kk_activation_sharpness,
        )

        # Time dimension: energy drives evolution
        momenta[0] = E_mean * self.scaling_factors[0]
        coords[0] += momenta[0] * dt

        # Spatial dimensions: speed + per-dim curvature
        for i in range(1, 4):
            momenta[i] = (
                speed * self.scaling_factors[i] * (0.5 ** i)
                + dim_curv.curvatures[i]
            )
            coords[i] += momenta[i] * dt

        # Compact dimensions
        if self.momentum_ops is not None:
            # Observable-based coordinates: coord_d = Tr(P_d @ rho) * R_d
            old_compact = coords[4:11].copy()
            for i in range(4, 11):
                d_compact = i - 4  # index into momentum_ops (0..6)
                P_d = self.momentum_ops[d_compact]
                R_d = self.kk_builder.radii[d_compact]
                # Absolute coordinate from quantum expectation value
                coords[i] = float(np.real(np.trace(P_d @ rho))) * R_d
                # Momentum = velocity = d(coord)/dt
                momenta[i] = (coords[i] - old_compact[i - 4]) / dt
        else:
            # Fallback: curvature-driven coordinates
            for i in range(4, 11):
                base_momentum = dim_curv.curvatures[i]

                # QFT boundary effects at dimensional transitions
                if self.nav_config.enable_boundary_effects:
                    activation = dim_curv.activation[i]
                    if 0.01 < activation < 0.99:
                        base_momentum += self._boundary_effect(i, activation)

                momenta[i] = base_momentum
                coords[i] += momenta[i] * dt

        self.coords = coords
        self.momenta = momenta
        return coords.copy(), momenta.copy(), dim_curv

    def _boundary_effect(self, dimension: int, activation: float) -> float:
        """
        QFT boundary effect at dimensional transition.

        Casimir energy from compactification geometry, weighted by
        transition strength (peaks at activation=0.5). Normalized
        relative to the system's energy scale so it contributes a
        perturbative correction rather than overwhelming the dynamics.
        """
        dim_info = self.dim_info[dimension]
        R = dim_info.compactification_radius

        if R == float('inf') or R <= 0:
            return 0.0

        # Casimir energy density in the compact dimension
        casimir = self.qft.casimir_energy(R, dimensions=dimension + 1)

        # Transition strength: 4*a*(1-a), peaks at 0.5
        transition_strength = 4.0 * activation * (1.0 - activation)

        # Normalize relative to system energy scale: boundary effect
        # should be a small perturbation (~1% of the curvature signal)
        E_scale = max(
            abs(self.spectral_mapper.E_range), 1e-10
        )
        if abs(casimir) > 1e-20:
            normalized = np.sign(casimir) * E_scale * 0.01
        else:
            normalized = 0.0

        return normalized * transition_strength

    def _holographic_analysis(
        self, rho: np.ndarray, energy: float
    ) -> Dict[str, float]:
        """Compute holographic quantities at this navigation step."""
        size = np.linalg.norm(self.coords[:4])
        size = max(size, 1e-10)

        bh_entropy = self.holographic.bekenstein_hawking_entropy(
            max(abs(energy), 1e-10), size
        )
        scrambling = self.holographic.information_scrambling_time(self.dim_q)

        # RT entropy: spacetime subsystem vs full system
        spacetime_size = np.linalg.norm(self.coords[:4])
        total_size = np.linalg.norm(self.coords)
        total_size = max(total_size, 1e-10)
        spacetime_size = max(spacetime_size, 1e-10)
        rt_entropy = self.holographic.entanglement_entropy_ryu_takayanagi(
            min(spacetime_size, total_size * 0.99), total_size
        )

        return {
            'bh_entropy': float(bh_entropy),
            'scrambling_time': float(scrambling),
            'rt_entanglement_entropy': float(rt_entropy),
        }

    def _qft_boundary_analysis(self, dimension: int) -> Dict[str, float]:
        """QFT effects at the boundary of a compact dimension."""
        dim_info = self.dim_info[dimension]
        R = dim_info.compactification_radius
        kk_mass = dim_info.kk_mass_scale

        casimir = self.qft.casimir_energy(
            max(R, 1e-35), dimensions=dimension + 1
        )
        hawking_T = self.qft.hawking_temperature(
            max(kk_mass, 1e-20), dimensions=dimension + 1
        )

        return {
            'casimir_energy': float(casimir),
            'hawking_temperature': float(hawking_T),
            'compactification_radius': float(R),
            'kk_mass_GeV': float(kk_mass),
        }

    def _compute_navigation_cost(
        self,
        nav_states: List[NavigationState],
        evo_result: EvolutionResult,
    ) -> NavigationCost:
        """
        Compute thermodynamic cost of the navigation trajectory.

        Work at each step = energy change of the quantum system.
        """
        energies = np.array([
            float(np.real(np.trace(self.H @ state.quantum_state)))
            for state in nav_states
        ])
        work_values = np.diff(energies)

        if len(work_values) == 0:
            return NavigationCost(
                free_energy_change=0.0,
                dissipated_work=0.0,
                entropy_production=0.0,
                landauer_cost=0.0,
                work_values=np.array([]),
                crooks_verification=None,
            )

        # Jarzynski free energy
        delta_F = self.noneq.jarzynski_equality(work_values)
        dissipated = float(np.mean(work_values)) - delta_F
        T = self.noneq.temperature
        entropy_prod = dissipated / T if T > 0 else 0.0

        # Landauer cost: information lost to decoherence
        initial_purity = float(np.real(np.trace(
            nav_states[0].quantum_state @ nav_states[0].quantum_state
        )))
        final_purity = float(np.real(np.trace(
            nav_states[-1].quantum_state @ nav_states[-1].quantum_state
        )))
        if final_purity > 1e-15 and initial_purity > 1e-15:
            bits_erased = max(
                0.0,
                np.log2(1.0 / final_purity) - np.log2(1.0 / initial_purity),
            )
        else:
            bits_erased = 0.0
        landauer_cost = T * np.log(2) * bits_erased

        # Crooks verification (forward/reverse split)
        crooks = None
        mid = len(work_values) // 2
        if mid > 5:
            try:
                crooks = self.noneq.crooks_fluctuation_theorem(
                    work_values[:mid], -work_values[mid:]
                )
            except Exception:
                pass

        return NavigationCost(
            free_energy_change=delta_F,
            dissipated_work=dissipated,
            entropy_production=entropy_prod,
            landauer_cost=landauer_cost,
            work_values=work_values,
            crooks_verification=crooks,
        )

    def _trajectory_chaos_analysis(
        self, coord_history: np.ndarray
    ) -> Dict[str, Any]:
        """
        Chaos analysis of the 11D navigation trajectory.

        Computes per-dimension Lyapunov exponents from finite differences
        of the coordinate time series.
        """
        n_steps = coord_history.shape[0]
        if n_steps < 20:
            return {
                'lyapunov_per_dim': np.zeros(11),
                'lyapunov_max': 0.0,
                'is_chaotic': False,
                'mss_bound_satisfied': True,
            }

        lyapunov_per_dim = np.zeros(11)
        for d in range(11):
            series = coord_history[:, d]
            if np.std(series) < 1e-30:
                continue
            diffs = np.abs(np.diff(series))
            diffs = diffs[diffs > 1e-30]
            if len(diffs) > 5:
                lyapunov_per_dim[d] = float(np.mean(np.log(diffs + 1e-30)))

        lyap_max = float(np.max(lyapunov_per_dim))

        # Check MSS bound: lambda <= 2*pi*T
        T = self.nav_config.navigation_temperature
        mss_bound = 2.0 * np.pi * T
        satisfies, _ = self.holographic.maldacena_shenker_stanford_bound(
            max(lyap_max, 0.0), T
        )

        return {
            'lyapunov_per_dim': lyapunov_per_dim,
            'lyapunov_max': lyap_max,
            'is_chaotic': lyap_max > 0,
            'mss_bound': float(mss_bound),
            'mss_bound_satisfied': satisfies,
        }

    def _update_transition_state(
        self, energy: float, dim_curv: DimensionalCurvatureTensor, time: float
    ):
        """Track dimensional activation/deactivation transitions."""
        # Active = activation > 0.5
        new_active = [
            d for d in range(11) if dim_curv.activation[d] > 0.5
        ]
        old_active = self.transition_state.active_dimensions

        for d in new_active:
            if d not in old_active:
                self.transition_state.transition_history.append(
                    (time, d, 'activate')
                )
        for d in old_active:
            if d not in new_active:
                self.transition_state.transition_history.append(
                    (time, d, 'deactivate')
                )

        self.transition_state.active_dimensions = new_active
        self.transition_state.current_energy = energy

    def navigate(
        self,
        rho_init: np.ndarray,
        dt: float,
        n_steps: int,
    ) -> NavigationResult:
        """
        Run 11D navigation with full physics integration.

        Args:
            rho_init: Initial quantum state density matrix
            dt: Time step
            n_steps: Number of evolution steps

        Returns:
            NavigationResult with trajectory, FCE metrics,
            per-dimension curvature tensors, transition tracking,
            thermodynamic cost, holographic metrics, and chaos analysis.
        """
        # Reset coordinates
        self.coords = np.zeros(11)
        self.momenta = np.zeros(11)
        self.transition_state.transition_history = []
        self.transition_state.active_dimensions = list(range(4))

        # Set initial compact coordinates from KK observables if available
        if self.momentum_ops is not None:
            for i in range(4, 11):
                d_compact = i - 4
                P_d = self.momentum_ops[d_compact]
                R_d = self.kk_builder.radii[d_compact]
                self.coords[i] = float(np.real(np.trace(P_d @ rho_init))) * R_d

        # Run FCE evolution
        evo_result = self.engine.evolve(rho_init, dt, n_steps)

        # Track 11D coordinates along the trajectory
        nav_states = []
        dim_curvatures = np.zeros((n_steps + 1, 11))
        curvature_tensors = []
        coord_history = np.zeros((n_steps + 1, 11))
        holographic_data = {}
        qft_effects_data = {}

        for step_idx in range(n_steps + 1):
            rho = evo_result.density_matrices[step_idx]
            t = evo_result.times[step_idx]
            E_mean = float(np.real(np.trace(self.H @ rho)))

            if step_idx > 0:
                coords, momenta, dim_curv = self._update_coordinates(rho, dt)
            else:
                coords = self.coords.copy()
                momenta = self.momenta.copy()
                # Initial curvature tensor
                dim_curv = self.spectral_mapper.compute_dimensional_curvature(
                    rho, self.fs, E_mean,
                    sharpness=self.nav_config.kk_activation_sharpness,
                )

            curvature_tensors.append(dim_curv)
            dim_curvatures[step_idx, :] = dim_curv.curvatures
            coord_history[step_idx, :] = coords

            # Update transition state
            self._update_transition_state(E_mean, dim_curv, t)

            # Periodic holographic analysis
            if step_idx % self.nav_config.holographic_check_interval == 0:
                holographic_data = self._holographic_analysis(rho, E_mean)

            # QFT boundary effects for transitioning dimensions
            for d in range(4, 11):
                if 0.01 < dim_curv.activation[d] < 0.99:
                    if d not in qft_effects_data:
                        qft_effects_data[d] = self._qft_boundary_analysis(d)

            nav_states.append(NavigationState(
                coordinates=coords.copy(),
                momenta=momenta.copy(),
                quantum_state=rho,
                time=t,
                dimension_info=self.dim_info,
            ))

        # Vulnerability analysis
        dim_vulnerabilities = np.zeros(11)
        kappa_errors_abs = np.abs(evo_result.kappa_errors)
        for d in range(11):
            dim_vulnerabilities[d] = float(
                np.mean(kappa_errors_abs) * self.scaling_factors[d]
            )
        v_max = np.max(dim_vulnerabilities)
        if v_max > 0:
            dim_vulnerabilities /= v_max

        # Thermodynamic cost
        nav_cost = self._compute_navigation_cost(nav_states, evo_result)

        # Chaos analysis
        chaos_data = self._trajectory_chaos_analysis(coord_history)

        # Build summary
        summary = self.engine.summary(evo_result)
        summary['dimension_vulnerabilities'] = {
            self.dim_info[d].name: float(dim_vulnerabilities[d])
            for d in range(11)
        }
        summary['compact_dim_suppression'] = {
            self.dim_info[d].name: float(self.scaling_factors[d])
            for d in range(4, 11)
        }
        summary['active_dimensions'] = self.transition_state.active_dimensions
        summary['n_transitions'] = len(
            self.transition_state.transition_history
        )

        return NavigationResult(
            states=nav_states,
            evolution_result=evo_result,
            dimension_curvatures=dim_curvatures,
            dimension_vulnerabilities=dim_vulnerabilities,
            total_fidelity=float(evo_result.fidelities[-1]),
            summary=summary,
            curvature_tensors=curvature_tensors,
            transition_state=self.transition_state,
            navigation_cost=nav_cost,
            holographic_metrics=holographic_data,
            chaos_metrics=chaos_data,
            qft_boundary_effects=qft_effects_data,
        )