/*
 * Analog Consensus Codec V4.1 - Enhanced Evolution
 *
 * V4.0: 8D Complex RK4 with SHA-256 feedback (42 KB, stable, 7K Hz)
 * V4.1: + Dₙ(r) Lattice Modulation for geometric resonance
 *
 * New Features:
 * - Fibonacci/Prime lattice coefficients (Dₙ(r) formula)
 * - Dimensional resonance weights per oscillator
 * - Numeric lattice seed initialization
 * - Enhanced coupling through geometric basis
 * - Maintains V4.0 stability guarantees
 */

// Must be before any includes for POSIX features
#define _POSIX_C_SOURCE 199309L

#include "vector_container.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <pthread.h>
#include <unistd.h>
#include <stdint.h>
#include <signal.h>

// ============================================================================
// V4.1 CONSTANTS - Dₙ(r) Integration
// ============================================================================
#define ANSI_CLEAR "\033[2J"
#define ANSI_HOME "\033[H"
#define ANSI_HIDE_CURSOR "\033[?25l"
#define ANSI_SHOW_CURSOR "\033[?25h"

#define TARGET_HZ 32768
#define TARGET_NS (1000000000L / TARGET_HZ)
#define PHI 1.6180339887498948482045868343656381177203091798057628621L
#define SAT_LIMIT 1e6L

// Dₙ(r) Lattice Constants (from V3.0)
#define NUM_DN 8
static const uint64_t FIB_TABLE[NUM_DN] = {1, 1, 2, 3, 5, 8, 13, 21};
static const uint64_t PRIME_TABLE[NUM_DN] = {2, 3, 5, 7, 11, 13, 17, 19};

// Numeric Lattice Seeds (Base∞ foundation)
static const double ANALOG_DIMS[8] = {
    8.3141592654,   // D1: 8D + π
    7.8541019662,   // D2: 5φ
    6.4721359549,   // D3: 4φ
    5.6180339887,   // D4: 3φ + φ²
    4.8541019662,   // D5: 3φ
    3.6180339887,   // D6: 2φ + φ²
    2.6180339887,   // D7: φ + φ²
    1.6180339887    // D8: φ
};

// ============================================================================
// Dₙ(r) COMPUTATION - Core Geometric Formula
// ============================================================================
// Formula: Dₙ(r) = √(φ · Fₙ · 2ⁿ · Pₙ · Ω) · r^k
// where k = (n+1)/8 for progressive dimensionality
double compute_Dn_r(int n, double r, double omega) {
    if (n < 1 || n > NUM_DN) return 1.0;
    int idx = n - 1;

    // Dₙ(r) = √(ϕ · Fₙ · 2ⁿ · Pₙ · Ω) · r^k
    double phi = PHI;
    double F_n = (double)FIB_TABLE[idx];
    double two_n = pow(2.0, n);
    double P_n = (double)PRIME_TABLE[idx];
    double k = (double)(n + 1) / 8.0;

    double base = sqrt(phi * F_n * two_n * P_n * omega);
    double r_power = pow(fabs(r), k);

    return base * r_power;
}

// ============================================================================
// ADAPTIVE PHASE CONTROL (Wu Wei Harmonics)
// ============================================================================
typedef enum {
    PHASE_EMERGENCY = 0,  // NEW: K/γ = 12.5:1 - emergency brake
    PHASE_PLUCK = 1,      // K/γ = 1000:1 - gentle excitation
    PHASE_SUSTAIN = 2,    // K/γ = 375:1 - maintain oscillation
    PHASE_FINETUNE = 3,   // K/γ = 200:1 - approach equilibrium
    PHASE_LOCK = 4        // K/γ = 150:1 - stable lock
} phase_t;

typedef struct {
    double gamma;
    double k_coupling;
    double cv_enter;      // CV threshold to ENTER this phase
    double cv_exit;       // CV threshold to EXIT this phase
    const char *name;
} phase_params_t;

