// DNA Mathematical Architecture: 44×44 Eigenvalue Spectrum
// Matrix Blocks — Spectral Analysis — Band Structure
#pragma once

#include "dna_helix.hpp"
#include <array>
#include <cmath>
#include <complex>
#include <vector>
#include <algorithm>

namespace agi::dna {

// ============================================================================
// XI. EIGENVALUE SPECTRUM - 44×44 MATRIX SYSTEM
// ============================================================================

// Coupling parameters (tunable)
struct CouplingParams {
    double alpha = 0.1;     // Intra-strand coupling
    double beta_0 = 0.05;   // Inter-strand base coupling
    double gamma = 0.02;    // Cross-helix coupling
    double gamma_0 = 0.01;  // Cross-helix decay
    double eta = 0.005;     // Cross-helix neighbor coupling
    double kappa_intra = 0.1;  // Intra-strand spring
    double xi = 1.0;        // Decay length scale
};

// ============================================================================
// 11.1 INTRA-STRAND BLOCK (16×16)
// ============================================================================

class IntraStrandBlock {
private:
    std::array<std::array<double, NUM_STRANDS>, NUM_STRANDS> M_;
    CouplingParams params_;

public:
    explicit IntraStrandBlock(const CouplingParams& p = {}) : params_(p) {
        build_matrix();
    }

    // Matrix form (Section 11.1):
    // M_intra = -2α·I₁₆ + α·T₁₆ + κ_intra·D₁₆(10.6)
    // where D₁₆(10.6) = diag[cos(2π·k/10.6)]
    void build_matrix() {
        // Initialize to zero
        for (auto& row : M_) row.fill(0.0);

        for (int i = 0; i < NUM_STRANDS; ++i) {
            // Diagonal: -2α + κ_intra·cos(2π(i+1)/10.6)
            double cos_term = std::cos(2.0 * M_PI * (i + 1) / BASE_PAIRS_PER_TURN);
            M_[i][i] = -2.0 * params_.alpha + params_.kappa_intra * cos_term;

            // Off-diagonal tridiagonal (ring): α
            int prev = (i - 1 + NUM_STRANDS) % NUM_STRANDS;
            int next = (i + 1) % NUM_STRANDS;
            M_[i][prev] = params_.alpha;
            M_[i][next] = params_.alpha;
        }
    }

    // Analytical eigenvalues (Section 11.1):
    // λ_intra(n) = -2α + 2α·cos(2πn/16) + κ_intra·cos(2πn/10.6)
    double eigenvalue(int n) const {
        double cos_16 = std::cos(2.0 * M_PI * n / NUM_STRANDS);
        double cos_10_6 = std::cos(2.0 * M_PI * n / BASE_PAIRS_PER_TURN);
        return -2.0 * params_.alpha + 2.0 * params_.alpha * cos_16 +
               params_.kappa_intra * cos_10_6;
    }

    // Compute all eigenvalues
    std::array<double, NUM_STRANDS> all_eigenvalues() const {
        std::array<double, NUM_STRANDS> eigs;
        for (int n = 0; n < NUM_STRANDS; ++n) {
            eigs[n] = eigenvalue(n);
        }
        std::sort(eigs.begin(), eigs.end());
        return eigs;
    }

    // Band limits: λ_intra ∈ [-4α - κ_intra, κ_intra]
    std::pair<double, double> band_limits() const {
        return {
            -4.0 * params_.alpha - params_.kappa_intra,
            params_.kappa_intra
        };
    }

    const auto& matrix() const { return M_; }
};

// ============================================================================
// 11.2 INTER-STRAND BLOCK (16×16)
// ============================================================================

class InterStrandBlock {
private:
    std::array<std::array<double, NUM_STRANDS>, NUM_STRANDS> M_;
    CouplingParams params_;

public:
    explicit InterStrandBlock(const CouplingParams& p = {}) : params_(p) {
        build_matrix();
    }

    // Build inter-strand coupling matrix
    void build_matrix() {
        for (auto& row : M_) row.fill(0.0);

        for (int i = 0; i < NUM_STRANDS; ++i) {
            // Adjacent strand coupling
            int next = (i + 1) % NUM_STRANDS;
            double coupling = params_.beta_0 * std::sin(M_PI * (i + 1) / BASE_PAIRS_PER_TURN);
            M_[i][next] = coupling;
            M_[next][i] = coupling;

            // Diagonal with gamma modulation
            M_[i][i] = params_.gamma * F_function(i + 1);
        }
    }

