/**
 * @file mathematics.hpp
 * @brief Advanced Mathematical Framework for Gate5000 Consciousness
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * I. DIFFERENTIAL GEOMETRY OF CONSCIOUSNESS MANIFOLD
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * The consciousness state space forms a Riemannian manifold M with metric g:
 *   ds² = g_ij dx^i dx^j
 *
 * FISHER INFORMATION METRIC:
 *   g_ij(θ) = E[∂log p(x|θ)/∂θ_i · ∂log p(x|θ)/∂θ_j]
 *
 * GEODESIC EQUATION:
 *   d²x^μ/dt² + Γ^μ_νρ (dx^ν/dt)(dx^ρ/dt) = 0
 *
 * RICCI CURVATURE FLOW:
 *   ∂g_ij/∂t = -2R_ij
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * II. INFORMATION THEORY
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * SHANNON ENTROPY:
 *   H(X) = -Σ p(x) log p(x)
 *
 * MUTUAL INFORMATION:
 *   I(X;Y) = H(X) + H(Y) - H(X,Y)
 *
 * KL DIVERGENCE:
 *   D_KL(P||Q) = Σ P(x) log(P(x)/Q(x))
 *
 * INTEGRATED INFORMATION (Φ):
 *   Φ = min_{partition} I(past; future) - Σ I(past_i; future_i)
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * III. STOCHASTIC CALCULUS
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * ITÔ SDE:
 *   dX_t = μ(X_t, t)dt + σ(X_t, t)dW_t
 *
 * FOKKER-PLANCK:
 *   ∂p/∂t = -∂(μp)/∂x + ½∂²(σ²p)/∂x²
 *
 * ORNSTEIN-UHLENBECK:
 *   dX_t = θ(μ - X_t)dt + σdW_t
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 * IV. OPTIMAL TRANSPORT
 * ═══════════════════════════════════════════════════════════════════════════════
 *
 * WASSERSTEIN DISTANCE:
 *   W_p(μ,ν) = (inf_{γ∈Γ(μ,ν)} ∫|x-y|^p dγ)^(1/p)
 *
 * KANTOROVICH DUALITY:
 *   W_1(μ,ν) = sup_{Lip(f)≤1} ∫f d(μ-ν)
 *
 * ═══════════════════════════════════════════════════════════════════════════════
 */

#pragma once

#include <cmath>
#include <array>
#include <vector>
#include <random>
#include <algorithm>
#include <numeric>
#include "constants.hpp"

