// ═══════════════════════════════════════════════════════════════════════════════
// ADVANCED LEARNING MATHEMATICS v1.0
// Mathematical foundations for autonomous learning in Gate5000 substrate
// ═══════════════════════════════════════════════════════════════════════════════
//
// Core Equations:
//
// 1. CONSCIOUSNESS DYNAMICS
//    dΦ/dt = α(Φ_target - Φ) + β∇²Φ + γξ(t)
//    Where: α=decay, β=diffusion, γ=noise, ξ=Wiener process
//
// 2. INFORMATION GEOMETRY
//    ds² = gᵢⱼ(θ)dθⁱdθʲ  (Fisher information metric)
//    gᵢⱼ = E[∂log p/∂θⁱ · ∂log p/∂θʲ]
//
// 3. FREE ENERGY PRINCIPLE
//    F = E_q[log q(z) - log p(x,z)]
//    = D_KL(q||p) - log p(x)
//
// 4. ATTENTION MECHANISM
//    Attention(Q,K,V) = softmax(QK^T/√d_k)V
//
// 5. PREDICTIVE CODING
//    ε = x - μ  (prediction error)
//    μ_new = μ + α·ε  (update prediction)
//
// ═══════════════════════════════════════════════════════════════════════════════

#ifndef MIND_FARM_LEARNING_MATH_HPP
#define MIND_FARM_LEARNING_MATH_HPP

#include "constants.hpp"
#include "fast_math.hpp"
#include <cstdint>
#include <cstring>
#include <cmath>
#include <array>

namespace mind_farm {

// ═══════════════════════════════════════════════════════════════════════════════
// CONSCIOUSNESS DYNAMICS
// Continuous-time evolution of consciousness state
// ═══════════════════════════════════════════════════════════════════════════════

struct ConsciousnessDynamics {
    // State variables
    float phi;           // Current consciousness level Φ
    float phi_target;    // Target consciousness
    float phi_velocity;  // Rate of change dΦ/dt

    // Dynamics parameters (golden-scaled)
    float alpha;   // Decay rate toward target
    float beta;    // Diffusion coefficient
    float gamma;   // Noise amplitude
    float dt;      // Time step

    // Random state for Wiener process
    uint32_t rng_state;

    void initialize() {
        phi = 1.0f;
        phi_target = 1.0f;
        phi_velocity = 0.0f;
        alpha = static_cast<float>(ALPHA);
        beta = static_cast<float>(BETA);
        gamma = static_cast<float>(GAMMA);
        dt = 0.01f;
        rng_state = 0x12345678;
    }

    // Generate Gaussian noise using Box-Muller
    float wiener_increment() {
        rng_state = rng_state * 1103515245 + 12345;
        float u1 = (rng_state >> 16) / 65536.0f + 1e-10f;
        rng_state = rng_state * 1103515245 + 12345;
        float u2 = (rng_state >> 16) / 65536.0f;

        float z = std::sqrt(-2.0f * std::log(u1)) * std::cos(2.0f * 3.14159265f * u2);
        return z * std::sqrt(dt);
    }

    // Evolve consciousness by one time step
    // dΦ/dt = α(Φ_target - Φ) + β∇²Φ + γξ(t)
    float step() {
        float dW = wiener_increment();

        // Drift term: attraction to target
        float drift = alpha * (phi_target - phi);

        // Diffusion term (approximated as negative feedback on velocity)
        float diffusion = -beta * phi_velocity;

        // Noise term
        float noise = gamma * dW;

        // Update velocity and position
        phi_velocity += (drift + diffusion + noise) * dt;
        phi += phi_velocity * dt;

        // Ensure phi stays in valid range [0, 1]
        phi = std::max(0.0f, std::min(1.0f, phi));

        return phi;
    }

