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

// example_linalg.cpp
// Demonstrates linear algebra features from the calx 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 calx;
using namespace calx::algorithms;

// ---------------------------------------------------------------------------
// Snippet 1: Solve a linear system  Ax = b  (LU + unified interface)
// ---------------------------------------------------------------------------
void demo_linear_solve() {
    std::cout << "=== demo_linear_solve ===" << std::endl;

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

    // LU solver
    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 with QR
    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
// ---------------------------------------------------------------------------
void demo_eigen_symmetric() {
    std::cout << "=== demo_eigen_symmetric ===" << std::endl;

    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)

    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)
// ---------------------------------------------------------------------------
void demo_expm() {
    std::cout << "=== demo_expm ===" << std::endl;

    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)
// ---------------------------------------------------------------------------
void demo_iterative_solver() {
    std::cout << "=== demo_iterative_solver ===" << std::endl;

    // Build a small SPD matrix: A = diag-dominant symmetric
    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
    Vector<double> b(n, 1.0);

    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;
}