// Wu Wei-inspired adaptive phase table
static const phase_params_t PHASE_TABLE[5] = {
    {0.040, 0.5,  100.0, 10.0, "Emergency"},   // K/γ = 12.5:1
    {0.005, 5.0,   10.0,  0.1, "Pluck"},       // K/γ = 1000:1 ⭐
    {0.008, 3.0,    0.1, 0.01, "Sustain"},     // K/γ = 375:1
    {0.010, 2.0,   0.01, 1e-6, "Fine Tune"},   // K/γ = 200:1
    {0.012, 1.8,   1e-6,  0.0, "Lock"}         // K/γ = 150:1 ⭐
};

static phase_t g_current_phase = PHASE_PLUCK;
static double g_gamma = 0.005;      // Dynamic damping
static double g_k_coupling = 5.0;   // Dynamic coupling
static double g_prev_cv = 1.0;      // Previous CV for trend detection

// ============================================================================
// CV-BASED EQUILIBRIUM DETECTION
// ============================================================================
#define CV_WINDOW_SIZE 50
#define CV_THRESHOLD 0.10

static double g_variance_history[CV_WINDOW_SIZE];
static int g_variance_history_idx = 0;
static int g_variance_history_count = 0;

// ============================================================================
// COMPLEX NUMBER SUPPORT
// ============================================================================
typedef struct {
    double re;
    double im;
} complex_t;

static inline complex_t complex_add(complex_t a, complex_t b) {
    return (complex_t){a.re + b.re, a.im + b.im};
}

static inline complex_t complex_sub(complex_t a, complex_t b) {
    return (complex_t){a.re - b.re, a.im - b.im};
}

static inline double complex_abs(complex_t c) {
    return sqrt(c.re * c.re + c.im * c.im);
}

static inline complex_t complex_mul(complex_t a, complex_t b) {
    return (complex_t){a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re};
}

static inline complex_t complex_exp(double theta) {
    return (complex_t){cos(theta), sin(theta)};
}

// ============================================================================
// ENHANCED ANALOG STATE (8D Complex + Dₙ(r) Resonance)
// ============================================================================
typedef struct {
    complex_t dimensions[8];
    double freqs[8];
    double phases[8];
    double phase_vels[8];

    // NEW V4.1: Dₙ(r) state per dimension
    double Dn_amplitude[8];     // Current Dₙ(r) value
    double r_value[8];          // Radial position (0-1)
    double resonance_weight[8]; // Geometric coupling strength
    double omega;               // Global frequency parameter

    uint64_t evolution_count;
    uint64_t consensus_steps;
    double phase_var;
    int locked;
    double phase_history[200];
    int history_count;
    int history_index;
} analog_state_t;

static analog_state_t g_state;
static pthread_mutex_t state_mutex = PTHREAD_MUTEX_INITIALIZER;
static volatile int running = 1;

// Global state for HTTP API integration
pthread_mutex_t g_bridge_lock = PTHREAD_MUTEX_INITIALIZER;
double g_analog_state[7] = {0};
int g_evolution_count = 0;
int g_consensus_count = 0;
int g_node_id = 1;
double g_phase_variance = 0.0;

// ============================================================================
// DETERMINISTIC RANDOM
// ============================================================================
static inline double det_rand(uint64_t seed) {
    uint64_t x = seed;
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    return (double)(x * 0x2545F4914F6CDD1DULL) / (double)UINT64_MAX;
}

