#include "graph.hpp"
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

using MatrixType =
    Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;

// A: Adjacency matrix
//    lambda_max <= d_max
//
// L: Laplacian
//    L = D - A
///
// P: Random walk matrix
//    lambda_max = 1

class GraphSpectrum {
  public:
    GraphSpectrum(const Graph& g) : graph(g) {}
    ~GraphSpectrum() {}

    std::vector<float> computeAdjacencySpectrum() const {
        // matrix stored as std::vector<std::vector<bool>>
        auto& badj = graph.getBooleanAdj();

        // Convert it to MatrixType
        auto n = badj.size();
        MatrixType m(n, n);
        for (auto i = 0u; i < n; ++i)
            for (auto j = 0u; j < n; ++j)
                m(i, j) = badj[i][j] ? 1.0f : 0.0f;

        return getEigenvalues_(m);
    }

    std::vector<float> computeLaplacianSpectrum() const {
        // matrix stored as std::vector<std::vector<bool>>
        auto& badj = graph.getBooleanAdj();
        auto& adj = graph.getAdj();

        // - A
        auto n = badj.size();
        MatrixType m(n, n);
        for (auto i = 0u; i < n; ++i)
            for (auto j = 0u; j < n; ++j)
                m(i, j) = badj[i][j] ? -1.0f : 0.0f;

        // + D
        for (auto i = 0u; i < n; ++i)
            m(i, i) = float(adj[i].size());

        return getEigenvalues_(m);
    }

  private:
    const Graph& graph;

    std::vector<float> getEigenvalues_(const MatrixType& m) const {
        Eigen::SelfAdjointEigenSolver<MatrixType> es(
            m, Eigen::DecompositionOptions::EigenvaluesOnly);
        auto ev = es.eigenvalues();
        return std::vector<float>(ev.data(), ev.data() + ev.rows() * ev.cols());
    }
};