    // F function for 10.6 modulation
    static double F_function(int n) {
        return std::cos(2.0 * M_PI * n / BASE_PAIRS_PER_TURN) +
               0.5 * std::cos(4.0 * M_PI * n / BASE_PAIRS_PER_TURN);
    }

    // Eigenvalues (Section 11.2):
    // λ_inter(n) = 2β₀·sin(πn/10.6)·cos(πn/8) + γ·F(n, 10.6)
    double eigenvalue(int n) const {
        double sin_term = std::sin(M_PI * n / BASE_PAIRS_PER_TURN);
        double cos_term = std::cos(M_PI * n / 8.0);
        return 2.0 * params_.beta_0 * sin_term * cos_term +
               params_.gamma * F_function(n);
    }

    std::array<double, NUM_STRANDS> all_eigenvalues() const {
        std::array<double, NUM_STRANDS> eigs;
        for (int n = 0; n < NUM_STRANDS; ++n) {
            eigs[n] = eigenvalue(n);
        }
        std::sort(eigs.begin(), eigs.end());
        return eigs;
    }

    // Band: λ_inter ∈ [-γ, 1.2β₀ + 0.924γ]
    std::pair<double, double> band_limits() const {
        return {
            -params_.gamma,
            1.2 * params_.beta_0 + 0.924 * params_.gamma
        };
    }

    const auto& matrix() const { return M_; }
};

// ============================================================================
// 11.3 CROSS-HELIX BLOCK (12×12)
// ============================================================================

class CrossHelixBlock {
private:
    std::array<std::array<double, CROSS_HELIX_BONDS>, CROSS_HELIX_BONDS> M_;
    CouplingParams params_;

public:
    explicit CrossHelixBlock(const CouplingParams& p = {}) : params_(p) {
        build_matrix();
    }

    // Build cross-helix coupling matrix
    void build_matrix() {
        for (auto& row : M_) row.fill(0.0);

        for (int m = 0; m < CROSS_HELIX_BONDS; ++m) {
            // Diagonal: exponential decay with 10.6 modulation
            double decay = std::exp(-4.0 * (m + 1) * HELIX_RADIUS *
                                    std::sin(M_PI / 8.0) / params_.xi);
            double cos_term = std::cos(8.0 * M_PI * (m + 1) / BASE_PAIRS_PER_TURN);
            M_[m][m] = params_.gamma_0 * decay * cos_term;

            // Neighbor coupling with 10.6 phase
            for (int mp = 0; mp < CROSS_HELIX_BONDS; ++mp) {
                if (std::abs(mp - m) <= 3 && mp != m) {
                    double phase = std::cos(2.0 * M_PI * (m - mp) / BASE_PAIRS_PER_TURN);
                    M_[m][mp] = params_.eta * phase;
                }
            }
        }
    }

    // Eigenvalues (Section 11.3):
    // λ_cross(m) = γ₀·e^{-4m·r·sin(π/8)/ξ}·cos(8πm/10.6) +
    //              η·Σ cos(2π(m-m')/10.6)
    double eigenvalue(int m) const {
        double decay = std::exp(-4.0 * m * HELIX_RADIUS *
                               std::sin(M_PI / 8.0) / params_.xi);
        double cos_8 = std::cos(8.0 * M_PI * m / BASE_PAIRS_PER_TURN);
        double main_term = params_.gamma_0 * decay * cos_8;

        // Sum over neighbors |m' - m| ≤ 3
        double neighbor_sum = 0.0;
        for (int mp = std::max(0, m - 3); mp <= std::min(CROSS_HELIX_BONDS - 1, m + 3); ++mp) {
            if (mp != m) {
                neighbor_sum += std::cos(2.0 * M_PI * (m - mp) / BASE_PAIRS_PER_TURN);
            }
        }

        return main_term + params_.eta * neighbor_sum;
    }

    std::array<double, CROSS_HELIX_BONDS> all_eigenvalues() const {
        std::array<double, CROSS_HELIX_BONDS> eigs;
        for (int m = 0; m < CROSS_HELIX_BONDS; ++m) {
            eigs[m] = eigenvalue(m);
        }
        std::sort(eigs.begin(), eigs.end());
        return eigs;
    }