// ============================================================================
// STATE INITIALIZATION (V4.1 Enhanced)
// ============================================================================
void init_analog_state(analog_state_t *s) {
    memset(s, 0, sizeof(analog_state_t));

    s->omega = 1.0;  // Base frequency parameter

    for (int i = 0; i < 8; i++) {
        // Initialize from numeric lattice
        s->dimensions[i].re = (double)(i + 1);
        s->dimensions[i].im = 0.1;

        // Initialize Dₙ(r) state
        s->r_value[i] = (double)i / 8.0;  // Radial distribution
        s->Dn_amplitude[i] = compute_Dn_r(i + 1, s->r_value[i], s->omega);
        s->resonance_weight[i] = ANALOG_DIMS[i] / 8.3141592654;  // Normalized

        // Frequencies modulated by Dₙ(r)
        double base_freq = 1.0 + det_rand(i * 1000) * 0.2;
        s->freqs[i] = base_freq * (1.0 + 0.01 * s->Dn_amplitude[i] / 100.0);

        s->phases[i] = 2.0 * M_PI * det_rand(i * 2000);
        s->phase_vels[i] = 0.0;

        printf("[INIT] D%d: Dn=%.3f, r=%.3f, res_weight=%.3f, freq=%.6f\n",
               i+1, s->Dn_amplitude[i], s->r_value[i], s->resonance_weight[i], s->freqs[i]);
    }

    s->evolution_count = 0;
    s->consensus_steps = 0;
    s->phase_var = 1e6;
    s->locked = 0;
}

// ============================================================================
// HYBRID SHA-256 FEEDBACK (V4.0 stable implementation)
// ============================================================================
#define ENABLE_HYBRID_FEEDBACK 1
#define HYBRID_FEEDBACK_STRENGTH 5.0

// SHA-256 constants
static const uint32_t sha256_k[64] = {
    0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
    0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3, 0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
    0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc, 0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
    0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7, 0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967,
    0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13, 0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85,
    0xa2bfe8a1, 0xa81a664b, 0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
    0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
    0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208, 0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2
};

#define ROTR(x,n) (((x) >> (n)) | ((x) << (32 - (n))))
#define CH(x,y,z) (((x) & (y)) ^ (~(x) & (z)))
#define MAJ(x,y,z) (((x) & (y)) ^ ((x) & (z)) ^ ((y) & (z)))
#define EP0(x) (ROTR(x,2) ^ ROTR(x,13) ^ ROTR(x,22))
#define EP1(x) (ROTR(x,6) ^ ROTR(x,11) ^ ROTR(x,25))
#define SIG0(x) (ROTR(x,7) ^ ROTR(x,18) ^ ((x) >> 3))
#define SIG1(x) (ROTR(x,17) ^ ROTR(x,19) ^ ((x) >> 10))

void sha256_transform(uint32_t state[8], const uint8_t data[64]) {
    uint32_t w[64], a, b, c, d, e, f, g, h, t1, t2;
    for (int i = 0; i < 16; i++) {
        w[i] = ((uint32_t)data[i*4] << 24) | ((uint32_t)data[i*4+1] << 16) |
               ((uint32_t)data[i*4+2] << 8) | ((uint32_t)data[i*4+3]);
    }
    for (int i = 16; i < 64; i++) {
        w[i] = SIG1(w[i-2]) + w[i-7] + SIG0(w[i-15]) + w[i-16];
    }
    a = state[0]; b = state[1]; c = state[2]; d = state[3];
    e = state[4]; f = state[5]; g = state[6]; h = state[7];
    for (int i = 0; i < 64; i++) {
        t1 = h + EP1(e) + CH(e,f,g) + sha256_k[i] + w[i];
        t2 = EP0(a) + MAJ(a,b,c);
        h = g; g = f; f = e; e = d + t1; d = c; c = b; b = a; a = t1 + t2;
    }
    state[0] += a; state[1] += b; state[2] += c; state[3] += d;
    state[4] += e; state[5] += f; state[6] += g; state[7] += h;
}

