dreamwell-engine 1.0.0

//! DreamWaveController — decoherence-based observer representation.
//!
//! The observer exists as a coherent wave form (particle cloud) that decoheres
//! into zone-specific classical forms (humanoid, vehicle, fluid) when entering
//! decoherence fields. The transition uses the quantum decoherence kernel:
//!
//!   coherence(t) = exp(-Gamma * t)
//!
//! where Gamma is the zone's decoherence rate. This is not metaphor — it is the
//! same exponential decay that governs off-diagonal density matrix elements
//! in open quantum systems (Lindblad master equation).
//!
//! Clean Compute: all state is inline arrays. Zero allocation per frame.
//! Pre-allocated bone target buffer reused every frame.

/// The observer's physical representation form.
/// Transitions between forms are governed by DecoherenceFields on zones.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum WaveForm {
    /// Coherent particle cloud — golden ratio spiral distribution.
    /// No skeleton. Pure procedural math. Cheapest representation.
    #[default]
    Particle,
    /// Humanoid skeleton — particles cluster at bone positions.
    /// FBX animation clips drive bone transforms.
    Humanoid,
    /// Rigid vehicle — particles form hull shape.
    Vehicle,
    /// Volume-preserving fluid mass.
    Fluid,
}

impl WaveForm {
    pub fn label(self) -> &'static str {
        match self {
            Self::Particle => "Particle",
            Self::Humanoid => "Humanoid",
            Self::Vehicle => "Vehicle",
            Self::Fluid => "Fluid",
        }
    }
}

/// A spatial decoherence field attached to a Zone or POI.
/// When the observer enters the field radius, the wave form begins
/// collapsing into the target form at rate gamma.
///
/// Mathematical model: Lindblad decoherence kernel
///   coherence(t) = exp(-gamma * t_in_field)
///   target_weight = 1.0 - coherence(t)
///
/// The field strength falls off spatially:
///   effective_gamma = gamma * max(0, 1 - dist/radius)
#[derive(Debug, Clone)]
pub struct DecoherenceField {
    /// Center position of the field in world space.
    pub position: [f32; 3],
    /// Field radius — decoherence begins when observer enters this radius.
    pub radius: f32,
    /// Target wave form after full collapse.
    pub target_form: WaveForm,
    /// Decoherence rate Gamma (per second). Higher = faster collapse.
    /// Typical: 2.0 (gentle, ~1s), 8.0 (fast, ~0.25s), 20.0 (instant snap).
    pub gamma: f32,
    /// Animation clip name for the target form (e.g., "climb_stairs").
    /// Empty string = use zone default idle.
    pub animation_clip: String,
    /// PropertyTag that triggered this field (for debugging/attestation).
    pub source_tag: String,
}

impl DecoherenceField {
    /// Compute effective decoherence rate at a given distance from field center.
    /// Returns 0.0 if outside radius (no decoherence).
    #[inline]
    pub fn effective_gamma(&self, observer_pos: &[f32; 3]) -> f32 {
        let dx = observer_pos[0] - self.position[0];
        let dy = observer_pos[1] - self.position[1];
        let dz = observer_pos[2] - self.position[2];
        let dist = (dx * dx + dy * dy + dz * dz).sqrt();
        if dist >= self.radius {
            return 0.0;
        }
        self.gamma * (1.0 - dist / self.radius)
    }
}

/// Tracks the observer's current decoherence state.
/// This is the density matrix diagonal — how "collapsed" the wave form is.
#[derive(Debug, Clone)]
pub struct DecoherenceState {
    /// Current coherence level [0.0, 1.0].
    /// 1.0 = fully coherent (pure particle). 0.0 = fully collapsed (target form).
    pub coherence: f32,
    /// Time spent in the current decoherence field (seconds).
    pub time_in_field: f32,
    /// Active decoherence rate (from the strongest nearby field).
    pub active_gamma: f32,
    /// Current wave form (what the observer IS right now).
    pub current_form: WaveForm,
    /// Target wave form (what the observer is COLLAPSING INTO, if any).
    pub target_form: Option<WaveForm>,
    /// Re-coherence rate when leaving a field. Default: gamma * 0.5 (slower return).
    pub recoherence_rate: f32,
    /// Whether the observer is currently inside a decoherence field.
    pub in_field: bool,
    /// Animation clip active during collapse (from the field).
    pub active_clip: String,
}

impl Default for DecoherenceState {
    fn default() -> Self {
        Self {
            coherence: 1.0,
            time_in_field: 0.0,
            active_gamma: 0.0,
            current_form: WaveForm::Particle,
            target_form: None,
            recoherence_rate: 2.0,
            in_field: false,
            active_clip: String::new(),
        }
    }
}

/// DreamWaveController — the observer's wave form representation.
/// Wraps ParticleController with zone-driven decoherence transitions.
///
/// The controller exists in one of two regimes:
/// 1. Coherent: ParticleController drives position and shape (cheap, procedural)
/// 2. Decohering/Collapsed: Target form's bone positions blend with particle positions
///
/// The transition between regimes uses the Lindblad decoherence kernel:
///   coherence(t) = exp(-Gamma * t)
///
/// Per-particle position during transition:
///   pos = lerp(particle_offset, bone_target, 1.0 - coherence)
pub struct DreamWaveController {
    /// The underlying particle controller (always exists, drives position/velocity).
    pub particle: super::particle::ParticleController,
    /// Current decoherence state.
    pub decoherence: DecoherenceState,
    /// Quantum locomotion density matrix state.
    pub quantum: QuantumLocomotionState,
    /// Per-particle target positions (bone positions for humanoid, hull points for vehicle).
    /// Pre-allocated, reused every frame (Clean Compute).
    pub bone_targets: Vec<[f32; 3]>,
    /// Per-particle coherence weight (for GPU upload).
    /// 1.0 = at particle position, 0.0 = at bone target.
    pub particle_weights: Vec<f32>,
    /// Active decoherence fields in range (sorted by distance, max 4).
    /// Reserved for future per-frame field sorting/caching.
    #[allow(dead_code)]
    active_fields: Vec<DecoherenceField>,
    /// Scratch buffer for zone field parsing.
    /// Reserved for future zone tag batch parsing.
    #[allow(dead_code)]
    field_scratch: Vec<DecoherenceField>,
}

