Coxeter_triangulation.h

/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Siargey Kachanovich
 *
 *    Copyright (C) 2019 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */
 
#ifndef COXETER_TRIANGULATION_H_
#define COXETER_TRIANGULATION_H_
 
#include <vector>
#include <cmath>  // for std::sqrt
 
#include <boost/range/iterator_range.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
 
#include <Eigen/Eigenvalues>
#include <Eigen/Sparse>
#include <Eigen/SVD>
 
#include <gudhi/Freudenthal_triangulation.h>
#include <gudhi/Permutahedral_representation.h>
 
namespace Gudhi {
 
namespace coxeter_triangulation {
 
template <class Permutahedral_representation_ =
              Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
class Coxeter_triangulation : public Freudenthal_triangulation<Permutahedral_representation_> {
  using Matrix = Eigen::MatrixXd;
 
  Matrix root_matrix(unsigned d) {
    Matrix cartan(Matrix::Identity(d, d));
    for (unsigned i = 1; i < d; i++) {
      cartan(i - 1, i) = -0.5;
      cartan(i, i - 1) = -0.5;
    }
    Eigen::SelfAdjointEigenSolver<Matrix> saes(cartan);
    Eigen::VectorXd sqrt_diag(d);
    for (unsigned i = 0; i < d; ++i) sqrt_diag(i) = std::sqrt(saes.eigenvalues()[i]);
 
    Matrix lower(Matrix::Ones(d, d));
    lower = lower.triangularView<Eigen::Lower>();
 
    Matrix result = (lower * saes.eigenvectors() * sqrt_diag.asDiagonal()).inverse();
    return result;
  }
 
 public:
  Coxeter_triangulation(std::size_t dimension)
      : Freudenthal_triangulation<Permutahedral_representation_>(dimension, root_matrix(dimension)) {}
};
 
}  // namespace coxeter_triangulation
 
}  // namespace Gudhi
 
#endif