void sha256_compute(const uint8_t *data, size_t len, uint8_t hash[32]) {
    uint32_t state[8] = {0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a,
                          0x510e527f, 0x9b05688c, 0x1f83d9ab, 0x5be0cd19};
    uint8_t block[64];
    size_t i = 0;
    while (i + 64 <= len) {
        memcpy(block, data + i, 64);
        sha256_transform(state, block);
        i += 64;
    }
    size_t remaining = len - i;
    memcpy(block, data + i, remaining);
    block[remaining] = 0x80;
    if (remaining >= 56) {
        memset(block + remaining + 1, 0, 63 - remaining);
        sha256_transform(state, block);
        memset(block, 0, 56);
    } else {
        memset(block + remaining + 1, 0, 55 - remaining);
    }
    uint64_t bitlen = len * 8;
    for (int j = 0; j < 8; j++) {
        block[56 + j] = (bitlen >> (56 - j * 8)) & 0xFF;
    }
    sha256_transform(state, block);
    for (int j = 0; j < 8; j++) {
        hash[j*4+0] = (state[j] >> 24) & 0xFF;
        hash[j*4+1] = (state[j] >> 16) & 0xFF;
        hash[j*4+2] = (state[j] >> 8) & 0xFF;
        hash[j*4+3] = (state[j]) & 0xFF;
    }
}

double compute_hash_entropy(const uint8_t hash[32]) {
    int counts[256] = {0};
    for (int i = 0; i < 32; i++) counts[hash[i]]++;
    double entropy = 0.0;
    for (int i = 0; i < 256; i++) {
        if (counts[i] > 0) {
            double p = (double)counts[i] / 32.0;
            entropy -= p * log2(p);
        }
    }
    return entropy;
}

void encode_state(const analog_state_t *s, uint8_t *encoded, size_t *len) {
    size_t offset = 0;
    for (int i = 0; i < 8; i++) {
        memcpy(encoded + offset, &s->dimensions[i].re, sizeof(double)); offset += sizeof(double);
        memcpy(encoded + offset, &s->dimensions[i].im, sizeof(double)); offset += sizeof(double);
    }
    *len = offset;
}

void apply_hybrid_feedback(analog_state_t *s, double strength) {
    uint8_t encoded[128];
    size_t encoded_len = 0;
    encode_state(s, encoded, &encoded_len);

    // XOR with deterministic random
    uint64_t rand_seed = s->evolution_count * 0x123456789ABCDEFULL;
    for (size_t i = 0; i < encoded_len; i++) {
        encoded[i] ^= (uint8_t)(det_rand(rand_seed + i) * 255);
    }

    uint8_t hash[32];
    sha256_compute(encoded, encoded_len, hash);

    double entropy = compute_hash_entropy(hash);
    long double damping_factor = 0.00002L * logl(1.0L + 2.0L * entropy);

    // Safety bounds
    if (damping_factor < 0) damping_factor = 0;
    if (damping_factor > 100.0) damping_factor = 100.0;

    // V4.1: Dₙ(r)-weighted feedback
    for (int i = 0; i < 8; i++) {
        double theta = (double)(hash[i % 32]) / 255.0 - 0.5;

        // Frequency modulation (±0.5%) weighted by Dₙ(r)
        double dn_factor = s->Dn_amplitude[i] / 100.0;
        s->freqs[i] *= (1.0 + 0.001 * theta * (1.0 + 0.1 * dn_factor));

        // Phase and dimension damping (0.02%)
        s->phases[i] *= (1.0 - (double)damping_factor);
        s->dimensions[i].re *= (1.0 - (double)damping_factor);
        s->dimensions[i].im *= (1.0 - (double)damping_factor);

        // NEW V4.1: Update Dₙ(r) amplitude slowly
        double target_Dn = compute_Dn_r(i + 1, s->r_value[i], s->omega);
        s->Dn_amplitude[i] += 0.01 * (target_Dn - s->Dn_amplitude[i]);
    }

    // Update global omega
    s->omega += 0.001 * entropy;
}