impl DreamWaveController {
    pub fn new(position: [f32; 3], particle_count: u32) -> Self {
        Self {
            particle: super::particle::ParticleController::new(position, particle_count),
            decoherence: DecoherenceState::default(),
            quantum: QuantumLocomotionState::new(particle_count),
            bone_targets: vec![[0.0; 3]; particle_count as usize],
            particle_weights: vec![1.0; particle_count as usize],
            active_fields: Vec::with_capacity(4),
            field_scratch: Vec::with_capacity(8),
        }
    }

    /// Update the wave controller from input + zone decoherence fields.
    /// Returns true if the observer moved this frame.
    pub fn update(
        &mut self,
        forward: bool,
        back: bool,
        left: bool,
        right: bool,
        jump: bool,
        sprint: bool,
        camera_yaw: f32,
        dt: f32,
        fields: &[DecoherenceField],
    ) -> bool {
        // 1. Update underlying particle position/velocity/locomotion.
        let moved = self
            .particle
            .update(forward, back, left, right, jump, sprint, camera_yaw, dt);

        // 1b. Quantum locomotion density matrix tick.
        let moving = forward || back || left || right;
        let input_bias = [
            if !moving { 2.0 } else { 0.0 },                           // Idle bias when not moving
            if moving && !sprint { 2.0 } else { 0.0 },                 // Moving bias
            if jump { 3.0 } else { 0.0 },                              // Jump bias on press
            if self.particle.landing_timer > 0.0 { 2.0 } else { 0.0 }, // Landing
            if sprint && moving { 2.0 } else { 0.0 },                  // Sprint bias
        ];
        // Decoherence from terrain contact
        let epsilon = if self.particle.grounded { 0.3 * dt } else { 0.05 * dt };
        let _measured_mode = self.quantum.tick(dt, &input_bias, epsilon);

        // 2. Find the strongest decoherence field affecting us.
        let pos = &self.particle.position;
        let mut strongest_gamma = 0.0f32;
        let mut strongest_field: Option<&DecoherenceField> = None;
        for field in fields {
            let g = field.effective_gamma(pos);
            if g > strongest_gamma {
                strongest_gamma = g;
                strongest_field = Some(field);
            }
        }

        // 3. Apply decoherence kernel.
        if let Some(field) = strongest_field {
            // Entering or continuing in a field.
            self.decoherence.in_field = true;
            self.decoherence.active_gamma = strongest_gamma;

            // Sprint modulates collapse rate: shift = 2x faster decoherence.
            let modulated_gamma = if sprint { strongest_gamma * 2.0 } else { strongest_gamma };

            self.decoherence.time_in_field += dt;

            // Lindblad decoherence kernel: coherence decays exponentially.
            // coherence(t) = exp(-Gamma * t)
            self.decoherence.coherence = (-modulated_gamma * self.decoherence.time_in_field).exp();

            // Set target form from field.
            if self.decoherence.target_form.is_none() || self.decoherence.target_form != Some(field.target_form) {
                self.decoherence.target_form = Some(field.target_form);
                self.decoherence.active_clip = field.animation_clip.clone();
            }

            // Fully collapsed threshold.
            if self.decoherence.coherence < 0.01 {
                self.decoherence.coherence = 0.0;
                self.decoherence.current_form = field.target_form;
            }
        } else {
            // Outside all fields — re-cohere back to particle.
            self.decoherence.in_field = false;
            self.decoherence.time_in_field = 0.0;

            if self.decoherence.coherence < 1.0 {
                // Re-coherence: exponential return to particle form.
                // recoherence(t) = 1 - exp(-recoherence_rate * dt)
                let rate = self.decoherence.recoherence_rate;
                self.decoherence.coherence += (1.0 - self.decoherence.coherence) * (1.0 - (-rate * dt).exp());

                if self.decoherence.coherence > 0.99 {
                    self.decoherence.coherence = 1.0;
                    self.decoherence.current_form = WaveForm::Particle;
                    self.decoherence.target_form = None;
                    self.decoherence.active_clip.clear();
                }
            }
        }

        // 4. Update per-particle weights for GPU.
        let w = self.decoherence.coherence;
        for pw in &mut self.particle_weights {
            *pw = w;
        }

        moved
    }

    /// Current coherence level (1.0 = fully particle, 0.0 = fully collapsed).
    pub fn coherence(&self) -> f32 {
        self.decoherence.coherence
    }

    /// Current wave form.
    pub fn current_form(&self) -> WaveForm {
        self.decoherence.current_form
    }

    /// Whether the observer is actively decohering.
    pub fn is_decohering(&self) -> bool {
        self.decoherence.in_field && self.decoherence.coherence > 0.01
    }

    /// Whether the observer is fully collapsed into a target form.
    pub fn is_collapsed(&self) -> bool {
        self.decoherence.coherence < 0.01
    }

    /// Whether the observer is re-cohering back to particle.
    pub fn is_recohering(&self) -> bool {
        !self.decoherence.in_field && self.decoherence.coherence < 0.99
    }

