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

// example_linalg.cpp
// Demonstrates linear algebra features from the sangi API reference.
//
// Build:
//   cd build && cmake --build . --config Release --target example-linalg

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/linalg/solvers.hpp>
#include <math/linalg/eigenvalues.hpp>
#include <math/linalg/matrix_functions.hpp>
#include <iostream>
#include <cmath>

using namespace sangi;
using namespace sangi::algorithms;

// ---------------------------------------------------------------------------
// Snippet 1: Solve a linear system  Ax = b  (LU + unified interface)
//
// Matrix A and right-hand side b form a classic test system (solution x = {2, 3, -1}).
// LU decomposition is O(n^3) and applies to general nonsingular matrices.
// QR decomposition also applies to overdetermined systems and is numerically more stable.
// ---------------------------------------------------------------------------
void demo_linear_solve() {
    std::cout << "=== demo_linear_solve ===" << std::endl;

    // Classic test matrix (solution: x = {2, 3, -1})
    Matrix<double> A({{2, 1, -1},
                      {-3, -1, 2},
                      {-2, 1, 2}});
    Vector<double> b({8, -11, -3});

    // LU decomposition: numerically stable with partial pivoting
    auto x = lu_solve(A, b);
    std::cout << "x (LU)  = [";
    for (std::size_t i = 0; i < x.size(); ++i)
        std::cout << " " << x[i];
    std::cout << " ]" << std::endl;
    // x should be approx {2, 3, -1}

    // Unified interface (SolverType::QR): also usable for overdetermined systems
    auto x2 = solve(A, b, SolverType::QR);
    std::cout << "x (QR)  = [";
    for (std::size_t i = 0; i < x2.size(); ++i)
        std::cout << " " << x2[i];
    std::cout << " ]" << std::endl;

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

// ---------------------------------------------------------------------------
// Snippet 2: Eigenvalue decomposition of a symmetric matrix
//
// Eigenvalue decomposition of a real symmetric matrix: S = V D V^T (D is the diagonal
// eigenvalue matrix, V is orthogonal).
// sangi::eigen_symmetric returns them in ascending order, compatible with LAPACK dsyev.
// ---------------------------------------------------------------------------
void demo_eigen_symmetric() {
    std::cout << "=== demo_eigen_symmetric ===" << std::endl;

    // 3×3 real symmetric matrix (trace = 9, det ≈ 14.62)
    Matrix<double> S({{4, 1, 1},
                      {1, 3, 0},
                      {1, 0, 2}});

    auto [eigenvalues, eigenvectors] = eigen_symmetric(S);
    // eigenvalues should be approx {1.38, 3.0, 4.62} (ascending order)

    std::cout << "eigenvalues = [";
    for (std::size_t i = 0; i < eigenvalues.size(); ++i)
        std::cout << " " << eigenvalues[i];
    std::cout << " ]" << std::endl;

    std::cout << "eigenvectors (columns):" << std::endl;
    for (std::size_t i = 0; i < eigenvectors.rows(); ++i) {
        std::cout << "  [";
        for (std::size_t j = 0; j < eigenvectors.cols(); ++j)
            std::cout << " " << eigenvectors(i, j);
        std::cout << " ]" << std::endl;
    }

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

// ---------------------------------------------------------------------------
// Snippet 3: Matrix exponential  (rotation matrix from skew-symmetric)
//
// For a skew-symmetric matrix A, expm(A) is a rotation matrix.
// A = [[0,-1],[1,0]] is the 2D rotation generator for θ=1rad.
// expm(A) = [[cos1, -sin1],[sin1, cos1]] holds.
// sangi::expm uses the Padé approximation + scaling & squaring method.
// ---------------------------------------------------------------------------
void demo_expm() {
    std::cout << "=== demo_expm ===" << std::endl;

    // 2D rotation generator (skew-symmetric): expm(A) is the rotation matrix for θ=1rad
    Matrix<double> A({{0, -1},
                      {1,  0}});

    auto E = expm(A);
    // E should be approx {{ cos(1), -sin(1) },
    //                      { sin(1),  cos(1) }}

    std::cout << "expm(A) =" << std::endl;
    for (std::size_t i = 0; i < E.rows(); ++i) {
        std::cout << "  [";
        for (std::size_t j = 0; j < E.cols(); ++j)
            std::cout << " " << E(i, j);
        std::cout << " ]" << std::endl;
    }
    std::cout << "expected: [[ " << std::cos(1.0) << " " << -std::sin(1.0)
              << " ] [ " << std::sin(1.0) << " " << std::cos(1.0) << " ]]"
              << std::endl;

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

// ---------------------------------------------------------------------------
// Snippet 4: Iterative solver (CG on a symmetric positive-definite matrix)
//
// Conjugate gradient (CG): an iterative method O(n * iter) suited to large sparse SPD systems.
// Convergence guarantee: if A is SPD, the exact solution is reached in at most n iterations
// (within rounding error).
// A tolerance of 1e-12 is about 5000 times the double machine epsilon (~2.2e-16).
// ---------------------------------------------------------------------------
void demo_iterative_solver() {
    std::cout << "=== demo_iterative_solver ===" << std::endl;

    // 1D Poisson difference matrix (tridiagonal SPD): A_{ii}=10, A_{i,i±1}=-1
    // Diagonally dominant, so the condition number O(n^2 / diagonal_ratio) is favorable
    const int n = 5;
    Matrix<double> A(n, n, 0.0);
    for (int i = 0; i < n; ++i) {
        A(i, i) = 10.0;
        if (i > 0)     { A(i, i - 1) = -1.0; A(i - 1, i) = -1.0; }
    }

    // Right-hand side: b = (1, 1, ..., 1)
    Vector<double> b(n, 1.0);

    // tol=1e-12, max_iter=5000 (since n=5, it should converge within 5 iterations)
    auto x = conjugate_gradient(A, b, 1e-12, 5000);

    std::cout << "x (CG) = [";
    for (std::size_t i = 0; i < x.size(); ++i)
        std::cout << " " << x[i];
    std::cout << " ]" << std::endl;

    // Verify residual
    auto r = A * x;
    double res = 0.0;
    for (std::size_t i = 0; i < b.size(); ++i)
        res += (r[i] - b[i]) * (r[i] - b[i]);
    std::cout << "||Ax - b|| = " << std::sqrt(res) << std::endl;

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

// ---------------------------------------------------------------------------
int main() {
    demo_linear_solve();
    demo_eigen_symmetric();
    demo_expm();
    demo_iterative_solver();
    return 0;
}