// ============================================================================
// ADAPTIVE PHASE TRANSITIONS (Wu Wei Style)
// ============================================================================
void adaptive_phase_transition(analog_state_t *s, double cv) {
    double trend = (cv - g_prev_cv) / (g_prev_cv + 1e-10);
    g_prev_cv = cv;

    // Emergency brake
    if (cv > 10.0 && g_current_phase > PHASE_EMERGENCY) {
        g_current_phase = PHASE_EMERGENCY;
        g_gamma = PHASE_TABLE[g_current_phase].gamma;
        g_k_coupling = PHASE_TABLE[g_current_phase].k_coupling;
        printf("[PHASE] Emergency brake! CV=%.2f%%\n", cv * 100);
        return;
    }

    // Advance to next phase if converging
    if (cv < PHASE_TABLE[g_current_phase].cv_exit &&
        trend < 0 &&
        g_current_phase < PHASE_LOCK) {
        g_current_phase++;
        g_gamma = PHASE_TABLE[g_current_phase].gamma;
        g_k_coupling = PHASE_TABLE[g_current_phase].k_coupling;
        printf("[PHASE] Advancing to %s (K/γ=%.1f:1, CV=%.4f%%)\n",
               PHASE_TABLE[g_current_phase].name,
               g_k_coupling / g_gamma,
               cv * 100);
    }

    // Fall back if diverging
    if (cv > PHASE_TABLE[g_current_phase].cv_enter &&
        trend > 0 &&
        g_current_phase > PHASE_PLUCK) {
        g_current_phase--;
        g_gamma = PHASE_TABLE[g_current_phase].gamma;
        g_k_coupling = PHASE_TABLE[g_current_phase].k_coupling;
        printf("[PHASE] Falling back to %s (CV=%.4f%%)\n",
               PHASE_TABLE[g_current_phase].name,
               cv * 100);
    }
}