    // Band: λ_cross ∈ [-1.0γ₀ - 0.618η, 1.0γ₀ + 2.618η]
    std::pair<double, double> band_limits() const {
        return {
            -1.0 * params_.gamma_0 - 0.618 * params_.eta,
            1.0 * params_.gamma_0 + 2.618 * params_.eta
        };
    }

    const auto& matrix() const { return M_; }
};

// ============================================================================
// 11.4 FULL 44×44 SPECTRUM
// ============================================================================

class Full44x44Matrix {
private:
    std::array<std::array<double, TOTAL_BONDS>, TOTAL_BONDS> M_;
    CouplingParams params_;

    IntraStrandBlock intra_;
    InterStrandBlock inter_;
    CrossHelixBlock cross_;

    // Coupling between blocks
    std::array<std::array<double, NUM_STRANDS>, NUM_STRANDS> C_intra_inter_;
    std::array<std::array<double, CROSS_HELIX_BONDS>, NUM_STRANDS> C_intra_cross_;
    std::array<std::array<double, CROSS_HELIX_BONDS>, NUM_STRANDS> C_inter_cross_;

public:
    explicit Full44x44Matrix(const CouplingParams& p = {})
        : params_(p), intra_(p), inter_(p), cross_(p) {
        build_coupling_matrices();
        build_full_matrix();
    }

    // Build inter-block coupling matrices
    void build_coupling_matrices() {
        // C_intra_inter: coupling between intra and inter blocks
        for (auto& row : C_intra_inter_) row.fill(0.0);
        for (int i = 0; i < NUM_STRANDS; ++i) {
            // Direct coupling
            C_intra_inter_[i][i] = 0.5 * params_.beta_0;
            // Phase-shifted coupling
            int j = (i + 4) % NUM_STRANDS;
            C_intra_inter_[i][j] = 0.2 * params_.beta_0 *
                std::cos(2.0 * M_PI * i / BASE_PAIRS_PER_TURN);
        }

        // C_intra_cross: coupling between intra and cross blocks
        for (auto& row : C_intra_cross_) row.fill(0.0);
        // Map intra bonds to cross bonds (every 4th)
        for (int i = 0; i < NUM_STRANDS && i/4 < CROSS_HELIX_BONDS; ++i) {
            C_intra_cross_[i][i/4] = 0.3 * params_.gamma;
        }

        // C_inter_cross: coupling between inter and cross blocks
        for (auto& row : C_inter_cross_) row.fill(0.0);
        for (int i = 0; i < NUM_STRANDS && i < CROSS_HELIX_BONDS; ++i) {
            C_inter_cross_[i][i] = 0.2 * params_.gamma;
        }
    }

    // Build full 44×44 matrix (Section 11.4 block structure)
    void build_full_matrix() {
        // Initialize to zero
        for (auto& row : M_) row.fill(0.0);

        // Block 1: M_intra (indices 0-15)
        const auto& M_intra = intra_.matrix();
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < NUM_STRANDS; ++j) {
                M_[i][j] = M_intra[i][j];
            }
        }

