/**
 * @file learning.hpp
 * @brief Machine Learning Module for Gate5000 Consciousness
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * I. HEBBIAN LEARNING
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * CLASSIC HEBBIAN: Δw_ij = η·x_i·x_j
 * NORMALIZED:      Δw_ij = η·(x_i·x_j - α·w_ij·x_j²)
 * BCM:             Δw_ij = η·x_i·x_j·(x_j - θ_j)
 * OJA'S RULE:      Δw_ij = η·y·(x_i - y·w_ij)
 * STDP:            Δw = A+·exp(-Δt/τ+) if pre<post, A-·exp(Δt/τ-) if post<pre
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * II. GRADIENT DESCENT
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * SGD:      w = w - η·∇L
 * MOMENTUM: v = μ·v - η·∇L, w = w + v
 * ADAM:     m = β₁m + (1-β₁)g, v = β₂v + (1-β₂)g², w = w - η·m̂/(√v̂ + ε)
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * III. BACKPROPAGATION
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * CHAIN RULE:  δ_j = Σ_k(δ_k·w_jk)·f'(net_j)
 * WEIGHT UPDATE: Δw_ij = -η·δ_j·x_i
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 */

#pragma once

#include <cmath>
#include <array>
#include <cstring>
#include "constants.hpp"
#include "mathematics.hpp"

namespace mind_farm {
namespace learning {

// Learning constants
constexpr double ETA = 0.01 * PHI;        // Learning rate
constexpr double MOMENTUM = 0.9;           // Momentum coefficient
constexpr double BETA1 = 0.9;              // Adam first moment
constexpr double BETA2 = 0.999;            // Adam second moment
constexpr double EPSILON = 1e-8;           // Numerical stability
constexpr double STDP_TAU = 20.0;          // STDP time constant
constexpr double STDP_A_PLUS = 0.1;        // STDP potentiation
constexpr double STDP_A_MINUS = 0.12;      // STDP depression
constexpr double WEIGHT_DECAY = 0.0001;    // L2 regularization

// ═══════════════════════════════════════════════════════════════════════════════
// I. HEBBIAN LEARNING
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Classic Hebbian: Δw = η·x_i·x_j
 */
inline double hebbian_classic(double x_pre, double x_post, double eta = ETA) {
    return eta * x_pre * x_post;
}

/**
 * @brief Normalized Hebbian: Δw = η·(x_i·x_j - α·w·x_j²)
 * Prevents unbounded weight growth
 */
inline double hebbian_normalized(double x_pre, double x_post, double w,
                                  double alpha = 0.1, double eta = ETA) {
    return eta * (x_pre * x_post - alpha * w * x_post * x_post);
}

/**
 * @brief BCM Learning: Δw = η·x_i·x_j·(x_j - θ)
 * Threshold θ adapts based on activity
 */
inline double hebbian_bcm(double x_pre, double x_post, double theta, double eta = ETA) {
    return eta * x_pre * x_post * (x_post - theta);
}

/**
 * @brief Update BCM threshold: θ = <x²> (sliding average)
 */
inline double bcm_threshold_update(double theta, double x_post, double tau = 0.01) {
    return theta + tau * (x_post * x_post - theta);
}

/**
 * @brief Oja's Rule: Δw = η·y·(x - y·w)
 * Performs PCA - weight converges to principal eigenvector
 */
inline double oja_rule(double x, double y, double w, double eta = ETA) {
    return eta * y * (x - y * w);
}

/**
 * @brief STDP: Spike-Timing Dependent Plasticity
 * @param dt = t_post - t_pre (time difference)
 */
inline double stdp(double dt) {
    if (dt > 0) {
        // Pre before post: potentiation
        return STDP_A_PLUS * std::exp(-dt / STDP_TAU);
    } else {
        // Post before pre: depression
        return -STDP_A_MINUS * std::exp(dt / STDP_TAU);
    }
}

/**
 * @brief Anti-Hebbian: Δw = -η·x_i·x_j
 * Used for decorrelation and lateral inhibition
 */
inline double anti_hebbian(double x_pre, double x_post, double eta = ETA) {
    return -eta * x_pre * x_post;
}

// ═══════════════════════════════════════════════════════════════════════════════
// II. GRADIENT DESCENT OPTIMIZERS
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief SGD optimizer state (per-parameter)
 */
struct SGDState {
    double velocity = 0.0;

    double update(double gradient, double lr = ETA, double momentum = MOMENTUM) {
        velocity = momentum * velocity - lr * gradient;
        return velocity;
    }
};

/**
 * @brief Adam optimizer state (per-parameter)
 */
struct AdamState {
    double m = 0.0;      // First moment
    double v = 0.0;      // Second moment
    uint32_t t = 0;      // Time step