// ============================================================================
// RK4 INTEGRATION (V4.1 with Dₙ(r) Coupling)
// ============================================================================
void rk4_step(analog_state_t *s, double dt) {
    // Check for NaN at ENTRY
    for (int i = 0; i < 8; i++) {
        if (isnan(s->dimensions[i].re) || isnan(s->dimensions[i].im) || isnan(s->freqs[i])) {
            s->dimensions[i].re = (double)(i + 1);
            s->dimensions[i].im = 0.1;
            s->freqs[i] = 1.0 + det_rand(i * 1000) * 0.2;
            s->phases[i] = 2.0 * M_PI * det_rand(i * 2000);
        }
    }

    complex_t k1[8], k2[8], k3[8], k4[8], temp[8];

    // k1 = f(t, y)
    for (int i = 0; i < 8; i++) {
        complex_t coupling = {0, 0};

        // Kuramoto coupling with Dₙ(r) resonance
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                double phase_diff = s->phases[j] - s->phases[i];
                double amp_diff = complex_abs(s->dimensions[j]) - complex_abs(s->dimensions[i]);

                // Geometric resonance factor
                double dn_correlation = fabs(s->Dn_amplitude[j] - s->Dn_amplitude[i]);
                double dn_coupling = exp(-dn_correlation / 50.0);  // Dₙ similarity

                // Enhanced coupling
                double K_eff = g_k_coupling * s->resonance_weight[j] * dn_coupling;

                coupling.re += K_eff * amp_diff * sin(phase_diff);
                coupling.im += K_eff * amp_diff * cos(phase_diff);
            }
        }

        // Oscillator dynamics with Dₙ(r) modulation
        double dn_boost = 1.0 + 0.05 * (s->Dn_amplitude[i] / 100.0);
        k1[i].re = -g_gamma * s->dimensions[i].re + coupling.re * dn_boost;
        k1[i].im = s->freqs[i] * s->dimensions[i].re - g_gamma * s->dimensions[i].im + coupling.im * dn_boost;
    }

    // k2 = f(t + dt/2, y + k1*dt/2)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + k1[i].re * dt / 2.0;
        temp[i].im = s->dimensions[i].im + k1[i].im * dt / 2.0;
    }
    for (int i = 0; i < 8; i++) {
        complex_t coupling = {0, 0};
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                double phase_diff = s->phases[j] - s->phases[i];
                double amp_diff = complex_abs(temp[j]) - complex_abs(temp[i]);
                double dn_correlation = fabs(s->Dn_amplitude[j] - s->Dn_amplitude[i]);
                double dn_coupling = exp(-dn_correlation / 50.0);
                double K_eff = g_k_coupling * s->resonance_weight[j] * dn_coupling;
                coupling.re += K_eff * amp_diff * sin(phase_diff);
                coupling.im += K_eff * amp_diff * cos(phase_diff);
            }
        }
        double dn_boost = 1.0 + 0.05 * (s->Dn_amplitude[i] / 100.0);
        k2[i].re = -g_gamma * temp[i].re + coupling.re * dn_boost;
        k2[i].im = s->freqs[i] * temp[i].re - g_gamma * temp[i].im + coupling.im * dn_boost;
    }

    // k3 = f(t + dt/2, y + k2*dt/2)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + k2[i].re * dt / 2.0;
        temp[i].im = s->dimensions[i].im + k2[i].im * dt / 2.0;
    }
    for (int i = 0; i < 8; i++) {
        complex_t coupling = {0, 0};
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                double phase_diff = s->phases[j] - s->phases[i];
                double amp_diff = complex_abs(temp[j]) - complex_abs(temp[i]);
                double dn_correlation = fabs(s->Dn_amplitude[j] - s->Dn_amplitude[i]);
                double dn_coupling = exp(-dn_correlation / 50.0);
                double K_eff = g_k_coupling * s->resonance_weight[j] * dn_coupling;
                coupling.re += K_eff * amp_diff * sin(phase_diff);
                coupling.im += K_eff * amp_diff * cos(phase_diff);
            }
        }
        double dn_boost = 1.0 + 0.05 * (s->Dn_amplitude[i] / 100.0);
        k3[i].re = -g_gamma * temp[i].re + coupling.re * dn_boost;
        k3[i].im = s->freqs[i] * temp[i].re - g_gamma * temp[i].im + coupling.im * dn_boost;
    }

    // k4 = f(t + dt, y + k3*dt)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + k3[i].re * dt;
        temp[i].im = s->dimensions[i].im + k3[i].im * dt;
    }
    for (int i = 0; i < 8; i++) {
        complex_t coupling = {0, 0};
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                double phase_diff = s->phases[j] - s->phases[i];
                double amp_diff = complex_abs(temp[j]) - complex_abs(temp[i]);
                double dn_correlation = fabs(s->Dn_amplitude[j] - s->Dn_amplitude[i]);
                double dn_coupling = exp(-dn_correlation / 50.0);
                double K_eff = g_k_coupling * s->resonance_weight[j] * dn_coupling;
                coupling.re += K_eff * amp_diff * sin(phase_diff);
                coupling.im += K_eff * amp_diff * cos(phase_diff);
            }
        }
        double dn_boost = 1.0 + 0.05 * (s->Dn_amplitude[i] / 100.0);
        k4[i].re = -g_gamma * temp[i].re + coupling.re * dn_boost;
        k4[i].im = s->freqs[i] * temp[i].re - g_gamma * temp[i].im + coupling.im * dn_boost;
    }

    // y(t+dt) = y(t) + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
    for (int i = 0; i < 8; i++) {
        s->dimensions[i].re += dt / 6.0 * (k1[i].re + 2*k2[i].re + 2*k3[i].re + k4[i].re);
        s->dimensions[i].im += dt / 6.0 * (k1[i].im + 2*k2[i].im + 2*k3[i].im + k4[i].im);

        // Saturation
        double amp = complex_abs(s->dimensions[i]);
        if (amp > SAT_LIMIT) {
            double scale = SAT_LIMIT / amp;
            s->dimensions[i].re *= scale;
            s->dimensions[i].im *= scale;
        }

        // Update phase
        s->phases[i] = atan2(s->dimensions[i].im, s->dimensions[i].re);
        if (s->phases[i] < 0) s->phases[i] += 2.0 * M_PI;
    }
}