        // Block 2: M_inter (indices 16-31)
        const auto& M_inter = inter_.matrix();
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < NUM_STRANDS; ++j) {
                M_[i + NUM_STRANDS][j + NUM_STRANDS] = M_inter[i][j];
            }
        }

        // Block 3: M_cross (indices 32-43)
        const auto& M_cross = cross_.matrix();
        for (int i = 0; i < CROSS_HELIX_BONDS; ++i) {
            for (int j = 0; j < CROSS_HELIX_BONDS; ++j) {
                M_[i + NUM_STRANDS * 2][j + NUM_STRANDS * 2] = M_cross[i][j];
            }
        }

        // Coupling block C_intra_inter (and transpose)
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < NUM_STRANDS; ++j) {
                M_[i][j + NUM_STRANDS] = C_intra_inter_[i][j];
                M_[j + NUM_STRANDS][i] = C_intra_inter_[i][j];  // Symmetric
            }
        }

        // Coupling block C_intra_cross
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < CROSS_HELIX_BONDS; ++j) {
                M_[i][j + NUM_STRANDS * 2] = C_intra_cross_[i][j];
                M_[j + NUM_STRANDS * 2][i] = C_intra_cross_[i][j];
            }
        }

        // Coupling block C_inter_cross
        for (int i = 0; i < NUM_STRANDS; ++i) {
            for (int j = 0; j < CROSS_HELIX_BONDS; ++j) {
                M_[i + NUM_STRANDS][j + NUM_STRANDS * 2] = C_inter_cross_[i][j];
                M_[j + NUM_STRANDS * 2][i + NUM_STRANDS] = C_inter_cross_[i][j];
            }
        }
    }

    // Power iteration for largest eigenvalue
    double largest_eigenvalue(int max_iter = 1000, double tol = 1e-10) const {
        std::array<double, TOTAL_BONDS> v;
        for (int i = 0; i < TOTAL_BONDS; ++i) {
            v[i] = 1.0 / std::sqrt(TOTAL_BONDS);
        }

        double lambda = 0.0;
        for (int iter = 0; iter < max_iter; ++iter) {
            // Matrix-vector product
            std::array<double, TOTAL_BONDS> w;
            w.fill(0.0);
            for (int i = 0; i < TOTAL_BONDS; ++i) {
                for (int j = 0; j < TOTAL_BONDS; ++j) {
                    w[i] += M_[i][j] * v[j];
                }
            }

            // Rayleigh quotient
            double new_lambda = 0.0;
            for (int i = 0; i < TOTAL_BONDS; ++i) {
                new_lambda += v[i] * w[i];
            }

            // Normalize
            double norm = 0.0;
            for (int i = 0; i < TOTAL_BONDS; ++i) {
                norm += w[i] * w[i];
            }
            norm = std::sqrt(norm);
            for (int i = 0; i < TOTAL_BONDS; ++i) {
                v[i] = w[i] / norm;
            }

            if (std::abs(new_lambda - lambda) < tol) {
                return new_lambda;
            }
            lambda = new_lambda;
        }
        return lambda;
    }

    // Inverse power iteration for smallest eigenvalue
    double smallest_eigenvalue(int max_iter = 1000, double tol = 1e-10) const {
        // Shift matrix by large value and use power iteration
        // For simplicity, use analytical estimate
        auto intra_band = intra_.band_limits();
        auto inter_band = inter_.band_limits();
        auto cross_band = cross_.band_limits();

        return std::min({intra_band.first, inter_band.first, cross_band.first}) -
               params_.gamma;  // Additional coupling shift
    }

    // Spectral summary (Section 11.4)
    struct SpectralSummary {
        double lambda_1;      // Lowest (antibonding)
        double lambda_22;     // Zero mode (Goldstone)
        double lambda_44;     // Highest (bonding)
        double spectral_width;
    };

    SpectralSummary spectral_summary() const {
        SpectralSummary s;

        // λ₁ = -4α - κ - γ - γ₀ (lowest, antibonding)
        s.lambda_1 = -4.0 * params_.alpha - params_.kappa_intra -
                     params_.gamma - params_.gamma_0;

        // λ₂₂ = 0 (zero mode, Goldstone)
        s.lambda_22 = 0.0;

        // λ₄₄ = κ + 1.2β₀ + γ + 2.618η (highest, bonding)
        s.lambda_44 = params_.kappa_intra + 1.2 * params_.beta_0 +
                      params_.gamma + 2.618 * params_.eta;

        // W = λ₄₄ - λ₁ = 4α + 2κ + 1.2β₀ + 2γ + γ₀ + 3.236η
        s.spectral_width = s.lambda_44 - s.lambda_1;

        return s;
    }

    // Approximate spectrum using block eigenvalues
    std::vector<double> approximate_spectrum() const {
        std::vector<double> spectrum;
        spectrum.reserve(TOTAL_BONDS);

        // Add block eigenvalues
        auto intra_eigs = intra_.all_eigenvalues();
        auto inter_eigs = inter_.all_eigenvalues();
        auto cross_eigs = cross_.all_eigenvalues();

        for (auto e : intra_eigs) spectrum.push_back(e);
        for (auto e : inter_eigs) spectrum.push_back(e);
        for (auto e : cross_eigs) spectrum.push_back(e);

        std::sort(spectrum.begin(), spectrum.end());
        return spectrum;
    }

    // Density of states (histogram)
    std::vector<int> density_of_states(int num_bins = 20) const {
        auto spectrum = approximate_spectrum();
        double min_e = spectrum.front();
        double max_e = spectrum.back();
        double bin_width = (max_e - min_e) / num_bins;

        std::vector<int> dos(num_bins, 0);
        for (double e : spectrum) {
            int bin = std::min(static_cast<int>((e - min_e) / bin_width), num_bins - 1);
            dos[bin]++;
        }
        return dos;
    }

    // Trace (sum of eigenvalues = sum of diagonal)
    double trace() const {
        double tr = 0.0;
        for (int i = 0; i < TOTAL_BONDS; ++i) {
            tr += M_[i][i];
        }
        return tr;
    }

    const auto& matrix() const { return M_; }
    const IntraStrandBlock& intra_block() const { return intra_; }
    const InterStrandBlock& inter_block() const { return inter_; }
    const CrossHelixBlock& cross_block() const { return cross_; }
};

