eigs_algorithms.hpp

Go to the documentation of this file.

#ifndef eigs_algs_H
#define eigs_algs_H
 
#include "../typedefs.hpp"
#include "Eigen/Eigenvalues"
#include "Eigen/SVD"
#include "../PersistentLaplacians.hpp"
namespace PersistentLaplacians{
 
    // wrapper classes for all algorithms to compute eigenvalues.
    // just need to implement at least one of the following signatures (or both if wrapping with python):
 
    // spectra_vec eigenvalues(DenseMatrix_PL &L);
    // std::pair<spectra_vec,DenseMatrix_spectra_PL> eigenpairs(DenseMatrix_PL &L);
 
    // The "eigenvalues" function is necessary if you call the "spectra" function of a PersistentLaplacian, 
    // The "eigenpairs" function is necessary if you call the "eigenpairs" function of a PersistentLaplacian.
    // The python wrapping code requires both, but for a c++ only version you can just use one.
 
 
    // Wrapper for Eigen's SelfAdjointEigenSolver
    class selfadjoint{
        public:
            spectra_vec eigenvalues(DenseMatrix_PL &L){
                Eigen::SelfAdjointEigenSolver<DenseMatrix_PL> es = Eigen::SelfAdjointEigenSolver<DenseMatrix_PL>(L, Eigen::EigenvaluesOnly);
                spectra_vec eigs = es.eigenvalues();
                round_zeros(eigs, 1e-3);
                return eigs;
            }
            std::pair<spectra_vec,DenseMatrix_spectra_PL> eigenpairs(DenseMatrix_PL &L){
                Eigen::SelfAdjointEigenSolver<DenseMatrix_PL> es = Eigen::SelfAdjointEigenSolver<DenseMatrix_PL>(L);
                spectra_vec eigs = es.eigenvalues();
                DenseMatrix_spectra_PL eigvs = es.eigenvectors();
                round_zeros(eigs, 1e-3);
                return std::pair<spectra_vec,DenseMatrix_spectra_PL>(eigs, eigvs);
            }
    };
 
    // Wrapper for Eigen's standard EigenSolver (does not assume self-adjoint)
    class eigensolver{
        public:
            spectra_vec eigenvalues(DenseMatrix_PL &L){
                Eigen::EigenSolver<DenseMatrix_PL> es = Eigen::EigenSolver<DenseMatrix_PL>(L, Eigen::EigenvaluesOnly);
                spectra_vec eigs = es.eigenvalues().real();
                round_zeros(eigs, 1e-3);
                std::sort(eigs.begin(), eigs.end()); // Eigen::Eigensolver returns unsorted
                return eigs;
            }
            // TODO: this is placeholder
            std::pair<spectra_vec,DenseMatrix_spectra_PL> eigenpairs(DenseMatrix_PL &L){
                Eigen::SelfAdjointEigenSolver<DenseMatrix_PL> es = Eigen::SelfAdjointEigenSolver<DenseMatrix_PL>(L);
                spectra_vec eigs = es.eigenvalues();
                DenseMatrix_spectra_PL eigvs = es.eigenvectors();
                round_zeros(eigs, 1e-3);
                return std::pair<spectra_vec,DenseMatrix_spectra_PL>(eigs, eigvs);
            }
    };
 
    // // Wrapper for Eigen's standard EigenSolver (does not assume self-adjoint) that gives eigenvectors too
    // class eigensolver_v{
    //     public:
    //         std::pair<spectra_vec,DenseMatrix_spectra_PL> operator()(DenseMatrix_PL &L){
    //             Eigen::EigenSolver<DenseMatrix_PL> es = Eigen::EigenSolver<DenseMatrix_PL>(L); // Eigen::Eigensolver returns unsorted
    //             spectra_vec eigs = es.eigenvalues().real();
    //             DenseMatrix_spectra_PL eigvs = es.eigenvectors().real();
                
    //             // sort eigenvectors and eigenvalues together - wrap in custom struct
    //             int num_eigs = eigs.size();
    //             std::vector<Eigenpair> eigenpairs;
    //             eigenpairs.reserve(num_eigs);
 
    //             for (int i = 0; i < num_eigs; i++){
    //                 Eigenpair p = {eigs[i], eigvs.col(i)};
    //                 eigenpairs.push_back(p);
    //             }
    //             round_zeros(eigs, 1e-3);
 
    //             std::sort(eigenpairs.begin(), eigenpairs.end(),
    //                         [](const auto& a, const auto& b) {return a.eig < b.eig;}); // custom comparator for sorting
                
    //             spectra_vec eigs_sorted(num_eigs);
    //             DenseMatrix_spectra_PL eigvs_sorted(eigvs.rows(),eigvs.cols());
                
    //             for (int i = 0; i < num_eigs; i++){
    //                 eigs_sorted[i] = eigenpairs[i].eig;
    //                 eigvs_sorted.col(i) = eigenpairs[i].eigv;
    //             }
 
 
    //             return std::pair<spectra_vec,DenseMatrix_spectra_PL>(eigs_sorted, eigvs_sorted);
    //         }
    //     private:
    //         struct Eigenpair{
    //             spectra_type eig;// one eigenvalue
    //             DenseMatrix_spectra_PL eigv; // one eigenvector
    //         };
 
    // };
 
    // Wrapper for Eigen's BDCSVD singular value decomposition 
    class bdcsvd{
        public:
            spectra_vec eigenvalues(DenseMatrix_PL &L){
                  Eigen::BDCSVD<DenseMatrix_PL> bdcsvd(L.rows(), L.cols());
                bdcsvd.compute(L);
                spectra_vec eigs = bdcsvd.singularValues();
                round_zeros(eigs, 1e-3);
                std::sort(eigs.begin(), eigs.end());
                return eigs;
            }
            // TODO: this is placeholder
            std::pair<spectra_vec,DenseMatrix_spectra_PL> eigenpairs(DenseMatrix_PL &L){
                Eigen::SelfAdjointEigenSolver<DenseMatrix_PL> es = Eigen::SelfAdjointEigenSolver<DenseMatrix_PL>(L);
                spectra_vec eigs = es.eigenvalues();
                DenseMatrix_spectra_PL eigvs = es.eigenvectors();
                round_zeros(eigs, 1e-3);
                return std::pair<spectra_vec,DenseMatrix_spectra_PL>(eigs, eigvs);
            }
    };
 
 
}
 
#endif