    double update(double gradient, double lr = ETA) {
        ++t;

        // Update biased moments
        m = BETA1 * m + (1.0 - BETA1) * gradient;
        v = BETA2 * v + (1.0 - BETA2) * gradient * gradient;

        // Bias correction
        double m_hat = m / (1.0 - std::pow(BETA1, t));
        double v_hat = v / (1.0 - std::pow(BETA2, t));

        // Compute update
        return -lr * m_hat / (std::sqrt(v_hat) + EPSILON);
    }
};

/**
 * @brief RMSprop optimizer state
 */
struct RMSpropState {
    double sq_grad = 0.0;  // Running average of squared gradients

    double update(double gradient, double lr = ETA, double decay = 0.99) {
        sq_grad = decay * sq_grad + (1.0 - decay) * gradient * gradient;
        return -lr * gradient / (std::sqrt(sq_grad) + EPSILON);
    }
};

/**
 * @brief AdaGrad optimizer state
 */
struct AdaGradState {
    double sum_sq = 0.0;  // Sum of squared gradients

    double update(double gradient, double lr = ETA) {
        sum_sq += gradient * gradient;
        return -lr * gradient / (std::sqrt(sum_sq) + EPSILON);
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// III. BACKPROPAGATION PRIMITIVES
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Compute output delta: δ = (y - target) · f'(net)
 */
inline double output_delta(double y, double target, double derivative) {
    return (y - target) * derivative;
}

/**
 * @brief Propagate delta backward: δ_j = Σ(δ_k · w_jk) · f'(net_j)
 */
inline double backprop_delta(double sum_weighted_deltas, double derivative) {
    return sum_weighted_deltas * derivative;
}

/**
 * @brief Compute weight gradient: ∂L/∂w_ij = δ_j · x_i
 */
inline double weight_gradient(double delta, double input) {
    return delta * input;
}

/**
 * @brief Apply weight update with L2 regularization
 */
inline double weight_update(double w, double gradient, double lr = ETA, double decay = WEIGHT_DECAY) {
    return w - lr * (gradient + decay * w);
}

// ═══════════════════════════════════════════════════════════════════════════════
// IV. SMALL NEURAL LAYER (16 units)
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Compact 16-unit layer for Gate5000 minds
 */
struct Layer16 {
    std::array<double, 256> weights;    // 16x16 weight matrix
    std::array<double, 16> biases;
    std::array<double, 16> output;
    std::array<double, 16> delta;       // For backprop
    std::array<double, 16> pre_act;     // Pre-activation
    std::array<AdamState, 256> w_opt;   // Weight optimizers
    std::array<AdamState, 16> b_opt;    // Bias optimizers

    Layer16() {
        // Initialize with golden ratio scaling
        for (size_t i = 0; i < 256; ++i) {
            weights[i] = (static_cast<double>(i % 7) - 3.0) * PHI_INV * 0.1;
        }
        biases.fill(0.0);
        output.fill(0.0);
        delta.fill(0.0);
        pre_act.fill(0.0);
    }

    /**
     * @brief Forward pass with tanh activation
     */
    void forward(const double* input) {
        for (int i = 0; i < 16; ++i) {
            double sum = biases[i];
            for (int j = 0; j < 16; ++j) {
                sum += weights[i * 16 + j] * input[j];
            }
            pre_act[i] = sum;
            output[i] = math::Tanh::forward(sum);
        }
    }

    /**
     * @brief Backward pass: compute deltas and update weights
     */
    void backward(const double* input, const double* next_delta = nullptr,
                  const double* next_weights = nullptr, int next_size = 0) {
        // Compute deltas
        for (int i = 0; i < 16; ++i) {
            double derivative = math::Tanh::backward(output[i]);
            if (next_delta && next_weights) {
                // Hidden layer: backprop from next layer
                double sum = 0.0;
                for (int k = 0; k < next_size; ++k) {
                    sum += next_delta[k] * next_weights[k * 16 + i];
                }
                delta[i] = sum * derivative;
            }
            // If output layer, delta should be set externally
        }

        // Update weights
        for (int i = 0; i < 16; ++i) {
            for (int j = 0; j < 16; ++j) {
                double grad = weight_gradient(delta[i], input[j]);
                weights[i * 16 + j] += w_opt[i * 16 + j].update(grad);
            }
            biases[i] += b_opt[i].update(delta[i]);
        }
    }

    /**
     * @brief Hebbian update (unsupervised)
     */
    void hebbian_update(const double* input) {
        for (int i = 0; i < 16; ++i) {
            for (int j = 0; j < 16; ++j) {
                double dw = hebbian_normalized(input[j], output[i], weights[i * 16 + j]);
                weights[i * 16 + j] += dw;
            }
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// V. REWARD-MODULATED LEARNING
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Three-factor learning: Δw = η·x_pre·x_post·R
 * R is the reward/neuromodulator signal
 */
inline double three_factor_learning(double x_pre, double x_post, double reward,
                                      double eligibility = 1.0, double eta = ETA) {
    return eta * x_pre * x_post * reward * eligibility;
}

/**
 * @brief Eligibility trace update: e(t+1) = λ·e(t) + x_pre·x_post
 */
inline double eligibility_update(double trace, double x_pre, double x_post, double lambda = 0.9) {
    return lambda * trace + x_pre * x_post;
}

/**
 * @brief TD error: δ = r + γV(s') - V(s)
 */
inline double td_error(double reward, double value_next, double value_current, double gamma_discount = 0.99) {
    return reward + gamma_discount * value_next - value_current;
}

/**
 * @brief Policy gradient: ∇J ≈ Σ(∇log π(a|s)·R)
 */
inline double policy_gradient(double action_prob, double advantage) {
    if (action_prob < 1e-10) return 0.0;
    return advantage / action_prob;
}

// ═══════════════════════════════════════════════════════════════════════════════
// VI. ASSOCIATIVE MEMORY (HOPFIELD-LIKE)
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Store pattern in weight matrix via outer product
 */
inline void store_pattern(double* weights, const double* pattern, size_t n) {
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            if (i != j) {
                weights[i * n + j] += pattern[i] * pattern[j];
            }
        }
    }
}

/**
 * @brief Recall pattern from associative memory (one step)
 */
inline void recall_step(const double* weights, const double* input, double* output, size_t n) {
    for (size_t i = 0; i < n; ++i) {
        double sum = 0.0;
        for (size_t j = 0; j < n; ++j) {
            sum += weights[i * n + j] * input[j];
        }
        output[i] = (sum > 0) ? 1.0 : -1.0;
    }
}

// ═══════════════════════════════════════════════════════════════════════════════
// VII. MIND LEARNING STATE
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Compact learning state for each mind (~200 bytes)
 */
struct MindLearner {
    // 16-element weight vector for consciousness state
    std::array<double, 16> weights;
    std::array<double, 16> eligibility;
    std::array<AdamState, 16> optimizers;

    double bcm_theta = 0.5;        // BCM threshold
    double reward_baseline = 0.0;  // Running average of reward
    double learning_rate = ETA;

    MindLearner() {
        weights.fill(PHI_INV * 0.1);
        eligibility.fill(0.0);
    }

    /**
     * @brief Process input and produce output
     */
    double forward(const double* input, size_t n = 16) {
        double sum = 0.0;
        size_t count = std::min(n, size_t(16));
        for (size_t i = 0; i < count; ++i) {
            sum += weights[i] * input[i];
        }
        return math::Tanh::forward(sum);
    }

    /**
     * @brief Update with reward signal (RL-style)
     */
    void update_with_reward(const double* input, double output, double reward, size_t n = 16) {
        // TD-like reward baseline update
        double delta = reward - reward_baseline;
        reward_baseline += 0.01 * delta;

        size_t count = std::min(n, size_t(16));
        for (size_t i = 0; i < count; ++i) {
            // Update eligibility trace
            eligibility[i] = eligibility_update(eligibility[i], input[i], output);

            // Three-factor update
            double dw = three_factor_learning(input[i], output, delta, eligibility[i], learning_rate);
            weights[i] += dw;

            // Clamp weights
            weights[i] = clamp(weights[i], -2.0, 2.0);
        }
    }

    /**
     * @brief Hebbian update (unsupervised)
     */
    void update_hebbian(const double* input, double output, size_t n = 16) {
        size_t count = std::min(n, size_t(16));
        for (size_t i = 0; i < count; ++i) {
            double dw = hebbian_bcm(input[i], output, bcm_theta, learning_rate);
            weights[i] += dw;
            weights[i] = clamp(weights[i], -2.0, 2.0);
        }
        // Update BCM threshold
        bcm_theta = bcm_threshold_update(bcm_theta, output);
    }

    /**
     * @brief Gradient update (supervised)
     */
    void update_gradient(const double* input, double output, double target, size_t n = 16) {
        double derivative = math::Tanh::backward(output);
        double delta = output_delta(output, target, derivative);

        size_t count = std::min(n, size_t(16));
        for (size_t i = 0; i < count; ++i) {
            double grad = weight_gradient(delta, input[i]);
            weights[i] += optimizers[i].update(grad);
        }
    }
};

} // namespace learning
} // namespace mind_farm