Freudenthal_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 FREUDENTHAL_TRIANGULATION_H_
#define FREUDENTHAL_TRIANGULATION_H_
 
#include <vector>
#include <algorithm>  // for std::sort
#include <cmath>      // for std::floor
#include <numeric>    // for std::iota
#include <cstdlib>    // for std::size_t
 
#include <Eigen/Eigenvalues>
#include <Eigen/SVD>
 
#include <gudhi/Permutahedral_representation.h>
#include <gudhi/Debug_utils.h>  // for GUDHI_CHECK
 
namespace Gudhi {
 
namespace coxeter_triangulation {
 
template <class Permutahedral_representation_ =
              Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
class Freudenthal_triangulation {
  using Matrix = Eigen::MatrixXd;
  using Vector = Eigen::VectorXd;
 
 public:
  using Simplex_handle = Permutahedral_representation_;
 
  using Vertex_handle = typename Permutahedral_representation_::Vertex;
 
  Freudenthal_triangulation(std::size_t dimension)
      : Freudenthal_triangulation(dimension, Matrix::Identity(dimension, dimension), Vector::Zero(dimension)) {
    is_freudenthal_ = true;
  }
 
  Freudenthal_triangulation(std::size_t dimension, const Matrix& matrix)
      : Freudenthal_triangulation(dimension, matrix, Vector::Zero(dimension)) {}
 
  Freudenthal_triangulation(unsigned dimension, const Matrix& matrix, const Vector& offset)
      : dimension_(dimension),
        matrix_(matrix),
        offset_(offset),
        colpivhouseholderqr_(matrix_.colPivHouseholderQr()),
        is_freudenthal_(false) {
    GUDHI_CHECK(dimension == offset_.size(), std::invalid_argument("Offset must be of size 'dimension'"));
  }
 
  unsigned dimension() const { return dimension_; }
 
  const Matrix& matrix() const { return matrix_; }
 
  const Vector& offset() const { return offset_; }
 
  void change_matrix(const Eigen::MatrixXd& matrix) {
    matrix_ = matrix;
    colpivhouseholderqr_ = matrix.colPivHouseholderQr();
    is_freudenthal_ = false;
  }
 
  void change_offset(const Eigen::VectorXd& offset) {
    offset_ = offset;
    is_freudenthal_ = false;
  }
 
  template <class Point_d>
  Simplex_handle locate_point(const Point_d& point, double scale = 1) const {
    using Ordered_set_partition = typename Simplex_handle::OrderedSetPartition;
    using Part = typename Ordered_set_partition::value_type;
    unsigned d = point.size();
    GUDHI_CHECK(d == dimension_,
                std::invalid_argument("The point must be of the same dimension as the triangulation"));
    double error = 1e-9;
    Simplex_handle output;
    std::vector<double> z;
    if (is_freudenthal_) {
      for (std::size_t i = 0; i < d; i++) {
        double x_i = scale * point[i];
        int y_i = std::floor(x_i);
        output.vertex().push_back(y_i);
        z.push_back(x_i - y_i);
      }
    } else {
      Eigen::VectorXd p_vect(d);
      for (std::size_t i = 0; i < d; i++) p_vect(i) = point[i];
      Eigen::VectorXd x_vect = colpivhouseholderqr_.solve(p_vect - offset_);
      for (std::size_t i = 0; i < d; i++) {
        double x_i = scale * x_vect(i);
        int y_i = std::floor(x_i);
        output.vertex().push_back(y_i);
        z.push_back(x_i - y_i);
      }
    }
    z.push_back(0);
    Part indices(d + 1);
    std::iota(indices.begin(), indices.end(), 0);
    std::sort(indices.begin(), indices.end(), [&z](std::size_t i1, std::size_t i2) { return z[i1] > z[i2]; });
 
    output.partition().push_back(Part(1, indices[0]));
    for (std::size_t i = 1; i <= d; ++i)
      if (z[indices[i - 1]] > z[indices[i]] + error)
        output.partition().push_back(Part(1, indices[i]));
      else
        output.partition().back().push_back(indices[i]);
    return output;
  }
 
  Eigen::VectorXd cartesian_coordinates(const Vertex_handle& vertex, double scale = 1) const {
    Eigen::VectorXd v_vect(dimension_);
    for (std::size_t j = 0; j < dimension_; j++) v_vect(j) = vertex[j] / scale;
    return matrix_ * v_vect + offset_;
  }
 
  Eigen::VectorXd barycenter(const Simplex_handle& simplex, double scale = 1) const {
    Eigen::VectorXd res_vector(dimension_);
    res_vector.setZero(dimension_, 1);
    for (auto v : simplex.vertex_range()) {
      res_vector += cartesian_coordinates(v, scale);
    }
    return (1. / (simplex.dimension() + 1)) * res_vector;
  }
 
 protected:
  unsigned dimension_;
  Matrix matrix_;
  Vector offset_;
  Eigen::ColPivHouseholderQR<Matrix> colpivhouseholderqr_;
  bool is_freudenthal_;
};
 
}  // namespace coxeter_triangulation
 
}  // namespace Gudhi
 
#endif