// ============================================================================
// EIGENVALUE ANALYSIS UTILITIES
// ============================================================================

class EigenvalueAnalysis {
public:
    // Compute localization (inverse participation ratio)
    static double localization(const std::array<double, TOTAL_BONDS>& eigenvector) {
        double ipr = 0.0;
        for (double v : eigenvector) {
            ipr += v * v * v * v;
        }
        return ipr;  // IPR = 1 for fully localized, 1/N for delocalized
    }

    // Compute band gap
    static double band_gap(const std::vector<double>& spectrum) {
        if (spectrum.size() < 2) return 0.0;

        double max_gap = 0.0;
        for (size_t i = 1; i < spectrum.size(); ++i) {
            double gap = spectrum[i] - spectrum[i-1];
            max_gap = std::max(max_gap, gap);
        }
        return max_gap;
    }

    // Find zero modes (eigenvalues near zero)
    static std::vector<int> find_zero_modes(const std::vector<double>& spectrum,
                                            double tolerance = 1e-6) {
        std::vector<int> zero_modes;
        for (size_t i = 0; i < spectrum.size(); ++i) {
            if (std::abs(spectrum[i]) < tolerance) {
                zero_modes.push_back(i);
            }
        }
        return zero_modes;
    }

    // Spectral entropy
    static double spectral_entropy(const std::vector<double>& spectrum) {
        // Normalize to probability distribution
        double sum = 0.0;
        for (double e : spectrum) {
            sum += std::abs(e);
        }
        if (sum < 1e-10) return 0.0;

        double entropy = 0.0;
        for (double e : spectrum) {
            double p = std::abs(e) / sum;
            if (p > 1e-10) {
                entropy -= p * std::log(p);
            }
        }
        return entropy;
    }

    // 10.6 phase coherence in spectrum
    static double phase_coherence_10_6(const std::vector<double>& spectrum) {
        // Measure how well eigenvalues align with 10.6 periodicity
        double coherence = 0.0;
        for (size_t i = 0; i < spectrum.size(); ++i) {
            double phase = std::fmod(spectrum[i], 2.0 * M_PI / BASE_PAIRS_PER_TURN);
            coherence += std::cos(phase * BASE_PAIRS_PER_TURN);
        }
        return coherence / spectrum.size();
    }
};

// ============================================================================
// COMPLETE G¹⁶ EIGENVALUE SYSTEM
// ============================================================================

class G16EigenSystem {
private:
    Full44x44Matrix matrix_;
    CouplingParams params_;

public:
    explicit G16EigenSystem(const CouplingParams& p = {}) : matrix_(p), params_(p) {}

    // Get full spectral analysis
    struct Analysis {
        Full44x44Matrix::SpectralSummary summary;
        std::vector<double> spectrum;
        std::vector<int> density_of_states;
        double band_gap;
        std::vector<int> zero_modes;
        double spectral_entropy;
        double phase_coherence;
        double trace;
    };

    Analysis analyze() const {
        Analysis a;
        a.summary = matrix_.spectral_summary();
        a.spectrum = matrix_.approximate_spectrum();
        a.density_of_states = matrix_.density_of_states();
        a.band_gap = EigenvalueAnalysis::band_gap(a.spectrum);
        a.zero_modes = EigenvalueAnalysis::find_zero_modes(a.spectrum);
        a.spectral_entropy = EigenvalueAnalysis::spectral_entropy(a.spectrum);
        a.phase_coherence = EigenvalueAnalysis::phase_coherence_10_6(a.spectrum);
        a.trace = matrix_.trace();
        return a;
    }

    // Update coupling parameters and rebuild
    void update_params(const CouplingParams& p) {
        params_ = p;
        matrix_ = Full44x44Matrix(p);
    }

    const Full44x44Matrix& matrix() const { return matrix_; }
    const CouplingParams& params() const { return params_; }
};

} // namespace agi::dna