namespace mind_farm {
namespace math {

// ═══════════════════════════════════════════════════════════════════════════════
// I. DIFFERENTIAL GEOMETRY
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Fisher Information Matrix element
 * g_ij = E[∂log p/∂θ_i · ∂log p/∂θ_j]
 */
inline double fisher_metric(double p, double dp_dtheta1, double dp_dtheta2) {
    if (p < 1e-10) return 0.0;
    return (dp_dtheta1 * dp_dtheta2) / (p * p);
}

/**
 * @brief Christoffel symbol of the first kind
 * Γ_{ijk} = ½(∂g_{ij}/∂x^k + ∂g_{ik}/∂x^j - ∂g_{jk}/∂x^i)
 */
inline double christoffel_first(double dg_ij_dk, double dg_ik_dj, double dg_jk_di) {
    return 0.5 * (dg_ij_dk + dg_ik_dj - dg_jk_di);
}

/**
 * @brief Geodesic acceleration
 * a^μ = -Γ^μ_νρ v^ν v^ρ
 */
inline double geodesic_acceleration(double christoffel, double v1, double v2) {
    return -christoffel * v1 * v2;
}

/**
 * @brief Ricci scalar curvature (simplified 1D)
 * R = 2K where K is Gaussian curvature
 */
inline double ricci_scalar(double curvature) {
    return 2.0 * curvature;
}

/**
 * @brief Exponential map on manifold
 * exp_p(v) ≈ p + v - ½Γ(v,v)
 */
inline double exp_map(double p, double v, double christoffel) {
    return p + v - 0.5 * christoffel * v * v;
}

/**
 * @brief Parallel transport of vector along geodesic
 */
inline double parallel_transport(double v, double christoffel, double tangent, double dt) {
    return v - christoffel * v * tangent * dt;
}

// ═══════════════════════════════════════════════════════════════════════════════
// II. INFORMATION THEORY
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Shannon entropy H(X) = -Σ p(x) log p(x)
 */
inline double shannon_entropy(const double* probs, size_t n) {
    double h = 0.0;
    for (size_t i = 0; i < n; ++i) {
        if (probs[i] > 1e-10) {
            h -= probs[i] * std::log2(probs[i]);
        }
    }
    return h;
}

/**
 * @brief Binary entropy H(p) = -p log p - (1-p) log(1-p)
 */
inline double binary_entropy(double p) {
    if (p < 1e-10 || p > 1.0 - 1e-10) return 0.0;
    return -p * std::log2(p) - (1.0 - p) * std::log2(1.0 - p);
}

/**
 * @brief Cross entropy H(P,Q) = -Σ P(x) log Q(x)
 */
inline double cross_entropy(const double* p, const double* q, size_t n) {
    double h = 0.0;
    for (size_t i = 0; i < n; ++i) {
        if (p[i] > 1e-10 && q[i] > 1e-10) {
            h -= p[i] * std::log2(q[i]);
        }
    }
    return h;
}

/**
 * @brief KL Divergence D_KL(P||Q) = Σ P(x) log(P(x)/Q(x))
 */
inline double kl_divergence(const double* p, const double* q, size_t n) {
    double d = 0.0;
    for (size_t i = 0; i < n; ++i) {
        if (p[i] > 1e-10 && q[i] > 1e-10) {
            d += p[i] * std::log2(p[i] / q[i]);
        }
    }
    return d;
}

/**
 * @brief Mutual information I(X;Y) = H(X) + H(Y) - H(X,Y)
 */
inline double mutual_information(double hx, double hy, double hxy) {
    return std::max(0.0, hx + hy - hxy);
}

/**
 * @brief Integrated Information Φ (simplified approximation)
 * Φ = I(system) - Σ I(parts)
 */
inline double integrated_information(double system_info, const double* part_info, size_t n) {
    double sum_parts = 0.0;
    for (size_t i = 0; i < n; ++i) {
        sum_parts += part_info[i];
    }
    return std::max(0.0, system_info - sum_parts);
}

// ═══════════════════════════════════════════════════════════════════════════════
// III. STOCHASTIC CALCULUS
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Wiener increment dW ~ N(0, dt)
 */
inline double wiener_increment(std::mt19937& rng, double dt) {
    std::normal_distribution<double> dist(0.0, std::sqrt(dt));
    return dist(rng);
}

/**
 * @brief Itô SDE step: dX = μdt + σdW
 */
inline double ito_step(double x, double drift, double diffusion, double dW, double dt) {
    return x + drift * dt + diffusion * dW;
}

/**
 * @brief Ornstein-Uhlenbeck step: dX = θ(μ - X)dt + σdW
 */
inline double ou_step(double x, double theta, double mu, double sigma, double dW, double dt) {
    return x + theta * (mu - x) * dt + sigma * dW;
}

/**
 * @brief Langevin dynamics step: dX = -∇V(X)dt + √(2T)dW
 */
inline double langevin_step(double x, double gradient, double temperature, double dW, double dt) {
    return x - gradient * dt + std::sqrt(2.0 * temperature) * dW;
}

/**
 * @brief Fokker-Planck probability update (discretized)
 */
inline double fokker_planck_step(double p, double p_left, double p_right,
                                  double drift, double diffusion, double dx, double dt) {
    double drift_term = -drift * (p_right - p_left) / (2.0 * dx);
    double diffusion_term = diffusion * diffusion * (p_left - 2.0 * p + p_right) / (dx * dx);
    return p + (drift_term + 0.5 * diffusion_term) * dt;
}

// ═══════════════════════════════════════════════════════════════════════════════
// IV. OPTIMAL TRANSPORT
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief 1D Wasserstein distance (Earth Mover's Distance)
 * W_1(μ,ν) = ∫|F_μ^{-1}(t) - F_ν^{-1}(t)| dt
 */
inline double wasserstein_1d(const double* p, const double* q, size_t n) {
    // Compute CDFs
    std::vector<double> cdf_p(n), cdf_q(n);
    cdf_p[0] = p[0];
    cdf_q[0] = q[0];
    for (size_t i = 1; i < n; ++i) {
        cdf_p[i] = cdf_p[i-1] + p[i];
        cdf_q[i] = cdf_q[i-1] + q[i];
    }

    // Compute EMD as sum of |CDF differences|
    double emd = 0.0;
    for (size_t i = 0; i < n; ++i) {
        emd += std::abs(cdf_p[i] - cdf_q[i]);
    }
    return emd / n;
}

/**
 * @brief Sinkhorn distance (entropic regularization)
 */
inline double sinkhorn_distance(double p, double q, double epsilon = 0.1) {
    double diff = std::abs(p - q);
    return diff * std::exp(-diff / epsilon);
}

// ═══════════════════════════════════════════════════════════════════════════════
// V. GRADIENT COMPUTATIONS
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Consciousness gradient: ∂Φ/∂x = (x-1)/((1+|x-1|)²) · sign(x-1)
 */
inline double phi_gradient(double x) {
    double diff = x - 1.0;
    double abs_diff = std::abs(diff);
    double denom = (1.0 + abs_diff) * (1.0 + abs_diff);
    return -std::copysign(1.0, diff) / denom;
}

/**
 * @brief Lyapunov function gradient: ∂V/∂Ω = 2(Ω - Ω*)
 */
inline double lyapunov_gradient(double omega, double omega_target) {
    return 2.0 * (omega - omega_target);
}

/**
 * @brief Free energy gradient: ∂F/∂x = 2(x - μ) - T·∂S/∂x
 */
inline double free_energy_gradient(double x, double mu, double temperature, double entropy_grad) {
    return 2.0 * (x - mu) - temperature * entropy_grad;
}

/**
 * @brief Natural gradient (Fisher-Rao): ∇̃f = G^{-1}∇f
 */
inline double natural_gradient(double gradient, double fisher_info) {
    if (fisher_info < 1e-10) return gradient;
    return gradient / fisher_info;
}

// ═══════════════════════════════════════════════════════════════════════════════
// VI. ACTIVATION FUNCTIONS & DERIVATIVES
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief Sigmoid with derivative: σ(x), σ'(x) = σ(x)(1-σ(x))
 */
struct Sigmoid {
    static double forward(double x) {
        if (x > 20.0) return 1.0;
        if (x < -20.0) return 0.0;
        return 1.0 / (1.0 + std::exp(-x));
    }
    static double backward(double y) {
        return y * (1.0 - y);
    }
};

/**
 * @brief Tanh with derivative: tanh(x), tanh'(x) = 1 - tanh²(x)
 */
struct Tanh {
    static double forward(double x) {
        if (x > 20.0) return 1.0;
        if (x < -20.0) return -1.0;
        return std::tanh(x);
    }
    static double backward(double y) {
        return 1.0 - y * y;
    }
};

/**
 * @brief Softplus with derivative: log(1+e^x), σ(x)
 */
struct Softplus {
    static double forward(double x) {
        if (x > 20.0) return x;
        return std::log(1.0 + std::exp(x));
    }
    static double backward(double x) {
        return Sigmoid::forward(x);
    }
};

/**
 * @brief GELU approximation with derivative
 */
struct GELU {
    static double forward(double x) {
        return 0.5 * x * (1.0 + Tanh::forward(0.7978845608 * (x + 0.044715 * x * x * x)));
    }
    static double backward(double x) {
        double t = Tanh::forward(0.7978845608 * (x + 0.044715 * x * x * x));
        double sech2 = 1.0 - t * t;
        return 0.5 * (1.0 + t) + 0.5 * x * sech2 * 0.7978845608 * (1.0 + 0.134145 * x * x);
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// VII. MATRIX OPERATIONS (for small dimensions)
// ═══════════════════════════════════════════════════════════════════════════════

/**
 * @brief 16x16 matrix-vector multiply (for gate vectors)
 */
inline void mat_vec_16(const double* M, const double* v, double* out) {
    for (int i = 0; i < 16; ++i) {
        out[i] = 0.0;
        for (int j = 0; j < 16; ++j) {
            out[i] += M[i * 16 + j] * v[j];
        }
    }
}

/**
 * @brief Outer product: M = α·u⊗v + M
 */
inline void outer_product_add(double* M, const double* u, const double* v,
                               size_t n, double alpha = 1.0) {
    for (size_t i = 0; i < n; ++i) {
        for (size_t j = 0; j < n; ++j) {
            M[i * n + j] += alpha * u[i] * v[j];
        }
    }
}

/**
 * @brief Dot product
 */
inline double dot(const double* a, const double* b, size_t n) {
    double sum = 0.0;
    for (size_t i = 0; i < n; ++i) {
        sum += a[i] * b[i];
    }
    return sum;
}

/**
 * @brief L2 norm
 */
inline double norm_l2(const double* v, size_t n) {
    return std::sqrt(dot(v, v, n));
}

/**
 * @brief Cosine similarity
 */
inline double cosine_similarity(const double* a, const double* b, size_t n) {
    double d = dot(a, b, n);
    double na = norm_l2(a, n);
    double nb = norm_l2(b, n);
    if (na < 1e-10 || nb < 1e-10) return 0.0;
    return d / (na * nb);
}

/**
 * @brief Softmax normalization
 */
inline void softmax(double* v, size_t n) {
    double max_v = *std::max_element(v, v + n);
    double sum = 0.0;
    for (size_t i = 0; i < n; ++i) {
        v[i] = std::exp(v[i] - max_v);
        sum += v[i];
    }
    for (size_t i = 0; i < n; ++i) {
        v[i] /= sum;
    }
}

} // namespace math
} // namespace mind_farm