Loading...
Searching...
No Matches
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
Gudhi::coxeter_triangulation::Freudenthal_triangulation::offset
const Vector & offset() const
Vector that defines the offset of the triangulation.
Definition Freudenthal_triangulation.h:99
Gudhi::coxeter_triangulation::Freudenthal_triangulation::Freudenthal_triangulation
Freudenthal_triangulation(unsigned dimension, const Matrix &matrix, const Vector &offset)
Constructor of the Freudenthal-Kuhn triangulation of a given dimension under an affine transformation...
Definition Freudenthal_triangulation.h:83
Gudhi::coxeter_triangulation::Freudenthal_triangulation::barycenter
Eigen::VectorXd barycenter(const Simplex_handle &simplex, double scale=1) const
Returns the Cartesian coordinates of the barycenter of a given simplex.
Definition Freudenthal_triangulation.h:198
Gudhi::coxeter_triangulation::Freudenthal_triangulation::locate_point
Simplex_handle locate_point(const Point_d &point, double scale=1) const
Returns the permutahedral representation of the simplex in the triangulation that contains a given qu...
Definition Freudenthal_triangulation.h:135
Gudhi::coxeter_triangulation::Freudenthal_triangulation::Vertex_handle
typename Permutahedral_representation_::Vertex Vertex_handle
Type of the vertices in the triangulation.
Definition Freudenthal_triangulation.h:55
Gudhi::coxeter_triangulation::Freudenthal_triangulation::Freudenthal_triangulation
Freudenthal_triangulation(std::size_t dimension)
Constructor of the Freudenthal-Kuhn triangulation of a given dimension.
Definition Freudenthal_triangulation.h:60
Gudhi::coxeter_triangulation::Freudenthal_triangulation::Freudenthal_triangulation
Freudenthal_triangulation(std::size_t dimension, const Matrix &matrix)
Constructor of the Freudenthal-Kuhn triangulation of a given dimension under a linear transformation ...
Definition Freudenthal_triangulation.h:71
Gudhi::coxeter_triangulation::Freudenthal_triangulation::Simplex_handle
Permutahedral_representation_ Simplex_handle
Type of the simplices in the triangulation.
Definition Freudenthal_triangulation.h:52
Gudhi::coxeter_triangulation::Freudenthal_triangulation::change_offset
void change_offset(const Eigen::VectorXd &offset)
Change the offset vector to a given value.
Definition Freudenthal_triangulation.h:113
Gudhi::coxeter_triangulation::Freudenthal_triangulation::matrix
const Matrix & matrix() const
Matrix that defines the linear transformation of the triangulation.
Definition Freudenthal_triangulation.h:96
Gudhi::coxeter_triangulation::Freudenthal_triangulation::cartesian_coordinates
Eigen::VectorXd cartesian_coordinates(const Vertex_handle &vertex, double scale=1) const
Returns the Cartesian coordinates of the given vertex.
Definition Freudenthal_triangulation.h:184
Gudhi::coxeter_triangulation::Freudenthal_triangulation::dimension
unsigned dimension() const
Dimension of the triangulation.
Definition Freudenthal_triangulation.h:93
Gudhi::coxeter_triangulation::Freudenthal_triangulation::change_matrix
void change_matrix(const Eigen::MatrixXd &matrix)
Change the linear transformation matrix to a given value.
Definition Freudenthal_triangulation.h:104
Gudhi::coxeter_triangulation::Permutahedral_representation
A class that stores the permutahedral representation of a simplex in a Coxeter triangulation or a Fre...
Definition Permutahedral_representation.h:38
Gudhi
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14
Gudhi::coxeter_triangulation::Ordered_set_partition
Class that represents an ordered set partition of a set {0,...,n-1} in k parts as a pair of an unorde...
Definition Ordered_set_partition_iterator.h:32