// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_matrix.cpp
// Matrix class demos: basic operations, Gaussian elimination with pivot
// strategies, SVD low-rank approximation, and exact Rational arithmetic.
//
// Build:
//   cd build && cmake --build . --config Release --target example-matrix

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/mp/Rational.hpp>
#include <math/linalg/decomposition.hpp>
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace calx;
using namespace calx::algorithms;

// ---------------------------------------------------------------------------
// Demo 1: Basic matrix operations
// ---------------------------------------------------------------------------
void demo_basic() {
    std::cout << "=== Basic matrix operations ===" << std::endl;

    Matrix<double> A({{1, 2, 3},
                      {4, 5, 6},
                      {7, 8, 10}});

    std::cout << "A:" << std::endl;
    print_matrix(A);

    std::cout << "det(A)     = " << A.determinant() << std::endl;
    std::cout << "trace(A)   = " << A.trace() << std::endl;
    std::cout << "norm_F(A)  = " << A.norm() << std::endl;

    auto At = A.transpose();
    std::cout << "A^T:" << std::endl;
    print_matrix(At);

    auto Ainv = A.inverse();
    std::cout << "A^{-1}:" << std::endl;
    print_matrix(Ainv);

    auto I = A * Ainv;
    std::cout << "A * A^{-1} (should be I):" << std::endl;
    print_matrix(I);

    std::cout << std::endl;
}

// ---------------------------------------------------------------------------
// Demo 2: Gaussian elimination with pivot strategies
// ---------------------------------------------------------------------------
void demo_gauss_pivot() {
    std::cout << "=== Gaussian elimination: pivot strategies ===" << std::endl;

    Matrix<double> A({{2, 1, -1},
                      {-3, -1, 2},
                      {-2, 1, 2}});
    Vector<double> b({8, -11, -3});

    auto x_none    = gaussian_elimination(A, b, PivotStrategy::None);
    auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);
    auto x_full    = gaussian_elimination(A, b, PivotStrategy::Full);

    auto print_vec = [](const char* label, const Vector<double>& v) {
        std::cout << label << " = [";
        for (std::size_t i = 0; i < v.size(); ++i)
            std::cout << " " << v[i];
        std::cout << " ]" << std::endl;
    };

    print_vec("x (no pivot)     ", x_none);
    print_vec("x (partial pivot)", x_partial);
    print_vec("x (full pivot)   ", x_full);

    // Residual
    Vector<double> r = A * x_partial - b;
    std::cout << "||Ax - b|| = " << r.norm() << std::endl;

    std::cout << std::endl;
}

// ---------------------------------------------------------------------------
// Demo 3: SVD low-rank approximation
// ---------------------------------------------------------------------------
void demo_svd_lowrank() {
    std::cout << "=== SVD low-rank approximation ===" << std::endl;

    // 5x4 matrix (approximately rank 2)
    Matrix<double> A({{1.02, 2.01, 3.00, 3.98},
                      {2.01, 3.99, 6.02, 8.01},
                      {0.49, 1.02, 1.48, 1.99},
                      {3.00, 5.98, 9.01, 12.0},
                      {1.50, 3.01, 4.49, 5.99}});

    auto [U, sigma, Vt] = svd_decomposition(A);

    std::cout << "Singular values:" << std::endl;
    for (std::size_t i = 0; i < sigma.size(); ++i)
        std::cout << "  sigma[" << i << "] = " << std::fixed
                  << std::setprecision(6) << sigma[i] << std::endl;

    // Rank-2 approximation
    std::size_t k = 2;
    Vector<double> s_trunc(sigma.size(), 0.0);
    for (std::size_t i = 0; i < k; ++i) s_trunc[i] = sigma[i];
    auto A_approx = reconstruct_matrix_from_svd(U, s_trunc, Vt);

    double err = 0;
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j) {
            double d = A(i, j) - A_approx(i, j);
            err += d * d;
        }
    std::cout << "Rank-" << k << " Frobenius error: "
              << std::scientific << std::setprecision(2) << std::sqrt(err) << std::endl;

    std::cout << std::endl;
}

// ---------------------------------------------------------------------------
// Demo 4: Exact Hilbert matrix inverse with Rational
// ---------------------------------------------------------------------------
void demo_rational_hilbert() {
    std::cout << "=== Exact Hilbert matrix inverse (Rational) ===" << std::endl;

    const int N = 5;
    Matrix<Rational> H(N, N);
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            H(i, j) = Rational(1, i + j + 1);

    std::cout << "H(5):" << std::endl;
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j)
            std::cout << std::setw(8) << H(i, j);
        std::cout << std::endl;
    }

    auto Hinv = lu_inverse(H);

    std::cout << "H^{-1}:" << std::endl;
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j)
            std::cout << std::setw(8) << Hinv(i, j);
        std::cout << std::endl;
    }

    // Verify H * H^{-1} == I exactly
    auto I = H * Hinv;
    bool exact = true;
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            if (I(i, j) != Rational(i == j ? 1 : 0)) exact = false;

    std::cout << "H * H^{-1} == I (exact): " << (exact ? "true" : "false") << std::endl;

    std::cout << std::endl;
}

// ---------------------------------------------------------------------------
// Demo 5: Rational Gaussian elimination with height-minimum pivot
// ---------------------------------------------------------------------------
void demo_rational_gauss() {
    std::cout << "=== Rational Gaussian elimination ===" << std::endl;

    // Solve Hx = b where H is 4x4 Hilbert, b chosen so x = (1,1,1,1)
    const int N = 4;
    Matrix<Rational> H(N, N);
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            H(i, j) = Rational(1, i + j + 1);

    Vector<Rational> b(N);
    for (int i = 0; i < N; ++i) {
        Rational s;
        for (int j = 0; j < N; ++j) s = s + H(i, j);
        b[i] = s;
    }

    auto x = gaussian_elimination(H, b, PivotStrategy::Partial);

    std::cout << "x = [";
    for (int i = 0; i < N; ++i)
        std::cout << " " << x[i];
    std::cout << " ]" << std::endl;
    std::cout << "(height-minimum pivot selection for Rational)" << std::endl;

    std::cout << std::endl;
}

// ---------------------------------------------------------------------------
int main() {
    demo_basic();
    demo_gauss_pivot();
    demo_svd_lowrank();
    demo_rational_hilbert();
    demo_rational_gauss();
    return 0;
}