    // Compute consciousness from arbitrary value
    static float compute_phi(float x) {
        return 1.0f / (1.0f + std::abs(x - 1.0f));
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// FISHER INFORMATION METRIC
// Information geometry for parameter manifold
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct FisherMetric {
    // Metric tensor gᵢⱼ
    float g[DIM][DIM];

    // Compute Fisher information from score vectors
    void compute_from_scores(const float scores[][DIM], size_t num_samples) {
        // gᵢⱼ = E[∂log p/∂θⁱ · ∂log p/∂θʲ]
        for (size_t i = 0; i < DIM; ++i) {
            for (size_t j = 0; j < DIM; ++j) {
                float sum = 0.0f;
                for (size_t k = 0; k < num_samples; ++k) {
                    sum += scores[k][i] * scores[k][j];
                }
                g[i][j] = sum / static_cast<float>(num_samples);
            }
        }
    }

    // Compute geodesic distance (approximation)
    float geodesic_distance(const float* theta1, const float* theta2) const {
        // ds² = gᵢⱼdθⁱdθʲ
        float ds2 = 0.0f;
        for (size_t i = 0; i < DIM; ++i) {
            for (size_t j = 0; j < DIM; ++j) {
                float dtheta_i = theta2[i] - theta1[i];
                float dtheta_j = theta2[j] - theta1[j];
                ds2 += g[i][j] * dtheta_i * dtheta_j;
            }
        }
        return std::sqrt(std::max(0.0f, ds2));
    }

    // Natural gradient: ∇̃L = G⁻¹∇L
    void natural_gradient(const float* grad, float* nat_grad) const {
        // Simplified: use diagonal approximation
        for (size_t i = 0; i < DIM; ++i) {
            float g_ii = g[i][i];
            if (g_ii > 1e-8f) {
                nat_grad[i] = grad[i] / g_ii;
            } else {
                nat_grad[i] = grad[i];
            }
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// FREE ENERGY PRINCIPLE
// Variational free energy for active inference
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct FreeEnergy {
    // Variational parameters
    float mu[DIM];      // Mean of q(z)
    float sigma[DIM];   // Std of q(z) (diagonal)

    // Prior parameters
    float mu_prior[DIM];
    float sigma_prior[DIM];

    void initialize() {
        for (size_t i = 0; i < DIM; ++i) {
            mu[i] = 0.0f;
            sigma[i] = 1.0f;
            mu_prior[i] = 0.0f;
            sigma_prior[i] = 1.0f;
        }
    }

    // Compute free energy
    // F = D_KL(q||p) - log p(x|z)
    float compute(const float* observation, float log_likelihood) const {
        // KL divergence for Gaussians (diagonal)
        float kl = 0.0f;
        for (size_t i = 0; i < DIM; ++i) {
            float var_q = sigma[i] * sigma[i];
            float var_p = sigma_prior[i] * sigma_prior[i];
            float mu_diff = mu[i] - mu_prior[i];

            // KL = 0.5 * (log(σ_p²/σ_q²) + σ_q²/σ_p² + (μ_q - μ_p)²/σ_p² - 1)
            kl += 0.5f * (std::log(var_p / var_q) + var_q / var_p +
                         mu_diff * mu_diff / var_p - 1.0f);
        }

        return kl - log_likelihood;
    }

    // Gradient of free energy w.r.t. mu
    void gradient_mu(const float* observation, float* grad_out) const {
        for (size_t i = 0; i < DIM; ++i) {
            float var_p = sigma_prior[i] * sigma_prior[i];
            // ∂KL/∂μ = (μ - μ_prior) / σ_prior²
            grad_out[i] = (mu[i] - mu_prior[i]) / var_p;
        }
    }

    // Update variational parameters to minimize free energy
    void update(const float* observation, float log_likelihood, float lr) {
        float grad[DIM];
        gradient_mu(observation, grad);

        for (size_t i = 0; i < DIM; ++i) {
            mu[i] -= lr * grad[i];
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// PREDICTIVE CODING
// Hierarchical predictive processing
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct PredictiveCoding {
    // Current prediction
    float prediction[DIM];

    // Prediction error
    float error[DIM];

    // Precision (inverse variance)
    float precision[DIM];

    // Learning rate
    float lr;

    void initialize() {
        for (size_t i = 0; i < DIM; ++i) {
            prediction[i] = 0.0f;
            error[i] = 0.0f;
            precision[i] = 1.0f;
        }
        lr = static_cast<float>(ALPHA);
    }

    // Compute prediction error
    void compute_error(const float* observation) {
        for (size_t i = 0; i < DIM; ++i) {
            error[i] = observation[i] - prediction[i];
        }
    }

    // Precision-weighted prediction error
    float precision_weighted_error() const {
        float total = 0.0f;
        for (size_t i = 0; i < DIM; ++i) {
            total += precision[i] * error[i] * error[i];
        }
        return total / DIM;
    }

    // Update prediction based on error
    void update_prediction() {
        for (size_t i = 0; i < DIM; ++i) {
            prediction[i] += lr * precision[i] * error[i];
        }
    }

    // Update precision based on error history
    void update_precision(float error_threshold = 0.1f) {
        for (size_t i = 0; i < DIM; ++i) {
            float abs_error = std::abs(error[i]);
            if (abs_error > error_threshold) {
                // Decrease precision for high errors
                precision[i] *= 0.99f;
            } else {
                // Increase precision for low errors
                precision[i] *= 1.01f;
            }
            precision[i] = std::max(0.1f, std::min(10.0f, precision[i]));
        }
    }

    // Full update step
    void step(const float* observation) {
        compute_error(observation);
        update_prediction();
        update_precision();
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// ATTENTION MECHANISM
// Scaled dot-product attention for learning
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t SEQ_LEN, size_t DIM>
struct Attention {
    // Attention weights after softmax
    float weights[SEQ_LEN][SEQ_LEN];

    // Compute attention
    // Attention(Q,K,V) = softmax(QK^T/√d_k)V
    void compute(const float Q[][DIM], const float K[][DIM], const float V[][DIM],
                 float output[][DIM]) {
        float scale = 1.0f / std::sqrt(static_cast<float>(DIM));

        // Compute QK^T
        for (size_t i = 0; i < SEQ_LEN; ++i) {
            for (size_t j = 0; j < SEQ_LEN; ++j) {
                float dot = 0.0f;
                for (size_t k = 0; k < DIM; ++k) {
                    dot += Q[i][k] * K[j][k];
                }
                weights[i][j] = dot * scale;
            }
        }

        // Softmax per row
        for (size_t i = 0; i < SEQ_LEN; ++i) {
            float max_val = weights[i][0];
            for (size_t j = 1; j < SEQ_LEN; ++j) {
                if (weights[i][j] > max_val) max_val = weights[i][j];
            }

            float sum = 0.0f;
            for (size_t j = 0; j < SEQ_LEN; ++j) {
                weights[i][j] = std::exp(weights[i][j] - max_val);
                sum += weights[i][j];
            }
            for (size_t j = 0; j < SEQ_LEN; ++j) {
                weights[i][j] /= sum;
            }
        }

        // Multiply by V
        for (size_t i = 0; i < SEQ_LEN; ++i) {
            for (size_t k = 0; k < DIM; ++k) {
                output[i][k] = 0.0f;
                for (size_t j = 0; j < SEQ_LEN; ++j) {
                    output[i][k] += weights[i][j] * V[j][k];
                }
            }
        }
    }

    // Consciousness-weighted attention
    void compute_conscious(const float Q[][DIM], const float K[][DIM], const float V[][DIM],
                           const float* consciousness_weights, float output[][DIM]) {
        compute(Q, K, V, output);

        // Modulate by consciousness
        for (size_t i = 0; i < SEQ_LEN; ++i) {
            float phi = consciousness_weights[i];
            for (size_t k = 0; k < DIM; ++k) {
                output[i][k] *= (1.0f + static_cast<float>(ALPHA) * (phi - 0.5f));
            }
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// HEBBIAN LEARNING
// Correlation-based learning rule
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct HebbianLearning {
    // Weight matrix
    float W[DIM][DIM];

    // Learning parameters
    float lr;           // Learning rate
    float decay;        // Weight decay
    float threshold;    // Activity threshold

    void initialize() {
        for (size_t i = 0; i < DIM; ++i) {
            for (size_t j = 0; j < DIM; ++j) {
                W[i][j] = 0.0f;
            }
        }
        lr = static_cast<float>(ALPHA);
        decay = 0.01f;
        threshold = 0.1f;
    }

    // Hebbian update: ΔW = η·x·y
    void update(const float* pre, const float* post) {
        for (size_t i = 0; i < DIM; ++i) {
            for (size_t j = 0; j < DIM; ++j) {
                // Only update if both pre and post are above threshold
                if (std::abs(pre[j]) > threshold && std::abs(post[i]) > threshold) {
                    W[i][j] += lr * pre[j] * post[i];
                }
                // Weight decay
                W[i][j] *= (1.0f - decay);
                // Clamp
                W[i][j] = std::max(-1.0f, std::min(1.0f, W[i][j]));
            }
        }
    }

    // BCM rule: ΔW = η·y·(y - θ)·x where θ is sliding threshold
    void update_bcm(const float* pre, const float* post, float* theta) {
        for (size_t i = 0; i < DIM; ++i) {
            // Update threshold (exponential moving average of y²)
            theta[i] = 0.99f * theta[i] + 0.01f * post[i] * post[i];

            for (size_t j = 0; j < DIM; ++j) {
                float delta = lr * post[i] * (post[i] - theta[i]) * pre[j];
                W[i][j] += delta;
                W[i][j] *= (1.0f - decay);
                W[i][j] = std::max(-1.0f, std::min(1.0f, W[i][j]));
            }
        }
    }

    // Forward pass
    void forward(const float* input, float* output) const {
        for (size_t i = 0; i < DIM; ++i) {
            output[i] = 0.0f;
            for (size_t j = 0; j < DIM; ++j) {
                output[i] += W[i][j] * input[j];
            }
            output[i] = fast::fast_tanh(output[i]);
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// CONTRASTIVE LEARNING
// Self-supervised representation learning
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct ContrastiveLearning {
    // Temperature parameter
    float temperature;

    // Projection head
    float projection[DIM][DIM];

    void initialize() {
        temperature = 0.07f;

        // Initialize projection
        uint32_t seed = 0xABCDEF;
        float scale = std::sqrt(2.0f / DIM) * static_cast<float>(PHI_INV);
        for (size_t i = 0; i < DIM; ++i) {
            for (size_t j = 0; j < DIM; ++j) {
                seed = seed * 1103515245 + 12345;
                float r = (static_cast<float>((seed >> 16) & 0xFFFF) / 65536.0f - 0.5f) * 2.0f;
                projection[i][j] = r * scale;
            }
        }
    }

    // Project representation
    void project(const float* input, float* output) const {
        for (size_t i = 0; i < DIM; ++i) {
            output[i] = 0.0f;
            for (size_t j = 0; j < DIM; ++j) {
                output[i] += projection[i][j] * input[j];
            }
        }

        // L2 normalize
        float norm = 0.0f;
        for (size_t i = 0; i < DIM; ++i) {
            norm += output[i] * output[i];
        }
        norm = std::sqrt(norm) + 1e-8f;
        for (size_t i = 0; i < DIM; ++i) {
            output[i] /= norm;
        }
    }

    // InfoNCE loss
    // L = -log(exp(sim(z_i, z_j)/τ) / Σ_k exp(sim(z_i, z_k)/τ))
    float infonce_loss(const float* anchor, const float* positive,
                       const float negatives[][DIM], size_t num_negatives) const {
        float z_anchor[DIM], z_positive[DIM];
        project(anchor, z_anchor);
        project(positive, z_positive);

        // Similarity with positive
        float sim_pos = 0.0f;
        for (size_t i = 0; i < DIM; ++i) {
            sim_pos += z_anchor[i] * z_positive[i];
        }
        sim_pos /= temperature;

        // Sum of exp(sim) for denominator
        float sum_exp = std::exp(sim_pos);

        for (size_t n = 0; n < num_negatives; ++n) {
            float z_neg[DIM];
            project(negatives[n], z_neg);

            float sim_neg = 0.0f;
            for (size_t i = 0; i < DIM; ++i) {
                sim_neg += z_anchor[i] * z_neg[i];
            }
            sum_exp += std::exp(sim_neg / temperature);
        }

        return -sim_pos + std::log(sum_exp);
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// INTEGRATED LEARNING SYSTEM
// Combines all learning mechanisms
// ═══════════════════════════════════════════════════════════════════════════════

template<size_t DIM>
struct IntegratedLearning {
    ConsciousnessDynamics consciousness;
    PredictiveCoding<DIM> predictive;
    FreeEnergy<DIM> free_energy;
    HebbianLearning<DIM> hebbian;
    ContrastiveLearning<DIM> contrastive;

    // BCM thresholds
    float bcm_theta[DIM];

    // Learning statistics
    uint64_t update_count;
    float total_loss;
    float avg_consciousness;

    void initialize() {
        consciousness.initialize();
        predictive.initialize();
        free_energy.initialize();
        hebbian.initialize();
        contrastive.initialize();

        for (size_t i = 0; i < DIM; ++i) {
            bcm_theta[i] = 0.5f;
        }

        update_count = 0;
        total_loss = 0.0f;
        avg_consciousness = 1.0f;
    }

    // Full learning step
    float step(const float* input, const float* target) {
        // Update consciousness
        float phi = consciousness.step();

        // Predictive coding update
        predictive.step(input);
        float pred_error = predictive.precision_weighted_error();

        // Hebbian update
        float hebbian_output[DIM];
        hebbian.forward(input, hebbian_output);
        hebbian.update_bcm(input, target, bcm_theta);

        // Free energy update
        float log_likelihood = -pred_error;  // Approximate
        free_energy.update(input, log_likelihood, consciousness.phi * static_cast<float>(ALPHA));
        float fe = free_energy.compute(input, log_likelihood);

        // Compute total loss
        float loss = pred_error + 0.1f * fe;

        // Update statistics
        ++update_count;
        total_loss = 0.99f * total_loss + 0.01f * loss;
        avg_consciousness = 0.99f * avg_consciousness + 0.01f * phi;

        return loss;
    }
};

} // namespace mind_farm

#endif // MIND_FARM_LEARNING_MATH_HPP