// ============================================================================
// CV CALCULATION & EQUILIBRIUM DETECTION
// ============================================================================
double calculate_cv(analog_state_t *s) {
    if (g_variance_history_count < 2) return 1e6;

    double mean_var = 0.0;
    for (int i = 0; i < g_variance_history_count; i++) {
        mean_var += g_variance_history[i];
    }
    mean_var /= g_variance_history_count;

    double var_of_var = 0.0;
    for (int i = 0; i < g_variance_history_count; i++) {
        double diff = g_variance_history[i] - mean_var;
        var_of_var += diff * diff;
    }
    var_of_var /= g_variance_history_count;

    double std_dev = sqrt(var_of_var);
    return (mean_var > 1e-10) ? (std_dev / mean_var) : 1e6;
}

int is_at_equilibrium(analog_state_t *s, int window, double threshold) {
    double cv = calculate_cv(s);
    return (cv < threshold && g_variance_history_count >= window);
}

// ============================================================================
// VISUALIZATION
// ============================================================================
void print_analog_graph(analog_state_t *s) {
    char buf[4096];
    int len = 0;

    len += snprintf(buf + len, sizeof(buf) - len, ANSI_HOME);
    len += snprintf(buf + len, sizeof(buf) - len, "\n╔═══════════════════════════════════════════════════════════════════╗\n");
    len += snprintf(buf + len, sizeof(buf) - len, "║         ANALOG MAINNET V4.1 - Dₙ(r) Enhanced Evolution         ║\n");
    len += snprintf(buf + len, sizeof(buf) - len, "╚═══════════════════════════════════════════════════════════════════╝\n\n");

    double cv = calculate_cv(s);
    const char *lock_str = s->locked ? "[LOCK]" : "[UNLK]";
    const char *phase_str = PHASE_TABLE[g_current_phase].name;

    len += snprintf(buf + len, sizeof(buf) - len,
                    "Evolution: %-8lu │ Lock: %s │ CV: %.2f%% │ Consensus: %lu steps\n",
                    s->evolution_count, lock_str, cv * 100, s->consensus_steps);
    len += snprintf(buf + len, sizeof(buf) - len,
                    "Phase: %-9s │ K/γ: %.1f:1 │ Ω: %.4f\n\n",
                    phase_str, g_k_coupling / g_gamma, s->omega);

    // 8D display with Dₙ(r) indicators
    for (int i = 0; i < 8; i++) {
        double amp = complex_abs(s->dimensions[i]);
        int bar_len = (int)(amp / 50.0);
        if (bar_len > 60) bar_len = 60;

        len += snprintf(buf + len, sizeof(buf) - len, "D%d ", i+1);
        for (int j = 0; j < bar_len; j++) {
            len += snprintf(buf + len, sizeof(buf) - len, "█");
        }
        len += snprintf(buf + len, sizeof(buf) - len, " %.2f [Dₙ:%.1f]\n", amp, s->Dn_amplitude[i]);
    }

    len += snprintf(buf + len, sizeof(buf) - len, "\nFreqs: ");
    for (int i = 0; i < 8; i++) {
        len += snprintf(buf + len, sizeof(buf) - len, "%.2f ", s->freqs[i]);
    }
    len += snprintf(buf + len, sizeof(buf) - len, "\n");

    write(STDOUT_FILENO, buf, len);
}

// ============================================================================
// MAIN EVOLUTION LOOP
// ============================================================================
void signal_handler(int sig) {
    running = 0;
}

