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

// example_matrix_demo.cpp
// Matrix demo: showcasing Matrix capabilities in 3 scenarios
//
// Build:
//   cl /std:c++latest /EHsc /O2 /MD /utf-8 /I<sangi>/include examples/example_matrix_demo.cpp /MACHINE:X64

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

using namespace sangi;
using namespace sangi::algorithms;

// ============================================================================
// Demo 1: Hilbert Matrix -- Exact vs Floating-Point
//   -- The Hilbert matrix is notoriously ill-conditioned.
//      Rational arithmetic gives the exact inverse; double accumulates error.
// ============================================================================
static void demo_hilbert() {
    std::cout << "============================================================\n";
    std::cout << " Demo 1: Hilbert Matrix -- Exact vs Floating-Point\n";
    std::cout << "============================================================\n\n";

    const int n = 5;

    // Build Hilbert matrix: H(i,j) = 1/(i+j+1)
    Matrix<Rational> Hr(n, n);
    Matrix<double> Hd(n, n);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            Hr(i, j) = Rational(1, i + j + 1);
            Hd(i, j) = 1.0 / (i + j + 1);
        }
    }

    std::cout << "  H(5x5) Hilbert matrix:\n";
    for (int i = 0; i < n; ++i) {
        std::cout << "    ";
        for (int j = 0; j < n; ++j)
            std::cout << std::setw(8) << Hr(i, j) << " ";
        std::cout << "\n";
    }

    // Exact inverse (Rational)
    auto Hr_inv = Hr ^ (-1);
    std::cout << "\n  Exact inverse H^{-1} (Rational):\n";
    for (int i = 0; i < n; ++i) {
        std::cout << "    ";
        for (int j = 0; j < n; ++j)
            std::cout << std::setw(8) << Hr_inv(i, j) << " ";
        std::cout << "\n";
    }

    // Verify: H * H^{-1} = I
    auto product = Hr * Hr_inv;
    bool is_identity = true;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            if (product(i, j) != Rational(i == j ? 1 : 0))
                is_identity = false;
    std::cout << "\n  H * H^{-1} = I? " << (is_identity ? "YES (exact)" : "NO") << "\n";

    // Double inverse -- check error
    auto Hd_inv = Hd ^ (-1);
    auto prod_d = Hd * Hd_inv;
    double max_err = 0;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            max_err = std::max(max_err, std::abs(prod_d(i, j) - (i == j ? 1.0 : 0.0)));
    std::cout << "  double: max|H*H^{-1} - I| = " << std::scientific
              << max_err << " (rounding error)\n\n";
}

// ============================================================================
// Demo 2: Image Compression with SVD
//   -- A synthetic 8x8 "image" compressed to rank-2 approximation
// ============================================================================
static void demo_svd_compression() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: Image Compression with SVD\n";
    std::cout << "============================================================\n\n";

    // Create an 8x8 matrix with a pattern (gradient + checkerboard)
    const int n = 8;
    Matrix<double> img(n, n);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            img(i, j) = 0.5 * (i + j) + 2.0 * ((i + j) % 2);

    std::cout << "  Original 8x8 image (gradient + checkerboard):\n";
    for (int i = 0; i < n; ++i) {
        std::cout << "    ";
        for (int j = 0; j < n; ++j)
            std::cout << std::fixed << std::setprecision(1) << std::setw(5) << img(i, j);
        std::cout << "\n";
    }

    // SVD
    auto [U, sigma, Vt] = svd_decomposition(img);

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

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

    std::cout << "\n  Rank-2 approximation:\n";
    for (int i = 0; i < n; ++i) {
        std::cout << "    ";
        for (int j = 0; j < n; ++j)
            std::cout << std::fixed << std::setprecision(1) << std::setw(5) << approx(i, j);
        std::cout << "\n";
    }

    // Compression ratio
    int original = n * n;
    int compressed = k * (n + 1 + n);  // k columns of U + k singular values + k rows of Vt
    std::cout << "\n  Storage: " << original << " -> " << compressed
              << " values (ratio " << std::setprecision(1)
              << (100.0 * compressed / original) << "%)\n";

    // Frobenius error
    double err = 0;
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            err += (img(i, j) - approx(i, j)) * (img(i, j) - approx(i, j));
    err = std::sqrt(err);
    std::cout << "  Frobenius error: " << std::setprecision(4) << err << "\n\n";
}

// ============================================================================
// Demo 3: LaTeX Output -- Publish-Ready Matrices
//   -- Generate LaTeX code for matrices and vectors
// ============================================================================
static void demo_latex() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: LaTeX Output\n";
    std::cout << "============================================================\n\n";

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

    std::cout << "  Matrix A (pmatrix):\n" << toLatex(A) << "\n\n";
    std::cout << "  Matrix A (bmatrix):\n" << toLatex(A, "bmatrix") << "\n\n";
    std::cout << "  Vector b:\n" << toLatex(b) << "\n\n";

    // Rational matrix
    Matrix<Rational> R({{Rational(1,2), Rational(1,3)},
                        {Rational(1,4), Rational(1,5)}});
    std::cout << "  Rational matrix:\n" << toLatex(R) << "\n\n";
}

// ============================================================================
int main() {
    std::cout << "###########################################################\n";
    std::cout << "#  sangi Matrix Demo                                       #\n";
    std::cout << "###########################################################\n\n";

    demo_hilbert();
    demo_svd_compression();
    demo_latex();

    std::cout << "=== All demos completed ===\n";
    return 0;
}