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

/**
 * @file sparse_ordering.hpp
 * @brief Fill-reducing ordering for sparse direct solvers (minimum-degree / AMD)
 *
 * Sparse LU / Cholesky / QR factorizations suffer from explosive fill-in unless
 * the elimination order is chosen carefully; with the natural order they become
 * catastrophically slow and memory-hungry on real problems such as FEM, circuits,
 * and graph Laplacians. This header computes a symmetric, fill-reducing order
 * using the minimum-degree heuristic on the elimination graph (the basis of the
 * AMD family).
 *
 *   SparseOrdering::Natural — no ordering (identity permutation)
 *   SparseOrdering::AMD     — minimum-degree ordering
 *
 * The returned permutation perm satisfies perm[k] = the original index eliminated
 * at step k. Because the symmetric elimination graph is built from pattern(A + Aᵀ),
 * it works for both symmetric (Cholesky) and unsymmetric (LU) cases. For LU it is
 * applied as a symmetric permutation P A Pᵀ together with partial pivoting.
 *
 * Note: this uses exact minimum-degree (an O(n) scan at each step). It targets
 * medium-sized problems (up to a few thousand dimensions); for large problems
 * Natural can be selected instead. Fill-reduction quality is at least as good as
 * AMD (AMD approximates for the sake of speed).
 */

#ifndef SANGI_SPARSE_ORDERING_HPP
#define SANGI_SPARSE_ORDERING_HPP

#include <vector>
#include <set>
#include <cstddef>
#include <numeric>

#include "../core/sparse_matrix.hpp"

namespace sangi {

/// Fill-reducing ordering options for sparse direct solvers
enum class SparseOrdering {
    Natural,  ///< no ordering (identity permutation)
    AMD       ///< minimum-degree ordering (AMD family)
};

namespace detail {
namespace ordering {

/**
 * @brief (Exact) minimum-degree ordering on the symmetric elimination graph
 *
 * @param n   number of vertices
 * @param adj symmetric adjacency sets (no self-loops). Destroyed by the call.
 * @return perm: perm[k] = vertex (original index) eliminated at step k
 */
inline std::vector<std::size_t> minimum_degree(
    std::size_t n, std::vector<std::set<std::size_t>>& adj)
{
    std::vector<std::size_t> perm;
    perm.reserve(n);
    std::vector<char> elim(n, 0);

    std::vector<std::size_t> nbrs;
    nbrs.reserve(n);

    for (std::size_t step = 0; step < n; ++step) {
        // Pick the uneliminated vertex of minimum degree (ties broken by smallest index)
        std::size_t v = n;
        std::size_t best_deg = static_cast<std::size_t>(-1);
        for (std::size_t u = 0; u < n; ++u) {
            if (elim[u]) continue;
            std::size_t deg = adj[u].size();
            if (deg < best_deg) { best_deg = deg; v = u; if (deg == 0) break; }
        }

        perm.push_back(v);
        elim[v] = 1;

        // Collect the uneliminated neighbors of v
        nbrs.clear();
        for (std::size_t u : adj[v])
            if (!elim[u]) nbrs.push_back(u);

        // Remove v from each neighbor
        for (std::size_t u : nbrs) adj[u].erase(v);
        adj[v].clear();

        // Turn the neighbors into a clique (models the fill-in)
        for (std::size_t a = 0; a < nbrs.size(); ++a) {
            for (std::size_t b = a + 1; b < nbrs.size(); ++b) {
                adj[nbrs[a]].insert(nbrs[b]);
                adj[nbrs[b]].insert(nbrs[a]);
            }
        }
    }
    return perm;
}

/// Build symmetric adjacency sets from pattern(A + Aᵀ) of a SparseMatrix
template<typename T, typename IndexType>
std::vector<std::set<std::size_t>> build_symmetric_adjacency(
    const SparseMatrix<T, IndexType>& A)
{
    const auto n = static_cast<std::size_t>(A.rows());
    std::vector<std::set<std::size_t>> adj(n);

    const auto& col_ptr = A.csc_col_ptr();
    const auto& row_ind = A.csc_row_indices();
    for (std::size_t j = 0; j < n; ++j) {
        const auto start = static_cast<std::size_t>(col_ptr[j]);
        const auto end   = static_cast<std::size_t>(col_ptr[j + 1]);
        for (std::size_t p = start; p < end; ++p) {
            const auto i = static_cast<std::size_t>(row_ind[p]);
            if (i != j) { adj[i].insert(j); adj[j].insert(i); }
        }
    }
    return adj;
}

} // namespace ordering
} // namespace detail

namespace algorithms {

/**
 * @brief Fill-reducing ordering (minimum-degree / AMD family)
 *
 * Builds the symmetric elimination graph from pattern(A + Aᵀ) and returns a
 * minimum-degree ordering.
 *
 * @return perm: perm[k] = the original index eliminated at step k.
 *         Empty when n == 0.
 */
template<typename T, typename IndexType>
std::vector<std::size_t> amd_ordering(const SparseMatrix<T, IndexType>& A)
{
    const auto n = static_cast<std::size_t>(A.rows());
    if (n == 0) return {};
    auto adj = detail::ordering::build_symmetric_adjacency(A);
    return detail::ordering::minimum_degree(n, adj);
}

/// Return the identity permutation (Natural ordering)
inline std::vector<std::size_t> natural_ordering(std::size_t n) {
    std::vector<std::size_t> perm(n);
    std::iota(perm.begin(), perm.end(), std::size_t(0));
    return perm;
}

/**
 * @brief Build the symmetric permutation B = P A Pᵀ (B(k,l) = A(perm[k], perm[l]))
 *
 * Preprocessing step for applying a fill-reducing order in Cholesky / symmetric
 * factorizations.
 */
template<typename T, typename IndexType>
SparseMatrix<T, IndexType> permute_symmetric(
    const SparseMatrix<T, IndexType>& A, const std::vector<std::size_t>& perm)
{
    const auto n = static_cast<std::size_t>(A.rows());
    SparseMatrix<T, IndexType> B(static_cast<IndexType>(n),
                                 static_cast<IndexType>(n),
                                 SparseStorageFormat::CSC);
    std::vector<std::size_t> inv(n);
    for (std::size_t k = 0; k < n; ++k) inv[perm[k]] = k;

    const auto& av = A.csc_values();
    const auto& ar = A.csc_row_indices();
    const auto& ac = A.csc_col_ptr();
    for (std::size_t j = 0; j < n; ++j) {
        for (auto p = static_cast<std::size_t>(ac[j]);
             p < static_cast<std::size_t>(ac[j + 1]); ++p) {
            const auto i = static_cast<std::size_t>(ar[p]);
            B.set_coeff(static_cast<IndexType>(inv[i]),
                        static_cast<IndexType>(inv[j]), av[p]);
        }
    }
    return B;
}

} // namespace algorithms
} // namespace sangi

#endif // SANGI_SPARSE_ORDERING_HPP