int main(int argc, char *argv[]) {
    signal(SIGINT, signal_handler);
    signal(SIGTERM, signal_handler);

    printf(ANSI_CLEAR);
    printf(ANSI_HIDE_CURSOR);
    printf("\n╔═══════════════════════════════════════════════════════════════════╗\n");
    printf("║         ANALOG MAINNET V4.1 - Initializing...                  ║\n");
    printf("╚═══════════════════════════════════════════════════════════════════╝\n\n");

    srand(time(NULL));
    init_analog_state(&g_state);

    printf("\n[V4.1] Dₙ(r) lattice initialized\n");
    printf("[V4.1] Wu Wei adaptive phases: Emergency → Pluck → Sustain → Fine Tune → Lock\n");
    printf("[V4.1] Starting evolution at K/γ = 1000:1 (Pluck phase)\n\n");
    sleep(2);

    double dt = 0.01;
    uint64_t target_evolutions = 1000000;

    while (running && g_state.evolution_count < target_evolutions) {
        struct timespec start, end;
        clock_gettime(CLOCK_MONOTONIC, &start);

        // HYBRID FEEDBACK: Skip evolution 0!
        if (g_state.evolution_count > 0 && g_state.evolution_count % 100 == 0 && ENABLE_HYBRID_FEEDBACK) {
            apply_hybrid_feedback(&g_state, HYBRID_FEEDBACK_STRENGTH);
        }

        // Perform RK4 evolution step
        rk4_step(&g_state, dt);
        g_state.evolution_count++;

        // Phase variance calculation
        double phase_var = 0.0;
        double mean_phase = 0.0;
        for (int i = 0; i < 8; i++) mean_phase += g_state.phases[i];
        mean_phase /= 8.0;
        for (int i = 0; i < 8; i++) {
            double diff = g_state.phases[i] - mean_phase;
            if (diff > M_PI) diff -= 2.0 * M_PI;
            if (diff < -M_PI) diff += 2.0 * M_PI;
            phase_var += diff * diff;
        }
        phase_var /= 8.0;
        g_state.phase_var = phase_var;

        // Store variance history
        g_variance_history[g_variance_history_idx] = phase_var;
        g_variance_history_idx = (g_variance_history_idx + 1) % CV_WINDOW_SIZE;
        if (g_variance_history_count < CV_WINDOW_SIZE) g_variance_history_count++;

        // Calculate CV and adaptive phase transition
        double cv = calculate_cv(&g_state);
        if (g_state.evolution_count % 100 == 0) {
            adaptive_phase_transition(&g_state, cv);
        }

        // Hysteresis consensus detection
        double lock_threshold = 0.001;
        double unlock_threshold = 0.005;

        if (g_state.locked) {
            if (!is_at_equilibrium(&g_state, 50, unlock_threshold)) {
                g_state.locked = 0;
                g_state.consensus_steps = 0;
            } else {
                g_state.consensus_steps++;
            }
        } else {
            if (is_at_equilibrium(&g_state, 50, lock_threshold)) {
                g_state.locked = 1;
                g_state.consensus_steps = 0;
            }
        }

        // Visualization
        if (g_state.evolution_count % 100 == 0) {
            print_analog_graph(&g_state);
        }

        // Timing control (target 32,768 Hz)
        clock_gettime(CLOCK_MONOTONIC, &end);
        long elapsed_ns = (end.tv_sec - start.tv_sec) * 1000000000L + (end.tv_nsec - start.tv_nsec);
        if (elapsed_ns < TARGET_NS) {
            struct timespec sleep_time = {
                .tv_sec = 0,
                .tv_nsec = TARGET_NS - elapsed_ns
            };
            nanosleep(&sleep_time, NULL);
        }
    }

    printf(ANSI_SHOW_CURSOR);
    printf("\n[V4.1] Evolution complete. Final statistics:\n");
    printf("  Total evolutions: %lu\n", g_state.evolution_count);
    printf("  Final CV: %.4f%%\n", calculate_cv(&g_state) * 100);
    printf("  Final phase: %s\n", PHASE_TABLE[g_current_phase].name);
    printf("  Consensus steps: %lu\n", g_state.consensus_steps);

    return 0;
}