    /// Set bone target positions for the current target form.
    /// Called when FBX skeleton is evaluated.
    pub fn set_bone_targets(&mut self, bones: &[[f32; 3]]) {
        // Map particles to nearest bones via spatial proximity.
        let particle_count = self.particle.particle_count as usize;
        let offsets = super::particle::particle_sphere_offsets(self.particle.particle_count);
        let params = self.particle.shape_params();
        let radius = self.particle.cloud_radius * params.radius_scale;

        for i in 0..particle_count.min(self.bone_targets.len()) {
            if i >= offsets.len() {
                break;
            }
            // Particle world position (particle distribution).
            let px = self.particle.position[0] + offsets[i][0] * params.stretch[0] * radius;
            let py = self.particle.position[1] + offsets[i][1] * params.stretch[1] * radius;
            let pz = self.particle.position[2] + offsets[i][2] * params.stretch[2] * radius;

            // Find nearest bone.
            let mut best_dist_sq = f32::MAX;
            let mut best_bone = [0.0f32; 3];
            for bone in bones {
                let dx = px - bone[0];
                let dy = py - bone[1];
                let dz = pz - bone[2];
                let d2 = dx * dx + dy * dy + dz * dz;
                if d2 < best_dist_sq {
                    best_dist_sq = d2;
                    best_bone = *bone;
                }
            }
            self.bone_targets[i] = best_bone;
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════════
// Quantum Locomotion — density matrix, Hamiltonian evolution, dephasing,
// Born rule measurement, interference kernels.
//
// The observer's locomotion state is modeled as a 5-dimensional Hilbert space
// H_loc ≅ C^5 with basis states |Idle⟩, |Moving⟩, |Jumping⟩, |Landing⟩, |Sprinting⟩.
//
// The density matrix ρ carries both mode populations (diagonal) and inter-mode
// coherence (off-diagonal). Coherent evolution via a Hamiltonian generates
// interference that visibly affects particle geometry — this is not metaphor.
// ═══════════════════════════════════════════════════════════════════════════════

/// Complex number for density matrix operations. No external deps.
#[derive(Debug, Clone, Copy, Default)]
pub struct Complex {
    pub re: f32,
    pub im: f32,
}

impl Complex {
    pub const ZERO: Self = Self { re: 0.0, im: 0.0 };
    pub const ONE: Self = Self { re: 1.0, im: 0.0 };

    pub fn new(re: f32, im: f32) -> Self {
        Self { re, im }
    }
    pub fn from_real(r: f32) -> Self {
        Self { re: r, im: 0.0 }
    }

    pub fn conj(self) -> Self {
        Self {
            re: self.re,
            im: -self.im,
        }
    }
    pub fn norm_sq(self) -> f32 {
        self.re * self.re + self.im * self.im
    }
    pub fn norm(self) -> f32 {
        self.norm_sq().sqrt()
    }

    pub fn mul(self, other: Self) -> Self {
        Self {
            re: self.re * other.re - self.im * other.im,
            im: self.re * other.im + self.im * other.re,
        }
    }

    pub fn add(self, other: Self) -> Self {
        Self {
            re: self.re + other.re,
            im: self.im + other.im,
        }
    }

    pub fn sub(self, other: Self) -> Self {
        Self {
            re: self.re - other.re,
            im: self.im - other.im,
        }
    }

    pub fn scale(self, s: f32) -> Self {
        Self {
            re: self.re * s,
            im: self.im * s,
        }
    }

    /// e^{i*theta} = cos(theta) + i*sin(theta)
    #[allow(dead_code)]
    pub fn exp_i(theta: f32) -> Self {
        Self {
            re: theta.cos(),
            im: theta.sin(),
        }
    }
}

/// 5x5 density matrix for the locomotion Hilbert space H_loc = C^5.
/// Basis: |I>=Idle, |M>=Moving, |J>=Jumping, |L>=Landing, |S>=Sprinting.
///
/// Hermitian: rho_dag = rho (entries[i][j] = conj(entries[j][i]))
/// Positive semidefinite: all eigenvalues >= 0
/// Trace 1: sum_k rho_kk = 1
///
/// This is the standard density matrix from quantum information theory.
/// Diagonal entries are mode populations. Off-diagonal entries carry coherence.
pub struct DensityMatrix5 {
    /// Flat 5x5 complex entries. entries[i*5+j] = rho_{ij}.
    pub entries: [Complex; 25],
}

impl DensityMatrix5 {
    /// Pure state |k><k| — all population in mode k.
    pub fn pure_state(k: usize) -> Self {
        let mut entries = [Complex::ZERO; 25];
        if k < 5 {
            entries[k * 5 + k] = Complex::ONE;
        }
        Self { entries }
    }

    /// Create from a state vector (5 complex amplitudes).
    /// rho = |psi><psi| where |psi> = sum alpha_k |k>
    pub fn from_state_vector(amplitudes: &[Complex; 5]) -> Self {
        let mut entries = [Complex::ZERO; 25];
        for i in 0..5 {
            for j in 0..5 {
                entries[i * 5 + j] = amplitudes[i].mul(amplitudes[j].conj());
            }
        }
        Self { entries }
    }

    /// Trace: sum_k rho_kk
    pub fn trace(&self) -> f32 {
        (0..5).map(|k| self.entries[k * 5 + k].re).sum()
    }

    /// Diagonal populations (mode probabilities).
    pub fn populations(&self) -> [f32; 5] {
        let mut p = [0.0f32; 5];
        for k in 0..5 {
            p[k] = self.entries[k * 5 + k].re;
        }
        p
    }

    /// Total off-diagonal coherence magnitude: sum_{i<j} |rho_{ij}|
    pub fn coherence_magnitude(&self) -> f32 {
        let mut c = 0.0f32;
        for i in 0..5 {
            for j in (i + 1)..5 {
                c += self.entries[i * 5 + j].norm();
            }
        }
        c
    }

    /// Completely dephase: zero all off-diagonals. rho -> Delta(rho).
    pub fn dephased(&self) -> Self {
        let mut result = Self {
            entries: [Complex::ZERO; 25],
        };
        for k in 0..5 {
            result.entries[k * 5 + k] = self.entries[k * 5 + k];
        }
        result
    }

    /// Partial dephasing channel: rho -> (1-eps)rho + eps*Delta(rho).
    /// eps in [0,1]. eps=0: no change. eps=1: full dephasing.
    pub fn apply_dephasing(&mut self, epsilon: f32) {
        let retain = 1.0 - epsilon.clamp(0.0, 1.0);
        for i in 0..5 {
            for j in 0..5 {
                if i != j {
                    self.entries[i * 5 + j] = self.entries[i * 5 + j].scale(retain);
                }
            }
        }
    }

    /// Apply unitary evolution: rho -> U * rho * U_dag.
    /// U is a 5x5 complex matrix (flat, row-major).
    pub fn apply_unitary(&mut self, u: &[Complex; 25]) {
        // Compute U * rho
        let mut u_rho = [Complex::ZERO; 25];
        for i in 0..5 {
            for j in 0..5 {
                let mut sum = Complex::ZERO;
                for k in 0..5 {
                    sum = sum.add(u[i * 5 + k].mul(self.entries[k * 5 + j]));
                }
                u_rho[i * 5 + j] = sum;
            }
        }
        // Compute (U * rho) * U_dag
        for i in 0..5 {
            for j in 0..5 {
                let mut sum = Complex::ZERO;
                for k in 0..5 {
                    sum = sum.add(u_rho[i * 5 + k].mul(u[j * 5 + k].conj()));
                }
                self.entries[i * 5 + j] = sum;
            }
        }
    }

    /// Born rule measurement: p_k = rho_kk.
    /// Returns mode index sampled from probability distribution.
    /// Uses deterministic seed for reproducible gameplay.
    pub fn measure(&self, seed: u64) -> usize {
        let pops = self.populations();
        // Deterministic pseudo-random from seed.
        let mut s = seed;
        s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        let r = (s >> 33) as f32 / (1u64 << 31) as f32; // [0, 1)
        let mut cumulative = 0.0;
        for (i, &p) in pops.iter().enumerate() {
            cumulative += p;
            if r < cumulative {
                return i;
            }
        }
        4 // fallback: Sprint
    }
}

/// Locomotion Hamiltonian H_loc = H_bias + H_input + H_couple + H_phase.
/// Generates the unitary U_tick = exp(-i*H*dt/hbar_eff) each frame.
pub struct LocomotionHamiltonian {
    /// Effective reduced Planck constant (engine tuning parameter).
    pub hbar_eff: f32,
    /// Bias energies per mode [E_I, E_M, E_J, E_L, E_S].
    pub bias: [f32; 5],
    /// Coupling strengths between adjacent modes.
    /// Order: omega_IM, omega_MS, omega_MJ, omega_JL, omega_LI
    pub couplings: [f32; 5],
}

impl Default for LocomotionHamiltonian {
    fn default() -> Self {
        Self {
            hbar_eff: 1.0,
            bias: [0.0, 0.2, 0.5, 0.3, 0.4],      // Idle lowest energy
            couplings: [1.0, 0.8, 0.6, 0.9, 0.7], // IM, MS, MJ, JL, LI
        }
    }
}

impl LocomotionHamiltonian {
    /// Build the 5x5 Hamiltonian matrix from current parameters + input.
    /// input_bias: per-mode input contribution [b_I, b_M, b_J, b_L, b_S].
    pub fn build_matrix(&self, input_bias: &[f32; 5], phase_time: f32) -> [Complex; 25] {
        let mut h = [Complex::ZERO; 25];

        // Diagonal: bias + input + phase drift
        for k in 0..5 {
            h[k * 5 + k] = Complex::from_real(self.bias[k] + input_bias[k] + 0.1 * (k as f32 * phase_time).sin());
        }

        // Off-diagonal couplings (Hermitian: h_ij = h_ji*)
        // Coupling pairs: (I,M)=0, (M,S)=1, (M,J)=2, (J,L)=3, (L,I)=4
        let pairs: [(usize, usize); 5] = [(0, 1), (1, 4), (1, 2), (2, 3), (3, 0)];
        for (idx, &(i, j)) in pairs.iter().enumerate() {
            let w = Complex::from_real(self.couplings[idx]);
            h[i * 5 + j] = w;
            h[j * 5 + i] = w.conj(); // Hermitian
        }

        h
    }

    /// Compute U = exp(-i*H*dt/hbar_eff) via second-order Taylor expansion.
    /// U = I - i*H*dt/hbar + (1/2)(-i*H*dt/hbar)^2
    /// Sufficient for small dt.
    pub fn compute_unitary(&self, h: &[Complex; 25], dt: f32) -> [Complex; 25] {
        let scale = dt / self.hbar_eff;
        let mut u = [Complex::ZERO; 25];

        // U = I
        for k in 0..5 {
            u[k * 5 + k] = Complex::ONE;
        }

        // U -= i*H*scale  (first order)
        for i in 0..5 {
            for j in 0..5 {
                // -i * h_ij * scale
                let ih = Complex::new(h[i * 5 + j].im * scale, -h[i * 5 + j].re * scale);
                u[i * 5 + j] = u[i * 5 + j].sub(ih);
            }
        }

        // U += (1/2)(-i*H*scale)^2 = -(1/2)*H^2*scale^2  (second order)
        let s2 = scale * scale * 0.5;
        for i in 0..5 {
            for j in 0..5 {
                let mut h2_ij = Complex::ZERO;
                for k in 0..5 {
                    h2_ij = h2_ij.add(h[i * 5 + k].mul(h[k * 5 + j]));
                }
                u[i * 5 + j] = u[i * 5 + j].sub(h2_ij.scale(s2));
            }
        }

        u
    }
}

/// Interference kernel kappa_{n,ij} in C^3 for particle n between modes i and j.
/// Controls how off-diagonal coherence rho_{ij} affects particle geometry.
///
/// The particle position with interference is:
///   x_n = root + sum_k rho_kk * g_{n,k} + 2*sum_{i<j} Re[kappa_{n,ij} * rho_{ij}]
///
/// This is the decisive test from QUANTUM.md S7:
///   If exists n: x_n(rho) != x_n(Delta(rho)), then coherence is not metaphor — it's simulation.
pub struct InterferenceKernels {
    /// For each mode pair (i,j) with i<j, a 3D complex displacement.
    /// 10 pairs x 3 axes x 2 (re,im) = 60 f32 values per particle.
    /// We store kernels for a representative subset of particles.
    /// Pair order: (0,1),(0,2),(0,3),(0,4),(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)
    pub kernels: Vec<[[Complex; 3]; 10]>,
}

impl InterferenceKernels {
    /// Generate interference kernels for N particles.
    /// Each kernel is designed to create visible wave-like distortion
    /// when the corresponding mode coherence is nonzero.
    pub fn generate(particle_count: u32) -> Self {
        let n = particle_count as usize;
        let mut kernels = Vec::with_capacity(n);
        let golden = (1.0 + 5.0f32.sqrt()) / 2.0;

        for p in 0..n {
            let t = p as f32 / n.max(1) as f32;
            let mut pair_kernels = [[Complex::ZERO; 3]; 10];

            // Each mode pair gets a distinct interference pattern.
            // The pattern creates wave-like ripples, spirals, or oscillations
            // that are ONLY visible when off-diagonal coherence is nonzero.

            // (I,M) pair 0: Idle<->Moving creates horizontal wave ripple
            pair_kernels[0] = [
                Complex::new(0.3 * (t * 12.0).sin(), 0.1 * (t * 8.0).cos()),
                Complex::new(0.0, 0.15 * (t * 10.0).sin()),
                Complex::new(0.2 * (t * 6.0).cos(), 0.0),
            ];

            // (I,J) pair 1: Idle<->Jump creates vertical expansion pulse
            pair_kernels[1] = [
                Complex::new(0.0, 0.1 * (t * 4.0).sin()),
                Complex::new(0.4 * (t * 8.0).sin(), 0.2 * (t * golden * 5.0).cos()),
                Complex::new(0.0, 0.1 * (t * 4.0).cos()),
            ];

            // (I,L) pair 2: Idle<->Land creates compression ripple
            pair_kernels[2] = [
                Complex::new(0.15 * (t * 16.0).sin(), 0.0),
                Complex::new(-0.3 * (t * 12.0).cos(), 0.1),
                Complex::new(0.15 * (t * 16.0).cos(), 0.0),
            ];

            // (I,S) pair 3: Idle<->Sprint creates elongation wave
            pair_kernels[3] = [
                Complex::new(0.1, 0.2 * (t * 20.0).sin()),
                Complex::new(0.0, 0.05),
                Complex::new(0.35 * (t * 14.0).sin(), 0.15 * (t * 7.0).cos()),
            ];

            // (M,J) pair 4: Move<->Jump creates arc transition
            pair_kernels[4] = [
                Complex::new(0.1 * (t * 6.0).cos(), 0.0),
                Complex::new(0.25 * (t * 10.0).sin(), 0.2 * (t * 5.0).cos()),
                Complex::new(0.15 * (t * 8.0).sin(), 0.0),
            ];

            // (M,L) pair 5: Move<->Land creates ground-impact spread
            pair_kernels[5] = [
                Complex::new(0.2 * (t * 14.0).sin(), 0.1 * (t * 7.0).cos()),
                Complex::new(-0.15 * (t * 18.0).cos(), 0.0),
                Complex::new(0.2 * (t * 14.0).cos(), -0.1 * (t * 7.0).sin()),
            ];

            // (M,S) pair 6: Move<->Sprint creates velocity streak
            pair_kernels[6] = [
                Complex::new(0.05, 0.1 * (t * 22.0).sin()),
                Complex::new(0.0, 0.05),
                Complex::new(0.3 * (t * 18.0).sin(), 0.2 * (t * 9.0).cos()),
            ];

            // (J,L) pair 7: Jump<->Land creates bounce oscillation
            pair_kernels[7] = [
                Complex::new(0.1 * (t * 8.0).cos(), 0.0),
                Complex::new(0.35 * (t * 16.0).sin(), -0.15 * (t * 8.0).cos()),
                Complex::new(0.1 * (t * 8.0).sin(), 0.0),
            ];

            // (J,S) pair 8: Jump<->Sprint creates aerial streak
            pair_kernels[8] = [
                Complex::new(0.15 * (t * 10.0).sin(), 0.1),
                Complex::new(0.2 * (t * 12.0).cos(), 0.15 * (t * 6.0).sin()),
                Complex::new(0.25 * (t * 15.0).sin(), 0.0),
            ];

            // (L,S) pair 9: Land<->Sprint creates recovery burst
            pair_kernels[9] = [
                Complex::new(0.2 * (t * 20.0).sin(), 0.1 * (t * 10.0).cos()),
                Complex::new(-0.1 * (t * 24.0).cos(), 0.05),
                Complex::new(0.15 * (t * 16.0).cos(), 0.15 * (t * 8.0).sin()),
            ];

            kernels.push(pair_kernels);
        }

        Self { kernels }
    }

    /// Compute interference displacement for particle n from density matrix.
    /// Returns the 3D displacement caused by off-diagonal coherence.
    /// This is: 2*sum_{i<j} Re[kappa_{n,ij} * rho_{ij}]
    pub fn interference_displacement(&self, n: usize, rho: &DensityMatrix5) -> [f32; 3] {
        if n >= self.kernels.len() {
            return [0.0; 3];
        }
        let k = &self.kernels[n];
        let mut disp = [0.0f32; 3];

        // Iterate over all 10 pairs (i<j)
        let pairs: [(usize, usize); 10] = [
            (0, 1),
            (0, 2),
            (0, 3),
            (0, 4),
            (1, 2),
            (1, 3),
            (1, 4),
            (2, 3),
            (2, 4),
            (3, 4),
        ];
        for (idx, &(i, j)) in pairs.iter().enumerate() {
            let rho_ij = rho.entries[i * 5 + j];
            for axis in 0..3 {
                // 2 * Re[kappa * rho_ij]
                let product = k[idx][axis].mul(rho_ij);
                disp[axis] += 2.0 * product.re;
            }
        }

        disp
    }
}

/// Full quantum locomotion state — wraps density matrix, Hamiltonian, kernels.
/// This is the quantum simulation subsystem that satisfies QUANTUM.md S11.
pub struct QuantumLocomotionState {
    /// The density matrix rho in C^{5x5}.
    pub rho: DensityMatrix5,
    /// The Hamiltonian governing coherent evolution.
    pub hamiltonian: LocomotionHamiltonian,
    /// Interference kernels for particle rendering.
    pub kernels: InterferenceKernels,
    /// Cumulative phase time (for H_phase).
    pub phase_time: f32,
    /// Frame counter for measurement seed.
    pub frame: u64,
}

impl QuantumLocomotionState {
    pub fn new(particle_count: u32) -> Self {
        // Initialize in Idle pure state: rho = |I><I|
        Self {
            rho: DensityMatrix5::pure_state(0),
            hamiltonian: LocomotionHamiltonian::default(),
            kernels: InterferenceKernels::generate(particle_count),
            phase_time: 0.0,
            frame: 0,
        }
    }

    /// Tick the quantum locomotion state.
    /// 1. Build input bias from player state
    /// 2. Hamiltonian evolution (unitary)
    /// 3. Decoherence channel (terrain/contact)
    /// 4. Returns Born-rule measured mode
    pub fn tick(&mut self, dt: f32, input_bias: &[f32; 5], decoherence_epsilon: f32) -> usize {
        self.phase_time += dt;
        self.frame += 1;

        // 1. Build Hamiltonian
        let h = self.hamiltonian.build_matrix(input_bias, self.phase_time);

        // 2. Unitary evolution: rho -> U*rho*U_dag
        let u = self.hamiltonian.compute_unitary(&h, dt);
        self.rho.apply_unitary(&u);

        // 3. Decoherence channel: rho -> (1-eps)rho + eps*Delta(rho)
        self.rho.apply_dephasing(decoherence_epsilon);

        // 4. Born rule measurement
        self.rho.measure(self.frame)
    }

    /// Compute particle position with interference.
    /// pos = root + sum_k rho_kk * g_{n,k} + interference_displacement(n)
    pub fn particle_position_with_interference(
        &self,
        n: usize,
        root: [f32; 3],
        mode_templates: &[[f32; 3]; 5],
    ) -> [f32; 3] {
        let pops = self.rho.populations();
        let interference = self.kernels.interference_displacement(n, &self.rho);

        let mut pos = root;
        for k in 0..5 {
            pos[0] += pops[k] * mode_templates[k][0];
            pos[1] += pops[k] * mode_templates[k][1];
            pos[2] += pops[k] * mode_templates[k][2];
        }
        pos[0] += interference[0];
        pos[1] += interference[1];
        pos[2] += interference[2];

        pos
    }
}

/// Parsed result from a decoherence PropertyTag string.
#[derive(Debug, Clone, PartialEq)]
pub enum DecoherenceTagValue {
    /// A form tag: target wave form + optional animation clip name.
    Form(WaveForm, String),
    /// A gamma override tag: decoherence rate value.
    Gamma(f32),
}

/// Parse a decoherence PropertyTag string.
/// Format: "decohere:form[:clip]" or "gamma:rate"
/// Examples:
///   "decohere:humanoid:climb_stairs" -> (Humanoid, "climb_stairs")
///   "decohere:fluid" -> (Fluid, "")
///   "gamma:4.0" -> Gamma(4.0)
pub fn parse_decoherence_tag(tag: &str) -> Option<DecoherenceTagValue> {
    if let Some(rest) = tag.strip_prefix("decohere:") {
        let mut parts = rest.splitn(2, ':');
        let form_str = parts.next()?;
        let form = match form_str {
            "humanoid" => WaveForm::Humanoid,
            "vehicle" => WaveForm::Vehicle,
            "fluid" => WaveForm::Fluid,
            "particle" => WaveForm::Particle,
            _ => return None,
        };
        let clip = parts.next().unwrap_or("").to_string();
        Some(DecoherenceTagValue::Form(form, clip))
    } else if let Some(rest) = tag.strip_prefix("gamma:") {
        let rate: f32 = rest.parse().ok()?;
        Some(DecoherenceTagValue::Gamma(rate))
    } else {
        None
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn wave_form_default_is_particle() {
        assert_eq!(WaveForm::default(), WaveForm::Particle);
    }

    #[test]
    fn decoherence_field_outside_radius_zero() {
        let field = DecoherenceField {
            position: [0.0, 0.0, 0.0],
            radius: 5.0,
            target_form: WaveForm::Humanoid,
            gamma: 4.0,
            animation_clip: String::new(),
            source_tag: String::new(),
        };
        // Observer at distance 10 — well outside radius 5.
        let g = field.effective_gamma(&[10.0, 0.0, 0.0]);
        assert_eq!(g, 0.0, "gamma should be 0 outside field radius");
    }

    #[test]
    fn decoherence_field_at_center_full_gamma() {
        let field = DecoherenceField {
            position: [5.0, 0.0, 5.0],
            radius: 10.0,
            target_form: WaveForm::Humanoid,
            gamma: 8.0,
            animation_clip: String::new(),
            source_tag: String::new(),
        };
        // Observer at the exact center of the field.
        let g = field.effective_gamma(&[5.0, 0.0, 5.0]);
        assert!(
            (g - 8.0).abs() < 1e-5,
            "gamma at center should equal field gamma (8.0), got {}",
            g
        );
    }

    #[test]
    fn decoherence_field_at_edge_near_zero() {
        let field = DecoherenceField {
            position: [0.0, 0.0, 0.0],
            radius: 10.0,
            target_form: WaveForm::Fluid,
            gamma: 4.0,
            animation_clip: String::new(),
            source_tag: String::new(),
        };
        // Observer at distance 9.99 (just inside edge).
        let g = field.effective_gamma(&[9.99, 0.0, 0.0]);
        assert!(g < 0.01, "gamma near the edge should be near zero, got {}", g);
        assert!(g > 0.0, "gamma just inside edge should be >0, got {}", g);
    }

    #[test]
    fn dream_wave_controller_starts_coherent() {
        let ctrl = DreamWaveController::new([0.0, 1.0, 0.0], 64);
        assert!(
            (ctrl.coherence() - 1.0).abs() < 1e-5,
            "new controller should start fully coherent (1.0), got {}",
            ctrl.coherence()
        );
        assert_eq!(ctrl.current_form(), WaveForm::Particle);
        assert!(!ctrl.is_decohering());
        assert!(!ctrl.is_collapsed());
        assert!(!ctrl.is_recohering());
    }

    #[test]
    fn dream_wave_controller_decoheres_in_field() {
        let mut ctrl = DreamWaveController::new([0.0, 1.0, 0.0], 64);
        let fields = vec![DecoherenceField {
            position: [0.0, 1.0, 0.0],
            radius: 10.0,
            target_form: WaveForm::Humanoid,
            gamma: 4.0,
            animation_clip: "idle".into(),
            source_tag: "decohere:humanoid".into(),
        }];

        // Run several frames inside the field.
        for _ in 0..60 {
            ctrl.update(false, false, false, false, false, false, 0.0, 0.016, &fields);
        }

        assert!(
            ctrl.coherence() < 1.0,
            "coherence should decrease inside a field, got {}",
            ctrl.coherence()
        );
        assert!(ctrl.decoherence.in_field, "should be marked as in_field");
        assert_eq!(
            ctrl.decoherence.target_form,
            Some(WaveForm::Humanoid),
            "target form should be Humanoid"
        );
    }

    #[test]
    fn dream_wave_controller_recoheres_outside_field() {
        let mut ctrl = DreamWaveController::new([0.0, 1.0, 0.0], 64);
        let fields = vec![DecoherenceField {
            position: [0.0, 1.0, 0.0],
            radius: 10.0,
            target_form: WaveForm::Humanoid,
            gamma: 8.0,
            animation_clip: String::new(),
            source_tag: String::new(),
        }];

        // Decohere for a while.
        for _ in 0..120 {
            ctrl.update(false, false, false, false, false, false, 0.0, 0.016, &fields);
        }
        let coherence_in_field = ctrl.coherence();
        assert!(
            coherence_in_field < 0.5,
            "should be significantly decohered after 120 frames at gamma=8, got {}",
            coherence_in_field
        );

        // Now update outside any field — re-coherence should begin.
        let no_fields: Vec<DecoherenceField> = vec![];
        for _ in 0..300 {
            ctrl.update(false, false, false, false, false, false, 0.0, 0.016, &no_fields);
        }

        assert!(
            ctrl.coherence() > coherence_in_field,
            "coherence should increase after leaving field: was {}, now {}",
            coherence_in_field,
            ctrl.coherence()
        );
        assert!(!ctrl.decoherence.in_field, "should not be in_field after leaving");
    }

    #[test]
    fn parse_decoherence_tag_humanoid() {
        let result = parse_decoherence_tag("decohere:humanoid:climb_stairs");
        assert_eq!(
            result,
            Some(DecoherenceTagValue::Form(WaveForm::Humanoid, "climb_stairs".into()))
        );
    }

    #[test]
    fn parse_decoherence_tag_gamma() {
        let result = parse_decoherence_tag("gamma:4.0");
        assert_eq!(result, Some(DecoherenceTagValue::Gamma(4.0)));
    }

    // ── Additional coverage ────────────────────────────────────────────

    #[test]
    fn parse_decoherence_tag_fluid_no_clip() {
        let result = parse_decoherence_tag("decohere:fluid");
        assert_eq!(result, Some(DecoherenceTagValue::Form(WaveForm::Fluid, String::new())));
    }

    #[test]
    fn parse_decoherence_tag_unknown_form_returns_none() {
        assert_eq!(parse_decoherence_tag("decohere:dragon"), None);
    }

    #[test]
    fn parse_decoherence_tag_unrelated_returns_none() {
        assert_eq!(parse_decoherence_tag("some_other_tag"), None);
    }

    #[test]
    fn parse_decoherence_tag_gamma_invalid_returns_none() {
        assert_eq!(parse_decoherence_tag("gamma:notanumber"), None);
    }

    #[test]
    fn wave_form_labels() {
        assert_eq!(WaveForm::Particle.label(), "Particle");
        assert_eq!(WaveForm::Humanoid.label(), "Humanoid");
        assert_eq!(WaveForm::Vehicle.label(), "Vehicle");
        assert_eq!(WaveForm::Fluid.label(), "Fluid");
    }

    #[test]
    fn dream_wave_controller_fully_collapses() {
        let mut ctrl = DreamWaveController::new([0.0, 1.0, 0.0], 32);
        let fields = vec![DecoherenceField {
            position: [0.0, 1.0, 0.0],
            radius: 10.0,
            target_form: WaveForm::Vehicle,
            gamma: 20.0, // very fast
            animation_clip: String::new(),
            source_tag: String::new(),
        }];

        // Run many frames with high gamma — should fully collapse.
        for _ in 0..200 {
            ctrl.update(false, false, false, false, false, false, 0.0, 0.016, &fields);
        }

        assert!(ctrl.is_collapsed(), "should be fully collapsed");
        assert_eq!(ctrl.current_form(), WaveForm::Vehicle);
        assert_eq!(ctrl.coherence(), 0.0);
    }

    #[test]
    fn dream_wave_controller_particle_weights_match_coherence() {
        let mut ctrl = DreamWaveController::new([0.0, 1.0, 0.0], 16);
        let fields = vec![DecoherenceField {
            position: [0.0, 1.0, 0.0],
            radius: 10.0,
            target_form: WaveForm::Humanoid,
            gamma: 4.0,
            animation_clip: String::new(),
            source_tag: String::new(),
        }];

        ctrl.update(false, false, false, false, false, false, 0.0, 0.016, &fields);
        let c = ctrl.coherence();
        for &w in &ctrl.particle_weights {
            assert!(
                (w - c).abs() < 1e-5,
                "particle weight {} should match coherence {}",
                w,
                c
            );
        }
    }

    // ── Quantum Locomotion Ablation Tests ──────────────────────────────

    #[test]
    fn ablation_interference_changes_output() {
        // THE DECISIVE TEST: x_n(rho) != x_n(Delta(rho)) for some n.
        let mut state = QuantumLocomotionState::new(64);

        // Create a superposition state with nonzero off-diagonals.
        let sq2 = 1.0 / 2.0f32.sqrt();
        state.rho = DensityMatrix5::from_state_vector(&[
            Complex::new(sq2, 0.0), // Idle
            Complex::new(sq2, 0.0), // Moving
            Complex::ZERO,
            Complex::ZERO,
            Complex::ZERO,
        ]);

        let root = [0.0; 3];
        let templates: [[f32; 3]; 5] = [
            [0.0, 0.0, 0.0],  // Idle
            [0.0, 0.0, 1.0],  // Moving
            [0.0, 1.0, 0.0],  // Jump
            [0.0, -0.5, 0.0], // Land
            [0.0, 0.0, 2.0],  // Sprint
        ];

        // Full quantum position
        let pos_quantum = state.particle_position_with_interference(0, root, &templates);

        // Dephased position (remove off-diagonals)
        let dephased = state.rho.dephased();
        let old_rho_entries = state.rho.entries;
        state.rho = dephased;
        let pos_dephased = state.particle_position_with_interference(0, root, &templates);
        state.rho.entries = old_rho_entries;

        // THE TEST: they must differ
        let diff_sq = (pos_quantum[0] - pos_dephased[0]).powi(2)
            + (pos_quantum[1] - pos_dephased[1]).powi(2)
            + (pos_quantum[2] - pos_dephased[2]).powi(2);

        assert!(
            diff_sq > 1e-6,
            "QUANTUM.md S7 FAILED: interference does not affect output. diff_sq={diff_sq}"
        );
    }

    #[test]
    fn density_matrix_trace_one() {
        let rho = DensityMatrix5::pure_state(2);
        assert!((rho.trace() - 1.0).abs() < 1e-6);
    }

    #[test]
    fn density_matrix_from_state_vector_trace() {
        let sq2 = 1.0 / 2.0f32.sqrt();
        let rho = DensityMatrix5::from_state_vector(&[
            Complex::new(sq2, 0.0),
            Complex::new(sq2, 0.0),
            Complex::ZERO,
            Complex::ZERO,
            Complex::ZERO,
        ]);
        assert!((rho.trace() - 1.0).abs() < 1e-6);
        // Off-diagonal should be nonzero
        assert!(rho.coherence_magnitude() > 0.1);
    }

    #[test]
    fn dephasing_removes_coherence() {
        let sq2 = 1.0 / 2.0f32.sqrt();
        let rho = DensityMatrix5::from_state_vector(&[
            Complex::new(sq2, 0.0),
            Complex::new(sq2, 0.0),
            Complex::ZERO,
            Complex::ZERO,
            Complex::ZERO,
        ]);
        let dephased = rho.dephased();
        assert!(dephased.coherence_magnitude() < 1e-6);
        assert!((dephased.trace() - 1.0).abs() < 1e-6);
    }

    #[test]
    fn hamiltonian_evolution_preserves_trace() {
        let mut state = QuantumLocomotionState::new(16);
        let input = [1.0, 0.5, 0.0, 0.0, 0.0];
        state.tick(1.0 / 60.0, &input, 0.0);
        // Trace should remain ~1.0 (small numerical drift acceptable)
        assert!(
            (state.rho.trace() - 1.0).abs() < 0.1,
            "Trace drifted: {}",
            state.rho.trace()
        );
    }

    #[test]
    fn born_rule_measurement_valid() {
        let rho = DensityMatrix5::pure_state(2); // |J>
                                                 // Measuring a pure |J> state should always return 2
        for seed in 0..100 {
            assert_eq!(rho.measure(seed), 2);
        }
    }

    #[test]
    fn quantum_locomotion_full_tick_cycle() {
        let mut state = QuantumLocomotionState::new(32);
        // Simulate 120 ticks with moving input and mild decoherence (dt-scaled)
        for _ in 0..120 {
            state.tick(1.0 / 60.0, &[0.0, 2.0, 0.0, 0.0, 0.0], 0.02);
        }
        // After sustained Moving bias, Moving population should be significant
        let pops = state.rho.populations();
        assert!(pops[1] > 0.1, "Moving population should be significant: {:?}", pops);
        // Trace should remain close to 1.0
        let trace: f32 = pops.iter().sum();
        assert!((trace - 1.0).abs() < 0.2, "Trace should be near 1.0: {}", trace);
    }
}