diff options
Diffstat (limited to 'src/Coxeter_triangulation/include')
37 files changed, 5037 insertions, 0 deletions
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h new file mode 100644 index 00000000..de68acb6 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h @@ -0,0 +1,77 @@ +/* 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 { + +/** + * \class Coxeter_triangulation + * \brief A class that stores Coxeter triangulation of type \f$\tilde{A}_d\f$. + * This triangulation has the greatest simplex quality out of all linear transformations + * of the Freudenthal-Kuhn triangulation. + * + * \ingroup coxeter_triangulation + * + * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation. + * Needs to be a model of SimplexInCoxeterTriangulation. + */ +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: + /** \brief Constructor of Coxeter triangulation of a given dimension. + * @param[in] dimension The dimension of the triangulation. + */ + Coxeter_triangulation(std::size_t dimension) + : Freudenthal_triangulation<Permutahedral_representation_>(dimension, root_matrix(dimension)) {} +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h new file mode 100644 index 00000000..de342ecc --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h @@ -0,0 +1,340 @@ +/* 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 CELL_COMPLEX_H_ +#define CELL_COMPLEX_H_ + +#include <Eigen/Dense> + +#include <vector> +#include <map> +#include <utility> // for std::make_pair + +#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES +#include <gudhi/Permutahedral_representation/Simplex_comparator.h> +#include <gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h> // for Hasse_cell + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Cell_complex + * \brief A class that constructs the cell complex from the output provided by the class + * \ref Gudhi::coxeter_triangulation::Manifold_tracing. + * + * The use and interfaces of this cell complex is limited to the \ref coxeter_triangulation implementation. + * + * \tparam Out_simplex_map_ The type of a map from a simplex type that is a + * model of SimplexInCoxeterTriangulation to Eigen::VectorXd. + */ +template <class Out_simplex_map_> +class Cell_complex { + public: + /** \brief Type of a simplex in the ambient triangulation. + * Is a model of the concept SimplexInCoxeterTriangulation. + */ + using Simplex_handle = typename Out_simplex_map_::key_type; + /** \brief Type of a cell in the cell complex. + * Always is Gudhi::Hasse_cell from the Hasse diagram module. + * The additional information is the boolean that is true if and only if the cell lies + * on the boundary. + */ + using Hasse_cell = Gudhi::Hasse_diagram::Hasse_diagram_cell<int, double, bool>; + /** \brief Type of a map from permutahedral representations of simplices in the + * ambient triangulation to the corresponding cells in the cell complex of some + * specific dimension. + */ + using Simplex_cell_map = std::map<Simplex_handle, Hasse_cell*, Simplex_comparator<Simplex_handle> >; + /** \brief Type of a vector of maps from permutahedral representations of simplices in the + * ambient triangulation to the corresponding cells in the cell complex of various dimensions. + */ + using Simplex_cell_maps = std::vector<Simplex_cell_map>; + + /** \brief Type of a map from cells in the cell complex to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + */ + using Cell_simplex_map = std::map<Hasse_cell*, Simplex_handle>; + + /** \brief Type of a map from vertex cells in the cell complex to the permutahedral representations + * of their Cartesian coordinates. + */ + using Cell_point_map = std::map<Hasse_cell*, Eigen::VectorXd>; + + private: + Hasse_cell* insert_cell(const Simplex_handle& simplex, std::size_t cell_d, bool is_boundary) { + Simplex_cell_maps& simplex_cell_maps = (is_boundary ? boundary_simplex_cell_maps_ : interior_simplex_cell_maps_); +#ifdef DEBUG_TRACES + CC_detail_list& cc_detail_list = + (is_boundary ? cc_boundary_detail_lists[cell_d] : cc_interior_detail_lists[cell_d]); + cc_detail_list.emplace_back(simplex); +#endif + Simplex_cell_map& simplex_cell_map = simplex_cell_maps[cell_d]; + auto map_it = simplex_cell_map.find(simplex); + if (map_it == simplex_cell_map.end()) { + hasse_cells_.push_back(new Hasse_cell(is_boundary, cell_d)); + Hasse_cell* new_cell = hasse_cells_.back(); + simplex_cell_map.emplace(simplex, new_cell); + cell_simplex_map_.emplace(new_cell, simplex); +#ifdef DEBUG_TRACES + cc_detail_list.back().status_ = CC_detail_info::Result_type::inserted; +#endif + return new_cell; + } +#ifdef DEBUG_TRACES + CC_detail_info& cc_info = cc_detail_list.back(); + cc_info.trigger_ = to_string(map_it->first); + cc_info.status_ = CC_detail_info::Result_type::self; +#endif + return map_it->second; + } + + void expand_level(std::size_t cell_d) { + bool is_manifold_with_boundary = boundary_simplex_cell_maps_.size() > 0; + for (auto& sc_pair : interior_simplex_cell_maps_[cell_d - 1]) { + const Simplex_handle& simplex = sc_pair.first; + Hasse_cell* cell = sc_pair.second; + for (Simplex_handle coface : simplex.coface_range(cod_d_ + cell_d)) { + Hasse_cell* new_cell = insert_cell(coface, cell_d, false); + new_cell->get_boundary().emplace_back(cell, 1); + } + } + + if (is_manifold_with_boundary) { + for (auto& sc_pair : boundary_simplex_cell_maps_[cell_d - 1]) { + const Simplex_handle& simplex = sc_pair.first; + Hasse_cell* cell = sc_pair.second; + if (cell_d != intr_d_) + for (Simplex_handle coface : simplex.coface_range(cod_d_ + cell_d + 1)) { + Hasse_cell* new_cell = insert_cell(coface, cell_d, true); + new_cell->get_boundary().emplace_back(cell, 1); + } + auto map_it = interior_simplex_cell_maps_[cell_d].find(simplex); + if (map_it == interior_simplex_cell_maps_[cell_d].end()) + std::cerr << "Cell_complex::expand_level error: A boundary cell does not have an interior counterpart.\n"; + else { + Hasse_cell* i_cell = map_it->second; + i_cell->get_boundary().emplace_back(cell, 1); + } + } + } + } + + void construct_complex_(const Out_simplex_map_& out_simplex_map) { +#ifdef DEBUG_TRACES + cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size()); + cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size()); + cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size()); +#endif + for (auto& os_pair : out_simplex_map) { + const Simplex_handle& simplex = os_pair.first; + const Eigen::VectorXd& point = os_pair.second; + Hasse_cell* new_cell = insert_cell(simplex, 0, false); + cell_point_map_.emplace(new_cell, point); + } + for (std::size_t cell_d = 1; + cell_d < interior_simplex_cell_maps_.size() && !interior_simplex_cell_maps_[cell_d - 1].empty(); ++cell_d) { + expand_level(cell_d); + } + } + + void construct_complex_(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map) { +#ifdef DEBUG_TRACES + cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size()); + cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size()); + cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size()); + cc_boundary_summary_lists.resize(boundary_simplex_cell_maps_.size()); + cc_boundary_prejoin_lists.resize(boundary_simplex_cell_maps_.size()); + cc_boundary_detail_lists.resize(boundary_simplex_cell_maps_.size()); +#endif + for (auto& os_pair : boundary_simplex_map) { + const Simplex_handle& simplex = os_pair.first; + const Eigen::VectorXd& point = os_pair.second; + Hasse_cell* new_cell = insert_cell(simplex, 0, true); + cell_point_map_.emplace(new_cell, point); + } + for (auto& os_pair : interior_simplex_map) { + const Simplex_handle& simplex = os_pair.first; + const Eigen::VectorXd& point = os_pair.second; + Hasse_cell* new_cell = insert_cell(simplex, 0, false); + cell_point_map_.emplace(new_cell, point); + } +#ifdef DEBUG_TRACES + for (const auto& sc_pair : interior_simplex_cell_maps_[0]) + cc_interior_summary_lists[0].push_back(CC_summary_info(sc_pair)); + for (const auto& sc_pair : boundary_simplex_cell_maps_[0]) + cc_boundary_summary_lists[0].push_back(CC_summary_info(sc_pair)); +#endif + + for (std::size_t cell_d = 1; + cell_d < interior_simplex_cell_maps_.size() && !interior_simplex_cell_maps_[cell_d - 1].empty(); ++cell_d) { + expand_level(cell_d); + +#ifdef DEBUG_TRACES + for (const auto& sc_pair : interior_simplex_cell_maps_[cell_d]) + cc_interior_summary_lists[cell_d].push_back(CC_summary_info(sc_pair)); + if (cell_d < boundary_simplex_cell_maps_.size()) + for (const auto& sc_pair : boundary_simplex_cell_maps_[cell_d]) + cc_boundary_summary_lists[cell_d].push_back(CC_summary_info(sc_pair)); +#endif + } + } + + public: + /** + * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold + * without boundary embedded in the \f$ d \f$-dimensional Euclidean space + * from the output of the class Gudhi::Manifold_tracing. + * + * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$ + * in the ambient triangulation that intersect the relative interior of the manifold + * to the intersection points. + */ + void construct_complex(const Out_simplex_map_& out_simplex_map) { + interior_simplex_cell_maps_.resize(intr_d_ + 1); + if (!out_simplex_map.empty()) cod_d_ = out_simplex_map.begin()->first.dimension(); + construct_complex_(out_simplex_map); + } + + /** + * \brief Constructs the skeleton of the cell complex that approximates + * an \f$m\f$-dimensional manifold without boundary embedded + * in the \f$d\f$-dimensional Euclidean space + * up to a limit dimension from the output of the class Gudhi::Manifold_tracing. + * + * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$ + * in the ambient triangulation that intersect the relative interior of the manifold + * to the intersection points. + * \param[in] limit_dimension The dimension of the constructed skeleton. + */ + void construct_complex(const Out_simplex_map_& out_simplex_map, std::size_t limit_dimension) { + interior_simplex_cell_maps_.resize(limit_dimension + 1); + if (!out_simplex_map.empty()) cod_d_ = out_simplex_map.begin()->first.dimension(); + construct_complex_(out_simplex_map); + } + + /** + * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold + * with boundary embedded in the \f$ d \f$-dimensional Euclidean space + * from the output of the class Gudhi::Manifold_tracing. + * + * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$ + * in the ambient triangulation that intersect the relative interior of the manifold + * to the intersection points. + * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$ + * in the ambient triangulation that intersect the boundary of the manifold + * to the intersection points. + */ + void construct_complex(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map) { + interior_simplex_cell_maps_.resize(intr_d_ + 1); + boundary_simplex_cell_maps_.resize(intr_d_); + if (!interior_simplex_map.empty()) cod_d_ = interior_simplex_map.begin()->first.dimension(); + construct_complex_(interior_simplex_map, boundary_simplex_map); + } + + /** + * \brief Constructs the skeleton of the cell complex that approximates + * an \f$m\f$-dimensional manifold with boundary embedded + * in the \f$d\f$-dimensional Euclidean space + * up to a limit dimension from the output of the class Gudhi::Manifold_tracing. + * + * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$ + * in the ambient triangulation that intersect the relative interior of the manifold + * to the intersection points. + * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$ + * in the ambient triangulation that intersect the boundary of the manifold + * to the intersection points. + * \param[in] limit_dimension The dimension of the constructed skeleton. + */ + void construct_complex(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map, + std::size_t limit_dimension) { + interior_simplex_cell_maps_.resize(limit_dimension + 1); + boundary_simplex_cell_maps_.resize(limit_dimension); + if (!interior_simplex_map.empty()) cod_d_ = interior_simplex_map.begin()->first.dimension(); + construct_complex_(interior_simplex_map, boundary_simplex_map); + } + + /** + * \brief Returns the dimension of the cell complex. + */ + std::size_t intrinsic_dimension() const { return intr_d_; } + + /** + * \brief Returns a vector of maps from the cells of various dimensions in the interior + * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + */ + const Simplex_cell_maps& interior_simplex_cell_maps() const { return interior_simplex_cell_maps_; } + + /** + * \brief Returns a vector of maps from the cells of various dimensions on the boundary + * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + */ + const Simplex_cell_maps& boundary_simplex_cell_maps() const { return boundary_simplex_cell_maps_; } + + /** + * \brief Returns a map from the cells of a given dimension in the interior + * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + * + * \param[in] cell_d The dimension of the cells. + */ + const Simplex_cell_map& interior_simplex_cell_map(std::size_t cell_d) const { + return interior_simplex_cell_maps_[cell_d]; + } + + /** + * \brief Returns a map from the cells of a given dimension on the boundary + * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + * + * \param[in] cell_d The dimension of the cells. + */ + const Simplex_cell_map& boundary_simplex_cell_map(std::size_t cell_d) const { + return boundary_simplex_cell_maps_[cell_d]; + } + + /** + * \brief Returns a map from the cells in the cell complex of type Gudhi::Hasse_cell + * to the permutahedral representations of the corresponding simplices in the + * ambient triangulation. + */ + const Cell_simplex_map& cell_simplex_map() const { return cell_simplex_map_; } + + /** + * \brief Returns a map from the vertex cells in the cell complex of type Gudhi::Hasse_cell + * to their Cartesian coordinates. + */ + const Cell_point_map& cell_point_map() const { return cell_point_map_; } + + /** + * \brief Constructor for the class Cell_complex. + * + * \param[in] intrinsic_dimension The dimension of the cell complex. + */ + Cell_complex(std::size_t intrinsic_dimension) : intr_d_(intrinsic_dimension) {} + + ~Cell_complex() { + for (Hasse_cell* hs_ptr : hasse_cells_) delete hs_ptr; + } + + private: + std::size_t intr_d_, cod_d_; + Simplex_cell_maps interior_simplex_cell_maps_, boundary_simplex_cell_maps_; + Cell_simplex_map cell_simplex_map_; + Cell_point_map cell_point_map_; + std::vector<Hasse_cell*> hasse_cells_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h new file mode 100644 index 00000000..9b57da3c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h @@ -0,0 +1,285 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2017 Swansea University UK + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef HASSE_DIAGRAM_CELL_H +#define HASSE_DIAGRAM_CELL_H + +#include <vector> +#include <utility> // for std::pair +#include <ostream> +#include <string> +#include <type_traits> // for std::is_same +#include <cstdlib> // for std::size_t + +namespace Gudhi { +namespace Hasse_diagram { + +template <typename Cell_type> +class Hasse_diagram; + +/** + * \class Hasse_diagram_cell + * \brief Data structure to store a cell in a Hasse diagram. + * + * \ingroup Hasse_diagram + * + * \details + * The use and interfaces of this Hasse diagram cell is limited to the \ref coxeter_triangulation implementation. + * + * This is a data structure to store a cell in a general Hasse diagram data structure. It stores the following + * information about the cell: References to boundary and coBoundary elements, dimension of a cell and its filtration. + * It also allow to store any additional information of a type Additional_information which is a template parameter of + * the class (set by default to void). + * + * The complex is a template class requiring the following parameters: + * Incidence_type_ - determine the type of incidence coefficients. Use integers in most general case. + * Filtration_type_ - type of filtration of cells. + * Additional_information_ (set by default to void) - allows to store any + * additional information in the cells of Hasse diagrams. + * + */ +template <typename Incidence_type_, typename Filtration_type_, typename Additional_information_ = void> +class Hasse_diagram_cell { + public: + typedef Incidence_type_ Incidence_type; + typedef Filtration_type_ Filtration_type; + typedef Additional_information_ Additional_information; + using Cell_range = std::vector<std::pair<Hasse_diagram_cell*, Incidence_type> >; + + /** + * Default constructor. + **/ + Hasse_diagram_cell() : dimension(0), position(0), deleted_(false) {} + + /** + * Constructor of a cell of dimension dim. + **/ + Hasse_diagram_cell(int dim) : dimension(dim), position(0), deleted_(false) {} + + /** + * Constructor of a cell of dimension dim. + **/ + Hasse_diagram_cell(int dim, Filtration_type filt_) + : dimension(dim), position(0), deleted_(false), filtration(filt_) {} + + /** + * Constructor of a cell of dimension dim with a given boundary. + **/ + Hasse_diagram_cell(const Cell_range& boundary_, int dim) + : dimension(dim), boundary(boundary_), position(0), deleted_(false) {} + + /** + * Constructor of a cell of dimension dim with a given boundary and coboundary. + **/ + Hasse_diagram_cell(const Cell_range& boundary_, const Cell_range& coboundary_, int dim) + : dimension(dim), boundary(boundary_), coBoundary(coboundary_), position(0), deleted_(false) {} + + /** + * Constructor of a cell of dimension dim with a given boundary, coboundary and + * additional information. + **/ + Hasse_diagram_cell(const Cell_range& boundary_, const Cell_range& coboundary_, const Additional_information& ai, + int dim) + : dimension(dim), + boundary(boundary_), + coBoundary(coboundary_), + additional_info(ai), + position(0), + deleted_(false) {} + + /** + * Constructor of a cell of dimension dim having given additional information. + **/ + Hasse_diagram_cell(Additional_information ai, int dim) + : dimension(dim), additional_info(ai), position(0), deleted_(false) {} + + /** + * Procedure to get the boundary of a fiven cell. The output format + * is a vector of pairs of pointers to boundary elements and incidence + * coefficients. + **/ + inline Cell_range& get_boundary() { return this->boundary; } + + /** + * Procedure to get the coboundary of a fiven cell. The output format + * is a vector of pairs of pointers to coboundary elements and incidence + * coefficients. + **/ + inline Cell_range& get_coBoundary() { return this->coBoundary; } + + /** + * Procedure to get the dimension of a cell. + **/ + inline int& get_dimension() { return this->dimension; } + + /** + * Procedure to get additional information about the cell.s + **/ + inline Additional_information& get_additional_information() { return this->additional_info; } + + /** + * Procedure to retrieve the position of the cell in the structure. It is used in + * the implementation of Hasse diagram and set by it. Note that removal of + * cell and subsequent call of clean_up_the_structure will change those + * positions. + **/ + inline unsigned& get_position() { return this->position; } + + /** + * Accessing the filtration of the cell. + **/ + inline Filtration_type& get_filtration() { + // std::cout << "Accessing the filtration of a cell : " << *this << std::endl; + return this->filtration; + } + + /** + * A procedure used to check if the cell is deleted. It is used by the + * subsequent implementation of Hasse diagram that is absed on lazy + * delete. + **/ + inline bool deleted() { return this->deleted_; } + + template <typename Cell_type> + friend class Hasse_diagram; + + template <typename Cell_type> + friend class is_before_in_filtration; + + template <typename Complex_type, typename Cell_type> + friend std::vector<Cell_type*> convert_to_vector_of_Cell_type(Complex_type& cmplx); + + /** + * Procedure to remove deleted boundary and coboundary elements from the + * vectors of boundary and coboundary elements of this cell. + **/ + void remove_deleted_elements_from_boundary_and_coboundary() { + Cell_range new_boundary; + new_boundary.reserve(this->boundary.size()); + for (std::size_t bd = 0; bd != this->boundary.size(); ++bd) { + if (!this->boundary[bd].first->deleted()) { + new_boundary.push_back(this->boundary[bd]); + } + } + this->boundary.swap(new_boundary); + + Cell_range new_coBoundary; + new_coBoundary.reserve(this->coBoundary.size()); + for (std::size_t cbd = 0; cbd != this->coBoundary.size(); ++cbd) { + if (!this->coBoundary[cbd].first->deleted()) { + new_coBoundary.push_back(this->coBoundary[cbd]); + } + } + this->coBoundary.swap(new_coBoundary); + } + + /** + * Writing to a stream operator. + **/ + friend std::ostream& operator<<( + std::ostream& out, const Hasse_diagram_cell<Incidence_type, Filtration_type, Additional_information>& c) { + // cout << "position : " << c.position << ", dimension : " << c.dimension << ", filtration: " << c.filtration << ", + // size of boundary : " << c.boundary.size() << "\n"; + out << c.position << " " << c.dimension << " " << c.filtration << std::endl; + for (std::size_t bd = 0; bd != c.boundary.size(); ++bd) { + // do not write out the cells that has been deleted + if (c.boundary[bd].first->deleted()) continue; + out << c.boundary[bd].first->position << " " << c.boundary[bd].second << " "; + } + out << std::endl; + return out; + } + + /** + * Procedure that return vector of pointers to boundary elements of a given cell. + **/ + inline std::vector<Hasse_diagram_cell*> get_list_of_boundary_elements() { + std::vector<Hasse_diagram_cell*> result; + std::size_t size_of_boundary = this->boundary.size(); + result.reserve(size_of_boundary); + for (std::size_t bd = 0; bd != size_of_boundary; ++bd) { + result.push_back(this->boundary[bd].first); + } + return result; + } + + /** + * Procedure that return vector of positios of boundary elements of a given cell. + **/ + inline std::vector<unsigned> get_list_of_positions_of_boundary_elements() { + std::vector<unsigned> result; + std::size_t size_of_boundary = this->boundary.size(); + result.reserve(size_of_boundary); + for (std::size_t bd = 0; bd != size_of_boundary; ++bd) { + result.push_back(this->boundary[bd].first->position); + } + return result; + } + + /** + * Function that display a string being a signature of a structure. + * Used mainly for debugging purposes. + **/ + std::string full_signature_of_the_structure() { + std::string result; + result += "dimension: "; + result += std::to_string(this->dimension); + result += " filtration: "; + result += std::to_string(this->filtration); + result += " position: "; + result += std::to_string(this->position); + result += " deleted_: "; + result += std::to_string(this->deleted_); + + // if the Additional_information is not void, add them to + // the signature as well. + if (std::is_same<Additional_information, void>::value) { + result += " Additional_information: "; + result += std::to_string(this->additional_info); + } + result += " boundary "; + for (std::size_t bd = 0; bd != this->boundary.size(); ++bd) { + result += "( " + std::to_string(this->boundary[bd].first->position); + result += " " + std::to_string(this->boundary[bd].second); + result += ") "; + } + + result += " coBoundary "; + for (std::size_t cbd = 0; cbd != this->coBoundary.size(); ++cbd) { + result += "( " + std::to_string(this->coBoundary[cbd].first->position); + result += " " + std::to_string(this->coBoundary[cbd].second); + result += ") "; + } + + return result; + } + + protected: + Cell_range boundary; + Cell_range coBoundary; + int dimension; + Additional_information additional_info; + unsigned position; + bool deleted_; + Filtration_type filtration; + + /** + * A procedure to delete a cell. It is a private function of the Hasse_diagram_cell + * class, since in the Hasse_diagram class I want to have a control + * of removal of cells. Therefore, to remove cell please use + * remove_cell in the Hasse_diagram structure. + **/ + void delete_cell() { this->deleted_ = true; } +}; // Hasse_diagram_cell + +} // namespace Hasse_diagram +} // namespace Gudhi + +#endif // CELL_H diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h new file mode 100644 index 00000000..5543c2fb --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h @@ -0,0 +1,40 @@ +/* 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 QUERY_RESULT_H_ +#define QUERY_RESULT_H_ + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Query_result + * \brief The result of a query by an oracle such as Implicit_manifold_intersection_oracle. + * + * \tparam Simplex_handle The class of the query simplex. + * + * \ingroup coxeter_triangulation + */ +template <class Simplex_handle> +struct Query_result { + /** \brief The potentially lower-dimensional face of the query simplex + * that contains the intersection point. OBSOLETE: as the snapping is removed. */ + // Simplex_handle face; + /** \brief The intersection point. */ + Eigen::VectorXd intersection; + /** \brief True if the query simplex intersects the manifold. */ + bool success; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h new file mode 100644 index 00000000..873c5c9b --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h @@ -0,0 +1,219 @@ +/* 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 { + +/** + * \class Freudenthal_triangulation + * \brief A class that stores any affine transformation of the Freudenthal-Kuhn + * triangulation. + * + * \ingroup coxeter_triangulation + * + * \details The data structure is a record that consists of a matrix + * that represents the linear transformation of the Freudenthal-Kuhn triangulation + * and a vector that represents the offset. + * + * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation. + * Needs to be a model of SimplexInCoxeterTriangulation. + */ +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: + /** \brief Type of the simplices in the triangulation. */ + using Simplex_handle = Permutahedral_representation_; + + /** \brief Type of the vertices in the triangulation. */ + using Vertex_handle = typename Permutahedral_representation_::Vertex; + + /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension. + * @param[in] dimension The dimension of the triangulation. + */ + Freudenthal_triangulation(std::size_t dimension) + : Freudenthal_triangulation(dimension, Matrix::Identity(dimension, dimension), Vector::Zero(dimension)) { + is_freudenthal_ = true; + } + + /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under + * a linear transformation by a given matrix. + * @param[in] dimension The dimension of the triangulation. + * @param[in] matrix The matrix that defines the linear transformation. + * Needs to be invertible. + */ + Freudenthal_triangulation(std::size_t dimension, const Matrix& matrix) + : Freudenthal_triangulation(dimension, matrix, Vector::Zero(dimension)) {} + + /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under + * an affine transformation by a given matrix and a translation vector. + * @param[in] dimension The dimension of the triangulation. + * @param[in] matrix The matrix that defines the linear transformation. + * Needs to be invertible. + * @param[in] offset The offset vector. + * + * @exception std::invalid_argument In debug mode, if offset size is different from 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'")); + } + + /** \brief Dimension of the triangulation. */ + unsigned dimension() const { return dimension_; } + + /** \brief Matrix that defines the linear transformation of the triangulation. */ + const Matrix& matrix() const { return matrix_; } + + /** \brief Vector that defines the offset of the triangulation. */ + const Vector& offset() const { return offset_; } + + /** \brief Change the linear transformation matrix to a given value. + * @param[in] matrix New value of the linear transformation matrix. + */ + void change_matrix(const Eigen::MatrixXd& matrix) { + matrix_ = matrix; + colpivhouseholderqr_ = matrix.colPivHouseholderQr(); + is_freudenthal_ = false; + } + + /** \brief Change the offset vector to a given value. + * @param[in] offset New value of the offset vector. + */ + void change_offset(const Eigen::VectorXd& offset) { + offset_ = offset; + is_freudenthal_ = false; + } + + /** \brief Returns the permutahedral representation of the simplex in the + * triangulation that contains a given query point. + * \details Using the additional parameter scale, the search can be done in a + * triangulation that shares the origin, but is scaled by a given factor. + * This parameter can be useful to simulate the point location in a subdivided + * triangulation. + * The returned simplex is always minimal by inclusion. + * + * \tparam Point_d A class that represents a point in d-dimensional Euclidean space. + * The coordinates should be random-accessible. Needs to provide the method size(). + * + * @param[in] point The query point. + * @param[in] scale The scale of the triangulation. + * + * @exception std::invalid_argument In debug mode, if point dimension is different from triangulation one. + */ + 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; + } + + /** \brief Returns the Cartesian coordinates of the given vertex. + * \details Using the additional parameter scale, the search can be done in a + * triangulation that shares the origin, but is scaled by a given factor. + * This parameter can be useful to simulate the computation of Cartesian coordinates + * of a vertex in a subdivided triangulation. + * @param[in] vertex The query vertex. + * @param[in] scale The scale of the triangulation. + */ + 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_; + } + + /** \brief Returns the Cartesian coordinates of the barycenter of a given simplex. + * \details Using the additional parameter scale, the search can be done in a + * triangulation that shares the origin, but is scaled by a given factor. + * This parameter can be useful to simulate the computation of Cartesian coordinates + * of the barycenter of a simplex in a subdivided triangulation. + * @param[in] simplex The query simplex. + * @param[in] scale The scale of the triangulation. + */ + 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 diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h new file mode 100644 index 00000000..0533bb83 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h @@ -0,0 +1,157 @@ +/* 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 FUNCTIONS_CARTESIAN_PRODUCT_H_ +#define FUNCTIONS_CARTESIAN_PRODUCT_H_ + +#include <cstdlib> +#include <tuple> +#include <type_traits> // for std::enable_if +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/* Get the domain dimension of the tuple of functions. + */ +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I == sizeof...(T), std::size_t>::type get_amb_d(const std::tuple<T...>& tuple) { + return 0; +} + +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I != sizeof...(T), std::size_t>::type get_amb_d(const std::tuple<T...>& tuple) { + return std::get<I>(tuple).amb_d() + get_amb_d<I + 1, T...>(tuple); +} + +/* Get the codomain dimension of the tuple of functions. + */ +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I == sizeof...(T), std::size_t>::type get_cod_d(const std::tuple<T...>& tuple) { + return 0; +} + +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I != sizeof...(T), std::size_t>::type get_cod_d(const std::tuple<T...>& tuple) { + return std::get<I>(tuple).cod_d() + get_cod_d<I + 1, T...>(tuple); +} + +/* Get the seed of the tuple of functions. + */ +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I == sizeof...(T), void>::type get_seed(const std::tuple<T...>& tuple, + Eigen::VectorXd& point, std::size_t i = 0) {} + +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I != sizeof...(T), void>::type get_seed(const std::tuple<T...>& tuple, + Eigen::VectorXd& point, std::size_t i = 0) { + const auto& f = std::get<I>(tuple); + std::size_t n = f.amb_d(); + Eigen::VectorXd seed = f.seed(); + for (std::size_t j = 0; j < n; ++j) point(i + j) = seed(j); + get_seed<I + 1, T...>(tuple, point, i + n); +} + +/* Get the seed of the tuple of functions. + */ +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I == sizeof...(T), void>::type get_value(const std::tuple<T...>& tuple, + const Eigen::VectorXd& x, + Eigen::VectorXd& point, std::size_t i = 0, + std::size_t j = 0) {} + +template <std::size_t I = 0, typename... T> +inline typename std::enable_if<I != sizeof...(T), void>::type get_value(const std::tuple<T...>& tuple, + const Eigen::VectorXd& x, + Eigen::VectorXd& point, std::size_t i = 0, + std::size_t j = 0) { + const auto& f = std::get<I>(tuple); + std::size_t n = f.amb_d(); + std::size_t k = f.cod_d(); + Eigen::VectorXd x_i(n); + for (std::size_t l = 0; l < n; ++l) x_i(l) = x(i + l); + Eigen::VectorXd res = f(x_i); + for (std::size_t l = 0; l < k; ++l) point(j + l) = res(l); + get_value<I + 1, T...>(tuple, x, point, i + n, j + k); +} + +/** + * \class Cartesian_product + * \brief Constructs the function the zero-set of which is the Cartesian product + * of the zero-sets of some given functions. + * + * \tparam Functions A pack template parameter for functions. All functions should be models of + * the concept FunctionForImplicitManifold. + */ +template <class... Functions> +struct Cartesian_product { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result(cod_d_); + get_value(function_tuple_, p, result, 0, 0); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return amb_d_; } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return cod_d_; } + + /** \brief Returns a point on the zero-set. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result(amb_d_); + get_seed(function_tuple_, result, 0); + return result; + } + + /** + * \brief Constructor of the Cartesian product function. + * + * @param[in] functions The functions the zero-sets of which are factors in the + * Cartesian product of the resulting function. + */ + Cartesian_product(const Functions&... functions) : function_tuple_(std::make_tuple(functions...)) { + amb_d_ = get_amb_d(function_tuple_); + cod_d_ = get_cod_d(function_tuple_); + } + + private: + std::tuple<Functions...> function_tuple_; + std::size_t amb_d_, cod_d_; +}; + +/** + * \brief Static constructor of a Cartesian product function. + * + * @param[in] functions The functions the zero-sets of which are factors in the + * Cartesian product of the resulting function. + * + * \tparam Functions A pack template parameter for functions. All functions should be models of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <typename... Functions> +Cartesian_product<Functions...> make_product_function(const Functions&... functions) { + return Cartesian_product<Functions...>(functions...); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h new file mode 100644 index 00000000..0603afd8 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h @@ -0,0 +1,64 @@ +/* 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 FUNCTIONS_CONSTANT_FUNCTION_H_ +#define FUNCTIONS_CONSTANT_FUNCTION_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Constant_function + * \brief A class that encodes a constant function from R^d to R^k. + * This class does not have any implicit manifold in correspondence. + */ +struct Constant_function { + /** \brief Value of the function at a specified point. The value is constant. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + return value_; + } + + /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */ + std::size_t amb_d() const { return d_; }; + + /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */ + std::size_t cod_d() const { return k_; }; + + /** \brief No seed point is available. Throws an exception on evocation. */ + Eigen::VectorXd seed() const { throw "Seed invoked on a constant function.\n"; } + + Constant_function() {} + + /** + * \brief Constructor of a constant function from R^d to R^m. + * + * @param[in] d The domain dimension. + * @param[in] k The codomain dimension. + * @param[in] value The constant value of the function. + */ + Constant_function(std::size_t d, std::size_t k, const Eigen::VectorXd& value) : d_(d), k_(k), value_(value) {} + + private: + std::size_t d_, k_; + Eigen::VectorXd value_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h new file mode 100644 index 00000000..e1fe868f --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h @@ -0,0 +1,93 @@ +/* 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 FUNCTIONS_EMBED_IN_RD_H_ +#define FUNCTIONS_EMBED_IN_RD_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Embed_in_Rd + * \brief Embedding of an implicit manifold in a higher dimension. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + */ +template <class Function_> +struct Embed_in_Rd { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd x = p; + Eigen::VectorXd x_k(fun_.amb_d()), x_rest(d_ - fun_.amb_d()); + for (std::size_t i = 0; i < fun_.amb_d(); ++i) x_k(i) = x(i); + for (std::size_t i = fun_.amb_d(); i < d_; ++i) x_rest(i - fun_.amb_d()) = x(i); + Eigen::VectorXd result = fun_(x_k); + result.conservativeResize(this->cod_d()); + for (std::size_t i = fun_.cod_d(); i < this->cod_d(); ++i) result(i) = x_rest(i - fun_.cod_d()); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return d_; } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return d_ - (fun_.amb_d() - fun_.cod_d()); } + + /** \brief Returns a point on the zero-set of the embedded function. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = fun_.seed(); + result.conservativeResize(d_); + for (std::size_t l = fun_.amb_d(); l < d_; ++l) result(l) = 0; + return result; + } + + /** + * \brief Constructor of the embedding function. + * + * @param[in] function The function to be embedded in higher dimension. + * @param[in] d Embedding dimension. + */ + Embed_in_Rd(const Function_& function, std::size_t d) : fun_(function), d_(d) {} + + private: + Function_ fun_; + std::size_t d_; +}; + +/** + * \brief Static constructor of an embedding function. + * + * @param[in] function The function to be embedded in higher dimension. + * @param[in] d Embedding dimension. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Embed_in_Rd<Function_> make_embedding(const Function_& function, std::size_t d) { + return Embed_in_Rd<Function_>(function, d); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h new file mode 100644 index 00000000..8911f990 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h @@ -0,0 +1,110 @@ +/* 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 FUNCTIONS_FUNCTION_SM_IN_RD_H_ +#define FUNCTIONS_FUNCTION_SM_IN_RD_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_Sm_in_Rd + * \brief A class for the function that defines an m-dimensional implicit sphere embedded + * in the d-dimensional Euclidean space. + */ +struct Function_Sm_in_Rd { + /** \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd x = p; + for (std::size_t i = 0; i < d_; ++i) x(i) -= center_[i]; + Eigen::VectorXd result = Eigen::VectorXd::Zero(k_); + for (std::size_t i = 0; i < m_ + 1; ++i) result(0) += x(i) * x(i); + result(0) -= r_ * r_; + for (std::size_t j = 1; j < k_; ++j) result(j) = x(m_ + j); + return result; + } + + /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */ + std::size_t amb_d() const { return d_; }; + + /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */ + std::size_t cod_d() const { return k_; }; + + /** \brief Returns a point on the sphere. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = Eigen::VectorXd::Zero(d_); + result(0) += r_; + for (std::size_t i = 0; i < d_; ++i) result(i) += center_[i]; + return result; + } + + /** + * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded + * in the d-dimensional Euclidean space. + * + * @param[in] r The radius of the sphere. + * @param[in] m The dimension of the sphere. + * @param[in] d The ambient dimension of the sphere. + * @param[in] center The center of the sphere. + */ + Function_Sm_in_Rd(double r, std::size_t m, std::size_t d, Eigen::VectorXd center) + : m_(m), k_(d - m), d_(d), r_(r), center_(center) {} + + /** + * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded + * in the d-dimensional Euclidean space centered at the origin. + * + * @param[in] r The radius of the sphere. + * @param[in] m The dimension of the sphere. + * @param[in] d The ambient dimension of the sphere. + */ + Function_Sm_in_Rd(double r, std::size_t m, std::size_t d) + : m_(m), k_(d - m), d_(d), r_(r), center_(Eigen::VectorXd::Zero(d_)) {} + + /** + * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded + * in the (m+1)-dimensional Euclidean space. + * + * @param[in] r The radius of the sphere. + * @param[in] m The dimension of the sphere. + * @param[in] center The center of the sphere. + */ + Function_Sm_in_Rd(double r, std::size_t m, Eigen::VectorXd center) + : m_(m), k_(1), d_(m_ + 1), r_(r), center_(center) {} + + /** + * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded + * in the (m+1)-dimensional Euclidean space centered at the origin. + * + * @param[in] r The radius of the sphere. + * @param[in] m The dimension of the sphere. + */ + Function_Sm_in_Rd(double r, std::size_t m) : m_(m), k_(1), d_(m_ + 1), r_(r), center_(Eigen::VectorXd::Zero(d_)) {} + + Function_Sm_in_Rd(const Function_Sm_in_Rd& rhs) : Function_Sm_in_Rd(rhs.r_, rhs.m_, rhs.d_, rhs.center_) {} + + private: + std::size_t m_, k_, d_; + double r_; + Eigen::VectorXd center_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h new file mode 100644 index 00000000..a9e2d507 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h @@ -0,0 +1,90 @@ +/* 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 FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_ +#define FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_affine_plane_in_Rd + * \brief A class for the function that defines an m-dimensional implicit affine plane + * embedded in d-dimensional Euclidean space. + */ +struct Function_affine_plane_in_Rd { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result = normal_matrix_.transpose() * (p - off_); + return result; + } + + /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */ + std::size_t amb_d() const { return d_; }; + + /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */ + std::size_t cod_d() const { return k_; }; + + /** \brief Returns a point on the affine plane. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = off_; + return result; + } + + /** + * \brief Constructor of the function that defines an m-dimensional implicit affine + * plane in the d-dimensional Euclidean space. + * + * @param[in] normal_matrix A normal matrix of the affine plane. The number of rows should + * correspond to the ambient dimension, the number of columns should correspond to + * the size of the normal basis (codimension). + * @param[in] offset The offset vector of the affine plane. + * The dimension of the vector should be the ambient dimension of the manifold. + */ + Function_affine_plane_in_Rd(const Eigen::MatrixXd& normal_matrix, const Eigen::VectorXd& offset) + : normal_matrix_(normal_matrix), d_(normal_matrix.rows()), k_(normal_matrix.cols()), off_(offset) { + normal_matrix_.colwise().normalize(); + } + + /** + * \brief Constructor of the function that defines an m-dimensional implicit affine + * plane in the d-dimensional Euclidean space that passes through origin. + * + * @param[in] normal_matrix A normal matrix of the affine plane. The number of rows should + * correspond to the ambient dimension, the number of columns should correspond to + * the size of the normal basis (codimension). + */ + Function_affine_plane_in_Rd(const Eigen::MatrixXd& normal_matrix) + : normal_matrix_(normal_matrix), + d_(normal_matrix.rows()), + k_(normal_matrix.cols()), + off_(Eigen::VectorXd::Zero(d_)) { + normal_matrix_.colwise().normalize(); + } + + private: + Eigen::MatrixXd normal_matrix_; + std::size_t d_, k_; + Eigen::VectorXd off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h new file mode 100644 index 00000000..620446da --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h @@ -0,0 +1,80 @@ +/* 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 FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_ +#define FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::pow + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_chair_in_R3 + * \brief A class that encodes the function, the zero-set of which is a so-called + * "chair" surface embedded in R^3. + */ +struct Function_chair_in_R3 { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2]; + Eigen::VectorXd result(cod_d()); + result(0) = std::pow(x * x + y * y + z * z - a_ * k_ * k_, 2) - + b_ * ((z - k_) * (z - k_) - 2 * x * x) * ((z + k_) * (z + k_) - 2 * y * y); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return 3; } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return 1; } + + /** \brief Returns a point on the surface. */ + Eigen::VectorXd seed() const { + double t1 = a_ - b_; + double discr = t1 * t1 - (1.0 - b_) * (a_ * a_ - b_); + double z0 = k_ * std::sqrt((t1 + std::sqrt(discr)) / (1 - b_)); + Eigen::Vector3d result(off_[0], off_[1], z0 + off_[2]); + return result; + } + + /** + * \brief Constructor of the function that defines the 'chair' surface + * embedded in R^3. + * + * @param[in] a A numerical parameter. + * @param[in] b A numerical parameter. + * @param[in] k A numerical parameter. + * @param[in] off Offset vector. + */ + Function_chair_in_R3(double a = 0.8, double b = 0.4, double k = 1.0, Eigen::Vector3d off = Eigen::Vector3d::Zero()) + : a_(a), b_(b), k_(k), off_(off) {} + + protected: + double a_, b_, k_; + Eigen::Vector3d off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif + +// (x^2 + y^2 + z^2 - a*k^2)^2 - b*((z-k)^2 - 2*x^2)*((z+k)^2 - 2*y^2) +// sqrt(k/(1-b))*sqrt(a-b + sqrt((a-b)^2 - (1-b)*(a^2 - b)*k^2)) diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h new file mode 100644 index 00000000..f73c4280 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h @@ -0,0 +1,69 @@ +/* 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 FUNCTIONS_FUNCTION_IRON_IN_R3_H_ +#define FUNCTIONS_FUNCTION_IRON_IN_R3_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::pow + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_iron_in_R3 + * \brief A class that encodes the function, the zero-set of which is a surface + * embedded in R^3 that ressembles an iron. + */ +struct Function_iron_in_R3 { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + double x = p(0), y = p(1), z = p(2); + Eigen::VectorXd result(cod_d()); + result(0) = -std::pow(x, 6) / 300. - std::pow(y, 6) / 300. - std::pow(z, 6) / 300. + x * y * y * z / 2.1 + y * y + + std::pow(z - 2, 4) - 1; + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return 3; }; + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return 1; }; + + /** \brief Returns a point on the surface. */ + Eigen::VectorXd seed() const { + Eigen::Vector3d result(std::pow(4500, 1. / 6), 0, 0); + return result; + } + + /** + * \brief Constructor of the function that defines a surface embedded in R^3 + * that ressembles an iron. + * + * @param[in] off Offset vector. + */ + Function_iron_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero()) : off_(off) {} + + private: + Eigen::Vector3d off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h new file mode 100644 index 00000000..beb41e00 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h @@ -0,0 +1,85 @@ +/* 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 FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_ +#define FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::sqrt + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_lemniscate_revolution_in_R3 + * \brief A class that encodes the function, the zero-set of which is a surface of revolution + * around the x axis based on the lemniscate of Bernoulli embedded in R^3. + */ +struct Function_lemniscate_revolution_in_R3 { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2]; + Eigen::VectorXd result(cod_d()); + double x2 = x * x, y2 = y * y, z2 = z * z, a2 = a_ * a_; + double t1 = x2 + y2 + z2; + result(0) = t1 * t1 - 2 * a2 * (x2 - y2 - z2); + return result; + } + + /** \brief Returns the (ambient) domain dimension.*/ + std::size_t amb_d() const { return 3; }; + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return 1; }; + + /** \brief Returns a point on the surface. This seed point is only one of + * two necessary seed points for the manifold tracing algorithm. + * See the method seed2() for the other point. + */ + Eigen::VectorXd seed() const { + Eigen::Vector3d result(std::sqrt(2 * a_) + off_[0], off_[1], off_[2]); + return result; + } + + /** \brief Returns a point on the surface. This seed point is only one of + * two necessary seed points for the manifold tracing algorithm. + * See the method seed() for the other point. + */ + Eigen::VectorXd seed2() const { + Eigen::Vector3d result(-std::sqrt(2 * a_) + off_[0], off_[1], off_[2]); + return result; + } + + /** + * \brief Constructor of the function that defines a surface of revolution + * around the x axis based on the lemniscate of Bernoulli embedded in R^3. + * + * @param[in] a A numerical parameter. + * @param[in] off Offset vector. + */ + Function_lemniscate_revolution_in_R3(double a = 1, Eigen::Vector3d off = Eigen::Vector3d::Zero()) + : a_(a), off_(off) {} + + private: + double a_; + Eigen::Vector3d off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h new file mode 100644 index 00000000..f315d794 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h @@ -0,0 +1,84 @@ +/* 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 FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_ +#define FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_moment_curve_in_Rd + * \brief A class for the function that defines an implicit moment curve + * in the d-dimensional Euclidean space. + */ +struct Function_moment_curve_in_Rd { + /** \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result(k_); + for (std::size_t i = 1; i < d_; ++i) result(i - 1) = p(i) - p(0) * p(i - 1); + return result; + } + + /** \brief Returns the domain (ambient) dimension.. */ + std::size_t amb_d() const { return d_; }; + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return k_; }; + + /** \brief Returns a point on the moment curve. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = Eigen::VectorXd::Zero(d_); + return result; + } + + /** @brief Returns the radius of the moment curve. */ + double get_radius() const{ + return r_; + } + + /** + * \brief Constructor of the function that defines an implicit moment curve + * in the d-dimensional Euclidean space. + * + * @param[in] r Numerical parameter. + * @param[in] d The ambient dimension. + */ + Function_moment_curve_in_Rd(double r, std::size_t d) : k_(d - 1), d_(d), r_(r) {} + + /** + * \brief Constructor of the function that defines an implicit moment curve + * in the d-dimensional Euclidean space. + * + * @param[in] r Numerical parameter. + * @param[in] d The ambient dimension. + * @param[in] offset The offset of the moment curve. + */ + Function_moment_curve_in_Rd(double r, std::size_t d, Eigen::VectorXd& offset) + : k_(d - 1), d_(d), r_(r), off_(offset) {} + + private: + std::size_t k_, d_; + double r_; + Eigen::VectorXd off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h new file mode 100644 index 00000000..b54d3c74 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h @@ -0,0 +1,71 @@ +/* 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 FUNCTIONS_FUNCTION_TORUS_IN_R3_H_ +#define FUNCTIONS_FUNCTION_TORUS_IN_R3_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::sqrt + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_torus_in_R3 + * \brief A class that encodes the function, the zero-set of which is a torus + * surface embedded in R^3. + */ +struct Function_torus_in_R3 { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2]; + Eigen::VectorXd result(cod_d()); + result(0) = (z * z + (std::sqrt(x * x + y * y) - r_) * (std::sqrt(x * x + y * y) - r_) - R_ * R_); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return 3; }; + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return 1; }; + + /** \brief Returns a point on the surface. */ + Eigen::VectorXd seed() const { + Eigen::Vector3d result(R_ + r_ + off_[0], off_[1], off_[2]); + return result; + } + + /** + * \brief Constructor of the function that defines a torus embedded in R^3. + * + * @param[in] R The outer radius of the torus. + * @param[in] r The inner radius of the torus. + * @param[in] off Offset vector. + */ + Function_torus_in_R3(double R = 1, double r = 0.5, Eigen::Vector3d off = Eigen::Vector3d::Zero()) + : R_(R), r_(r), off_(off) {} + + private: + double R_, r_; + Eigen::Vector3d off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h new file mode 100644 index 00000000..df1f1eec --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h @@ -0,0 +1,78 @@ +/* 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 FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_ +#define FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function_whitney_umbrella_in_R3 + * \brief A class that encodes the function, the zero-set of which is the Whitney umbrella + * surface embedded in R^3. + */ +struct Function_whitney_umbrella_in_R3 { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2]; + Eigen::VectorXd result(cod_d()); + result(0) = x * x - y * y * z; + return result; + } + + /** \brief Returns the (ambient) domain dimension.*/ + std::size_t amb_d() const { return 3; }; + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return 1; }; + + /** \brief Returns a point on the surface. This seed point is only one of + * two necessary seed points for the manifold tracing algorithm. + * See the method seed2() for the other point. + */ + Eigen::VectorXd seed() const { + Eigen::Vector3d result(1 + off_[0], 1 + off_[1], 1 + off_[2]); + return result; + } + + /** \brief Returns a point on the surface. This seed point is only one of + * two necessary seed points for the manifold tracing algorithm. + * See the method seed() for the other point. + */ + Eigen::VectorXd seed2() const { + Eigen::Vector3d result(-1 + off_[0], -1 + off_[1], 1 + off_[2]); + return result; + } + + /** + * \brief Constructor of the function that defines the Whitney umbrella in R^3. + * + * @param[in] off Offset vector. + */ + Function_whitney_umbrella_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero()) : off_(off) {} + + private: + Eigen::Vector3d off_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h new file mode 100644 index 00000000..82e25bb9 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h @@ -0,0 +1,88 @@ +/* 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 FUNCTIONS_LINEAR_TRANSFORMATION_H_ +#define FUNCTIONS_LINEAR_TRANSFORMATION_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Linear_transformation + * \brief Transforms the zero-set of the function by a given linear transformation. + * The underlying function corresponds to f(M*x), where M is the transformation matrix. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + */ +template <class Function_> +struct Linear_transformation { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result = fun_(matrix_.householderQr().solve(p)); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief Returns a point on the zero-set. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = fun_.seed(); + result = matrix_ * result; + return result; + } + + /** + * \brief Constructor of a linearly transformed function. + * + * @param[in] function The function to be linearly transformed. + * @param[in] matrix The transformation matrix. Its dimension should be d*d, + * where d is the domain (ambient) dimension of 'function'. + */ + Linear_transformation(const Function_& function, const Eigen::MatrixXd& matrix) : fun_(function), matrix_(matrix) {} + + private: + Function_ fun_; + Eigen::MatrixXd matrix_; +}; + +/** + * \brief Static constructor of a linearly transformed function. + * + * @param[in] function The function to be linearly transformed. + * @param[in] matrix The transformation matrix. Its dimension should be d*d, + * where d is the domain (ambient) dimension of 'function'. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Linear_transformation<Function_> make_linear_transformation(const Function_& function, const Eigen::MatrixXd& matrix) { + return Linear_transformation<Function_>(function, matrix); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h new file mode 100644 index 00000000..fdf07f27 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h @@ -0,0 +1,84 @@ +/* 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 FUNCTIONS_NEGATION_H_ +#define FUNCTIONS_NEGATION_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + *\class Negation + * \brief Constructs the "minus" function. The zero-set is the same, but + * the values at other points are the negative of their original value. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + */ +template <class Function_> +struct Negation { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result = -fun_(p); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief Returns a point on the zero-set. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = fun_.seed(); + return result; + } + + /** + * \brief Constructor of the negative function. + * + * @param[in] function The function to be negated. + */ + Negation(const Function_& function) : fun_(function) {} + + private: + Function_ fun_; +}; + +/** + * \brief Static constructor of the negative function. + * + * @param[in] function The function to be translated. + * domain (ambient) dimension of 'function'. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Negation<Function_> negation(const Function_& function) { + return Negation<Function_>(function); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h new file mode 100644 index 00000000..22071d6d --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h @@ -0,0 +1,111 @@ +/* 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 FUNCTIONS_PL_APPROXIMATION_H_ +#define FUNCTIONS_PL_APPROXIMATION_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class PL_approximation + * \brief Constructs a piecewise-linear approximation of a function induced by + * an ambient triangulation. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * \tparam Triangulation The triangulation template parameter. Should be a model of + * the concept TriangulationForManifoldTracing. + */ +template <class Function_, class Triangulation_> +struct PL_approximation { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + std::size_t cod_d = this->cod_d(); + std::size_t amb_d = this->amb_d(); + auto s = tr_.locate_point(p); + Eigen::MatrixXd matrix(cod_d, s.dimension() + 1); + Eigen::MatrixXd vertex_matrix(amb_d + 1, s.dimension() + 1); + for (std::size_t i = 0; i < s.dimension() + 1; ++i) vertex_matrix(0, i) = 1; + std::size_t j = 0; + for (auto v : s.vertex_range()) { + Eigen::VectorXd pt_v = tr_.cartesian_coordinates(v); + Eigen::VectorXd fun_v = fun_(pt_v); + for (std::size_t i = 1; i < amb_d + 1; ++i) vertex_matrix(i, j) = pt_v(i - 1); + for (std::size_t i = 0; i < cod_d; ++i) matrix(i, j) = fun_v(i); + j++; + } + assert(j == s.dimension() + 1); + Eigen::VectorXd z(amb_d + 1); + z(0) = 1; + for (std::size_t i = 1; i < amb_d + 1; ++i) z(i) = p(i - 1); + Eigen::VectorXd lambda = vertex_matrix.colPivHouseholderQr().solve(z); + Eigen::VectorXd result = matrix * lambda; + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief Returns a point on the zero-set. */ + Eigen::VectorXd seed() const { + // TODO: not finished. Should use an oracle. + return Eigen::VectorXd(amb_d()); + } + + /** + * \brief Constructor of the piecewise-linear approximation of a function + * induced by an ambient triangulation. + * + * @param[in] function The function. + * @param[in] triangulation The ambient triangulation. + */ + PL_approximation(const Function_& function, const Triangulation_& triangulation) + : fun_(function), tr_(triangulation) {} + + private: + Function_ fun_; + Triangulation_ tr_; +}; + +/** + * \brief Static constructor of the piecewise-linear approximation of a function + * induced by an ambient triangulation. + * + * @param[in] function The function. + * @param[in] triangulation The ambient triangulation. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_, class Triangulation_> +PL_approximation<Function_, Triangulation_> make_pl_approximation(const Function_& function, + const Triangulation_& triangulation) { + return PL_approximation<Function_, Triangulation_>(function, triangulation); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h new file mode 100644 index 00000000..cbe65abe --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h @@ -0,0 +1,89 @@ +/* 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 FUNCTIONS_TRANSLATE_H_ +#define FUNCTIONS_TRANSLATE_H_ + +#include <cstdlib> // for std::size_t + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Translate + * \brief Translates the zero-set of the function by a vector. + * The underlying function corresponds to f(x-off), where off is the offset vector. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + */ +template <class Function_> +struct Translate { + /** + * \brief Value of the function at a specified point. + * @param[in] p The input point. The dimension needs to coincide with the ambient dimension. + */ + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + Eigen::VectorXd result = fun_(p - off_); + return result; + } + + /** \brief Returns the domain (ambient) dimension. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Returns the codomain dimension. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief Returns a point on the zero-set. */ + Eigen::VectorXd seed() const { + Eigen::VectorXd result = fun_.seed(); + result += off_; + return result; + } + + /** + * \brief Constructor of the translated function. + * + * @param[in] function The function to be translated. + * @param[in] off The offset vector. The dimension should correspond to the + * domain (ambient) dimension of 'function'. + */ + Translate(const Function_& function, const Eigen::VectorXd& off) : fun_(function), off_(off) {} + + private: + Function_ fun_; + Eigen::VectorXd off_; +}; + +/** + * \brief Static constructor of a translated function. + * + * @param[in] function The function to be translated. + * @param[in] off The offset vector. The dimension should correspond to the + * domain (ambient) dimension of 'function'. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Translate<Function_> translate(const Function_& function, Eigen::VectorXd off) { + return Translate<Function_>(function, off); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h new file mode 100644 index 00000000..6a896e94 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h @@ -0,0 +1,72 @@ +/* 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 FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_ +#define FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_ + +#include <cstdlib> // for std::size_t +#include <cmath> // for std::cos, std::sin +#include <random> // for std::uniform_real_distribution, std::random_device + +#include <Eigen/Dense> +#include <Eigen/Sparse> +#include <Eigen/SVD> + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> + +#include <boost/math/constants/constants.hpp> // for PI value + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \brief Generates a uniform random orthogonal matrix using the "subgroup algorithm" by + * Diaconis & Shashahani. + * \details Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Uniform_random_rotation_matrices. + * The idea: take a random rotation matrix of dimension d-1, embed it + * as a d*d matrix M with the last column (0,...,0,1). + * Pick a random vector v on a sphere S^d. rotate the matrix M so that its last column is v. + * The determinant of the matrix can be either 1 or -1 + */ +// Note: the householderQR operation at the end seems to take a lot of time at compilation. +// The CGAL headers are another source of long compilation time. +Eigen::MatrixXd random_orthogonal_matrix(std::size_t d) { + typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; + typedef typename Kernel::Point_d Point_d; + if (d == 1) return Eigen::VectorXd::Constant(1, 1.0); + if (d == 2) { + // 0. < alpha < 2 Pi + std::uniform_real_distribution<double> unif(0., 2 * boost::math::constants::pi<double>()); + std::random_device rand_dev; + std::mt19937 rand_engine(rand_dev()); + double alpha = unif(rand_engine); + + Eigen::Matrix2d rot; + rot << std::cos(alpha), -std::sin(alpha), std::sin(alpha), cos(alpha); + return rot; + } + Eigen::MatrixXd low_dim_rot = random_orthogonal_matrix(d - 1); + Eigen::MatrixXd rot(d, d); + Point_d v = *CGAL::Random_points_on_sphere_d<Point_d>(d, 1); + for (std::size_t i = 0; i < d; ++i) rot(i, 0) = v[i]; + for (std::size_t i = 0; i < d - 1; ++i) + for (std::size_t j = 1; j < d - 1; ++j) rot(i, j) = low_dim_rot(i, j - 1); + for (std::size_t j = 1; j < d; ++j) rot(d - 1, j) = 0; + rot = rot.householderQr() + .householderQ(); // a way to do Gram-Schmidt, see https://forum.kde.org/viewtopic.php?f=74&t=118568#p297246 + return rot; +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h new file mode 100644 index 00000000..ca08f629 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h @@ -0,0 +1,60 @@ +/* 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 IO_MESH_MEDIT_H_ +#define IO_MESH_MEDIT_H_ + +#include <Eigen/Dense> + +#include <vector> +#include <utility> // for std::pair + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Mesh_medit + * \brief Structure to store a mesh that can be output in Medit .mesh file format + * using the output_meshes_to_medit method. + * + * \ingroup coxeter_triangulation + */ +struct Mesh_medit { + /** \brief Type of a range of vertices. */ + typedef std::vector<Eigen::VectorXd> Vertex_points; + /** \brief Type of a mesh element. + * A pair consisting of a vector of vertex indices of type std::size_t + * and of an integer that represents the common reference number for + * the mesh elements of this type. */ + typedef std::pair<std::vector<std::size_t>, std::size_t> Mesh_element; + /** \brief Type of a range of mesh elements. */ + typedef std::vector<Mesh_element> Mesh_elements; + /** \brief Type of a range of scalar field . */ + typedef std::vector<double> Scalar_field_range; + + /** \brief Range of vertices of type Eigen::VectorXd to output. */ + Vertex_points vertex_points; + /** \brief Range of edges. */ + Mesh_elements edges; + /** \brief Range of triangles. */ + Mesh_elements triangles; + /** \brief Range of tetrahedra. */ + Mesh_elements tetrahedra; + /** \brief Range of scalar values over triangles. */ + Scalar_field_range triangles_scalar_range; + /** \brief Range of scalar values over tetrahedra. */ + Scalar_field_range tetrahedra_scalar_range; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h new file mode 100644 index 00000000..9750f366 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h @@ -0,0 +1,171 @@ +/* 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 IO_BUILD_MESH_FROM_CELL_COMPLEX_H_ +#define IO_BUILD_MESH_FROM_CELL_COMPLEX_H_ + +#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES +#include <gudhi/IO/Mesh_medit.h> + +#include <Eigen/Dense> + +#include <cstdlib> // for std::size_t +#include <map> +#include <set> +#include <string> +#include <utility> // for std::make_pair +#include <algorithm> // for std::min + +namespace Gudhi { + +namespace coxeter_triangulation { + +struct Configuration { + Configuration(bool t_edges, bool t_triangles, bool t_tetrahedra, std::size_t r_edges, std::size_t r_triangles, + std::size_t r_tetrahedra) + : toggle_edges(t_edges), + toggle_triangles(t_triangles), + toggle_tetrahedra(t_tetrahedra), + ref_edges(r_edges), + ref_triangles(r_triangles), + ref_tetrahedra(r_tetrahedra) {} + + Configuration() {} + + bool toggle_edges = true, toggle_triangles = true, toggle_tetrahedra = true; + std::size_t ref_edges = 1, ref_triangles = 1, ref_tetrahedra = 1; +}; + +template <class Hasse_cell, class Simplex_cell_map> +void populate_mesh(Mesh_medit& output, Simplex_cell_map& sc_map, Configuration configuration, std::size_t amb_d, + std::map<Hasse_cell*, std::size_t> vi_map) { + using Mesh_element_vertices = Mesh_medit::Mesh_elements::value_type::first_type; + std::map<Hasse_cell*, std::size_t> ci_map; + std::size_t index = vi_map.size() + 1; // current size of output.vertex_points + if (sc_map.size() >= 3) + for (const auto& sc_pair : sc_map[2]) { + Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d); + std::set<std::size_t> vertex_indices; + Hasse_cell* cell = sc_pair.second; + for (const auto& ei_pair : cell->get_boundary()) + for (const auto& vi_pair : ei_pair.first->get_boundary()) vertex_indices.emplace(vi_map[vi_pair.first]); + for (const std::size_t& v : vertex_indices) barycenter += output.vertex_points[v - 1]; + ci_map.emplace(cell, index++); + output.vertex_points.emplace_back((1. / vertex_indices.size()) * barycenter); +#ifdef DEBUG_TRACES + std::string vlist = " (" + std::to_string(index - 1) + ")"; + for (const std::size_t& v : vertex_indices) vlist += " " + std::to_string(v); + cell_vlist_map.emplace(to_string(cell), vlist); +#endif + } + + if (configuration.toggle_edges && sc_map.size() >= 2) + for (const auto& sc_pair : sc_map[1]) { + Hasse_cell* edge_cell = sc_pair.second; + Mesh_element_vertices edge; + for (const auto& vi_pair : edge_cell->get_boundary()) edge.push_back(vi_map[vi_pair.first]); + output.edges.emplace_back(edge, configuration.ref_edges); +#ifdef DEBUG_TRACES + std::string vlist; + for (const std::size_t& v : edge) vlist += " " + std::to_string(v); + cell_vlist_map.emplace(to_string(edge_cell), vlist); +#endif + } + + if (configuration.toggle_triangles && sc_map.size() >= 3) + for (const auto& sc_pair : sc_map[2]) { + for (const auto& ei_pair : sc_pair.second->get_boundary()) { + Mesh_element_vertices triangle(1, ci_map[sc_pair.second]); + for (const auto& vi_pair : ei_pair.first->get_boundary()) triangle.push_back(vi_map[vi_pair.first]); + output.triangles.emplace_back(triangle, configuration.ref_triangles); + } + } + + if (configuration.toggle_tetrahedra && sc_map.size() >= 4) + for (const auto& sc_pair : sc_map[3]) { + Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d); + std::set<std::size_t> vertex_indices; + Hasse_cell* cell = sc_pair.second; + for (const auto& ci_pair : cell->get_boundary()) + for (const auto& ei_pair : ci_pair.first->get_boundary()) + for (const auto& vi_pair : ei_pair.first->get_boundary()) vertex_indices.emplace(vi_map[vi_pair.first]); + for (const std::size_t& v : vertex_indices) barycenter += output.vertex_points[v - 1]; + output.vertex_points.emplace_back((1. / vertex_indices.size()) * barycenter); +#ifdef DEBUG_TRACES + std::string vlist = " (" + std::to_string(index) + ")"; + for (const std::size_t& v : vertex_indices) vlist += " " + std::to_string(v); + cell_vlist_map.emplace(to_string(cell), vlist); +#endif + + for (const auto& ci_pair : cell->get_boundary()) + for (const auto& ei_pair : ci_pair.first->get_boundary()) { + Mesh_element_vertices tetrahedron = {index, ci_map[sc_pair.second]}; + for (const auto& vi_pair : ei_pair.first->get_boundary()) tetrahedron.push_back(vi_map[vi_pair.first]); + output.tetrahedra.emplace_back(tetrahedron, configuration.ref_tetrahedra); + } + index++; + } +} + +/** @brief Builds a Gudhi::coxeter_triangulation::Mesh_medit from a Gudhi::coxeter_triangulation::Cell_complex + * + * @ingroup coxeter_triangulation + */ +template <class Cell_complex> +Mesh_medit build_mesh_from_cell_complex(const Cell_complex& cell_complex, + Configuration i_configuration = Configuration(), + Configuration b_configuration = Configuration()) { + using Hasse_cell = typename Cell_complex::Hasse_cell; + Mesh_medit output; + std::map<Hasse_cell*, std::size_t> vi_map; // one for vertices, other for 2d-cells + std::size_t index = 1; // current size of output.vertex_points + + if (cell_complex.cell_point_map().empty()) return output; + std::size_t amb_d = std::min((int)cell_complex.cell_point_map().begin()->second.size(), 3); + + for (const auto& cp_pair : cell_complex.cell_point_map()) { +#ifdef DEBUG_TRACES + std::string vlist; + vlist += " " + std::to_string(index); + cell_vlist_map.emplace(to_string(cp_pair.first), vlist); +#endif + vi_map.emplace(cp_pair.first, index++); + output.vertex_points.push_back(cp_pair.second); + output.vertex_points.back().conservativeResize(amb_d); + } + + populate_mesh(output, cell_complex.interior_simplex_cell_maps(), i_configuration, amb_d, vi_map); +#ifdef DEBUG_TRACES + for (const auto& sc_map : cell_complex.interior_simplex_cell_maps()) + for (const auto& sc_pair : sc_map) { + std::string simplex = "I" + to_string(sc_pair.first); + std::string cell = to_string(sc_pair.second); + std::string vlist = cell_vlist_map.at(cell).substr(1); + simplex_vlist_map.emplace(simplex, vlist); + } +#endif + populate_mesh(output, cell_complex.boundary_simplex_cell_maps(), b_configuration, amb_d, vi_map); +#ifdef DEBUG_TRACES + for (const auto& sc_map : cell_complex.boundary_simplex_cell_maps()) + for (const auto& sc_pair : sc_map) { + std::string simplex = "B" + to_string(sc_pair.first); + std::string cell = to_string(sc_pair.second); + std::string vlist = cell_vlist_map.at(cell).substr(1); + simplex_vlist_map.emplace(simplex, vlist); + } +#endif + return output; +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h new file mode 100644 index 00000000..a2995738 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h @@ -0,0 +1,550 @@ +#ifndef IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_ +#define IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_ + +#ifdef DEBUG_TRACES // All this part of code can be skipped if DEBUG_TRACES are not ON - cmake -DDEBUG_TRACES=ON . + +#include <sstream> +#include <fstream> +#include <vector> +#include <list> +#include <string> +#include <regex> + +#include <Eigen/Dense> + +namespace Gudhi { + +namespace coxeter_triangulation { + +template <class T> +std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector) { + os << "("; + if (vector.empty()) { + os << ")"; + return os; + } + auto v_it = vector.begin(); + os << *v_it++; + for (; v_it != vector.end(); ++v_it) os << ", " << *v_it; + os << ")"; + return os; +} + +/* A class to make the vector horizontal instead of vertical */ +struct Straighten { + Straighten(const Eigen::VectorXd& vector) : vector_(vector) {} + const Eigen::VectorXd& vector_; +}; + +std::ostream& operator<<(std::ostream& os, const Straighten& str) { + std::size_t size = str.vector_.size(); + os << "(" << str.vector_(0); + if (size == 0) { + os << ")"; + return os; + } + for (std::size_t i = 1; i < size; ++i) os << ", " << str.vector_(i); + os << ")"; + return os; +} + +std::string id_from_simplex(const std::string& simplex) { + std::regex r("\\s+"), r2("\\(|\\)|\\{|\\}"), r3(","), r4("\\["), r5("\\]"); + std::string output = std::regex_replace(simplex, r, ""); + output = std::regex_replace(output, r2, ":"); + output = std::regex_replace(output, r3, "."); + output = std::regex_replace(output, r4, "_"); + output = std::regex_replace(output, r5, ""); + return output; +} + +template <typename T> +std::string to_string(const T& t) { + std::ostringstream oss; + oss << t; + return oss.str(); +} + +struct MT_inserted_info { + std::string qr_face_, init_face_, qr_intersection_; + bool qr_success_, is_boundary_; + template <class Query_result, class Simplex_handle> + MT_inserted_info(const Query_result& qr, const Simplex_handle& face, bool is_boundary) + : qr_face_(to_string(face)), + init_face_(to_string(face)), + qr_intersection_(to_string(qr.intersection)), + qr_success_(qr.success), + is_boundary_(is_boundary) {} +}; +std::list<MT_inserted_info> mt_seed_inserted_list, mt_inserted_list; + +struct CC_summary_info { + std::string face_, cell_; + template <class SC_pair> + CC_summary_info(const SC_pair& sc_pair) : face_(to_string(sc_pair.first)), cell_(to_string(sc_pair.second)) {} +}; +using CC_summary_list = std::list<CC_summary_info>; +std::vector<CC_summary_list> cc_interior_summary_lists, cc_boundary_summary_lists; + +struct CC_detail_info { + enum class Result_type { self, face, coface, inserted, join_single, join_is_face }; + std::string simplex_, trigger_, init_simplex_; + Result_type status_; + bool join_trigger_ = false; + std::list<std::string> faces_, post_faces_, cofaces_; + template <class Simplex_handle> + CC_detail_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)) {} +}; +using CC_detail_list = std::list<CC_detail_info>; +std::vector<CC_detail_list> cc_interior_detail_lists, cc_boundary_detail_lists; +std::vector<CC_detail_list> cc_interior_insert_detail_lists, cc_boundary_insert_detail_lists; + +struct CC_prejoin_info { + enum class Result_type { join_single, join_is_face, join_different, join_same }; + std::string simplex_, join_; + std::vector<std::string> faces_; + std::size_t dimension_; + Result_type status_; + template <class Simplex_handle> + CC_prejoin_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)), dimension_(simplex.dimension()) {} +}; +using CC_prejoin_list = std::list<CC_prejoin_info>; +std::vector<CC_prejoin_list> cc_interior_prejoin_lists, cc_boundary_prejoin_lists; + +struct CC_join_info { + enum class Result_type { self, face, coface, inserted, join_single, join_is_face }; + std::string simplex_, join_, trigger_; + Result_type status_; + std::list<std::string> boundary_faces_; + std::list<std::string> faces_, post_faces_, cofaces_; + template <class Simplex_handle> + CC_join_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)) {} +}; +bool join_switch = false; +std::vector<CC_detail_list> cc_interior_join_detail_lists, cc_boundary_join_detail_lists; + +std::map<std::string, std::string> cell_vlist_map; +std::map<std::string, std::string> simplex_vlist_map; + +std::ostringstream mt_ostream, vis_ostream; +std::vector<std::ostringstream> cc_summary_ostream, cc_traces_ostream; + +std::string simplex_format(const std::string& simplex, bool is_boundary) { + std::string b_simplex = (is_boundary ? "B" : "I") + simplex; + std::string tooltiptext; + auto it = simplex_vlist_map.find(b_simplex); + if (it == simplex_vlist_map.end()) + tooltiptext = "deleted"; + else + tooltiptext = simplex_vlist_map.at(b_simplex); + return (std::string) "<a class=\"" + (is_boundary ? "boundary" : "interior") + "\" href=\"#" + + id_from_simplex(b_simplex) + "\">" + b_simplex + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>"; +} + +std::string simplex_format(const std::string& b_simplex) { + bool is_boundary = b_simplex[0] == 'B'; + std::string tooltiptext; + auto it = simplex_vlist_map.find(b_simplex); + if (it == simplex_vlist_map.end()) + tooltiptext = "deleted"; + else + tooltiptext = simplex_vlist_map.at(b_simplex); + return (std::string) "<a class=\"" + (is_boundary ? "boundary" : "interior") + "\" href=\"#" + + id_from_simplex(b_simplex) + "\">" + b_simplex + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>"; +} + +void write_head(std::ofstream& ofs) { + ofs << " <head>\n" + << " <title>Cell complex debug trace</title>\n" + << " <style>\n" + << " a.boundary {\n" + << " position: relative;\n" + << " display: inline-block;\n" + << " color: darkred;\n" + << " background-color: lightgreen\n" + << " }\n" + << " a.interior {\n" + << " position: relative;\n" + << " display: inline-block;\n" + << " color: navy;\n" + << " background-color: yellow\n" + << " }\n" + << " .tooltiptext {\n" + << " visibility: hidden;\n" + << " width: 120px;\n" + << " background-color: #555;\n" + << " color: #fff;\n" + << " text-align: center;\n" + << " padding: 5px 0;\n" + << " border-radius: 6px;\n" + << " position: absolute;\n" + << " z-index: 1;\n" + << " bottom: 125%;\n" + << " left: 50%;\n" + << " margin-left: -60px;\n" + << " opacity: 0;\n" + << " transition: opacity 0.3s;\n" + << " }\n" + << " .boundary .tooltiptext::after {\n" + << " content: \"\";\n" + << " position: absolute;\n" + << " top: 100%;\n" + << " left: 50%;\n" + << " margin-left: -5px;\n" + << " border-width: 5px;\n" + << " border-style: solid;\n" + << " border-color: #555 transparent transparent transparent;\n" + << " }\n" + << " .interior .tooltiptext::after {\n" + << " content: \"\";\n" + << " position: absolute;\n" + << " top: 100%;\n" + << " left: 50%;\n" + << " margin-left: -5px;\n" + << " border-width: 5px;\n" + << " border-style: solid;\n" + << " border-color: #555 transparent transparent transparent;\n" + << " }\n" + << " .boundary:hover .tooltiptext {\n" + << " visibility: visible;\n" + << " opacity: 1;\n" + << " }\n" + << " .interior:hover .tooltiptext {\n" + << " visibility: visible;\n" + << " opacity: 1;\n" + << " }\n" + << " ul.nav {\n" + << " list-style-type: none;\n" + << " margin: 0;\n" + << " padding: 0;\n" + << " overflow: auto;\n" + << " background-color: #333;\n" + << " position: fixed;\n" + << " height: 100%;\n" + << " width: 15%;\n" + << " }\n" + << " ul.nav li a {\n" + << " display: block;\n" + << " color: white;\n" + << " text-align: left;\n" + << " padding: 14px 16px;\n" + << " text-decoration: none;\n" + << " }\n" + << " .active {\n" + << " background-color: #4CAF50;\n" + << " }\n" + << " div {\n" + << " margin-left: 15%;\n" + << " padding: 1px 16px\n" + << " }\n" + << " div.navi {\n" + << " margin-left: 0%;\n" + << " padding: 0px 0px\n" + << " }\n" + << " h1 {\n" + << " margin-left: 15%;\n" + << " padding: 1px 16px\n" + << " }\n" + << " </style>\n" + << " </head>\n"; +} + +void write_nav(std::ofstream& ofs) { + ofs << " <div class=\"navi\" style=\"margin-top:30px;background-color:#1abc9c;\">\n" + << " <ul class=\"nav\">\n" + << " <li><a href=\"#mant\">Manifold tracing</a></li>\n" + << " <li><a href=\"#cell\">Cell complex</a>\n" + << " <ul>\n"; + for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) { + ofs << " <li><a href=\"#dim" << i << "\">Dimension " << i << "</a>\n"; + ofs << " <ul>\n"; + ofs << " <li><a href=\"#dim" << i << "i\">Interior</a></li>\n"; + if (i < cc_boundary_summary_lists.size()) { + ofs << " <li><a href=\"#dim" << i << "b\">Boundary</a></li>\n"; + } + ofs << " </ul>\n"; + ofs << " </li>\n"; + } + ofs << " </ul>\n" + << " </li>\n" + << " <li><a href=\"#visu\">Visualization details</a></li>\n" + << " </ul>\n" + << " </div>\n"; +} + +void write_mt(std::ofstream& ofs) { + ofs << " <div id=\"mant\">\n"; + ofs << " <h2> Manifold debug trace </h2>\n"; + ofs << " <h3> Simplices inserted during the seed phase </h3>\n"; + ofs << " <ul>\n"; + for (const MT_inserted_info& mt_info : mt_seed_inserted_list) { + if (mt_info.qr_success_) { + ofs << " <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_); + if (mt_info.qr_face_ != mt_info.init_face_) + ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")"; + ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n"; + } else + ofs << " <li>Failed to insert " << mt_info.init_face_ << "</li>\n"; + } + ofs << " </ul>\n"; + ofs << " <h3> Simplices inserted during the while loop phase </h3>\n"; + ofs << " <ul>\n"; + for (const MT_inserted_info& mt_info : mt_inserted_list) { + if (mt_info.qr_success_) { + ofs << " <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_); + if (mt_info.qr_face_ != mt_info.init_face_) + ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")"; + ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n"; + } else + ofs << " <li>Failed to insert " << mt_info.init_face_ << ")</li>\n"; + } + ofs << " </ul>\n"; + ofs << " </div>\n"; +} + +void write_cc(std::ofstream& ofs) { + ofs << " <div id=\"cell\">\n" + << " <h2> Cell complex debug trace </h2>\n" + << " <p>Go to:</p>\n" + << " <ul>\n"; + for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) { + ofs << " <li><a href=\"#dim" << i << "\">Dimension " << i << "</a></li>\n"; + } + ofs << " </ul>\n"; + for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) { + ofs << " <h3 id=\"dim" << i << "\"> Dimension " << i << "</h3>\n"; + ofs << " <h4 id=\"dim" << i << "i\"> Summary for interior simplices</h4>\n"; + if (i < cc_boundary_summary_lists.size()) ofs << " <p><a href=\"#dim" << i << "b\">Go to boundary</a></p>\n"; + ofs << " <ul>\n"; + for (const CC_summary_info& cc_info : cc_interior_summary_lists[i]) + ofs << " <li id = \"" << id_from_simplex("I" + cc_info.face_) << "\">" + << simplex_format(cc_info.face_, false) << " cell =" << cc_info.cell_ << "</li>\n"; + ofs << " </ul>\n"; + ofs << " <h4> Prejoin state of the interior cells of dimension " << i << "</h4>\n"; + auto prejoin_it = cc_interior_prejoin_lists[i].begin(); + while (prejoin_it != cc_interior_prejoin_lists[i].end()) { + std::size_t j = prejoin_it->dimension_; + ofs << " <h5>" << j << "-dimensional ambient simplices</h5>\n"; + ofs << " <ul>\n"; + for (; prejoin_it->dimension_ == j; ++prejoin_it) { + ofs << " <li>" << simplex_format(prejoin_it->simplex_, false) + << " join = " << simplex_format(prejoin_it->join_, false) << " boundary:\n" + << " <ul>\n"; + for (const auto& face : prejoin_it->faces_) ofs << " <li>" << simplex_format(face) << "</li>"; + ofs << " </ul>\n"; + switch (prejoin_it->status_) { + case (CC_prejoin_info::Result_type::join_single): + ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, false) + << " as it has a single face.</p>"; + break; + case (CC_prejoin_info::Result_type::join_is_face): + ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, false) + << " as its join " << simplex_format(prejoin_it->join_, false) << " is one of the faces.</p>"; + break; + case (CC_prejoin_info::Result_type::join_different): + ofs << " <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, false) + << " and replaced by its join " << simplex_format(prejoin_it->join_, false) << ".</p>"; + break; + case (CC_prejoin_info::Result_type::join_same): + ofs << " <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, false) + << ".</p>"; + } + ofs << " </li>"; + } + ofs << " </ul>\n"; + } + ofs << " <h4> Details for interior simplices</h4>\n"; + ofs << " <ul>\n"; + for (const CC_detail_info& cc_info : cc_interior_detail_lists[i]) { + if (cc_info.status_ == CC_detail_info::Result_type::join_single) { + ofs << " <li style=\"color:magenta\" id = \"" << id_from_simplex("I" + cc_info.simplex_) + << "\"> Simplex " << simplex_format(cc_info.simplex_, false) << " has only one face (" + << simplex_format(cc_info.trigger_, false) << ") and is deleted."; + continue; + } + if (cc_info.status_ == CC_detail_info::Result_type::join_single) { + ofs << " <li style=\"color:darkmagenta\" id = \"" << id_from_simplex("I" + cc_info.simplex_) + << "\"> The join of the simplex " << simplex_format(cc_info.simplex_, false) << " is one of its faces (" + << simplex_format(cc_info.trigger_, false) << "), hence it is is deleted."; + continue; + } + ofs << " <li> Insert_cell called for " << simplex_format(cc_info.simplex_, false) << "\n"; + ofs << " <ul>\n"; + // for (const std::string& cof: cc_info.faces_) + // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, false) + // << " is a face of " << simplex_format(cof, false) << "\n"; + ofs << " </ul>\n"; + ofs << " <ul>\n"; + if (cc_info.status_ == CC_detail_info::Result_type::self) { + ofs << " <p><span style=\"color:blue\">The simplex " << simplex_format(cc_info.simplex_, false) + << " already exists in the cell complex!</span></p>\n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::face) { + ofs << " <p><span style=\"color:red\">The simplex " << simplex_format(cc_info.simplex_, false) + << " is a face of the simplex " << simplex_format(cc_info.trigger_, false) << "!</span><br>\n"; + ofs << " <ul>\n"; + for (const std::string post_face : cc_info.post_faces_) + ofs << " <li id = \"" << id_from_simplex("I" + post_face) << "\">" + << "Post deleting " << simplex_format(post_face, false) << "</li>\n"; + ofs << " </ul>\n"; + ofs << " </p>\n"; + ofs << " <p id = \"" << id_from_simplex("I" + cc_info.trigger_) << "\">" + << "Deleting " << simplex_format(cc_info.trigger_, false) << "</p>\n"; + } + // for (const std::string& fac: cc_info.cofaces_) + // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, false) + // << " is a coface of " << simplex_format(fac, false) << "\n"; + if (cc_info.status_ == CC_detail_info::Result_type::coface) { + ofs << " <p><span style=\"color:darkorange\">The simplex " << simplex_format(cc_info.simplex_, false) + << " is a coface of the simplex " << simplex_format(cc_info.trigger_, false) << "!</span><p>\n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::inserted) { + ofs << " <p><span style=\"color:green\">Successfully inserted " + << simplex_format(cc_info.simplex_, false) << "!</span><p>\n"; + } + ofs << " </ul>\n"; + ofs << " </li>\n"; + } + ofs << " </ul>\n"; + + if (i < cc_boundary_summary_lists.size()) { + ofs << " <h4 id=\"dim" << i << "b\"> Summary for boundary simplices</h4>\n"; + ofs << " <p><a href=\"#dim" << i << "i\">Go to interior</a></p>\n"; + ofs << " <ul>\n"; + for (const CC_summary_info& cc_info : cc_boundary_summary_lists[i]) + ofs << " <li id = \"" << id_from_simplex("B" + cc_info.face_) << "\">" + << simplex_format(cc_info.face_, true) << " cell =" << cc_info.cell_ << "</li>\n"; + ofs << " </ul>\n"; + ofs << " <h4> Prejoin state of the boundary cells of dimension " << i << "</h4>\n"; + auto prejoin_it = cc_boundary_prejoin_lists[i].begin(); + while (prejoin_it != cc_boundary_prejoin_lists[i].end()) { + std::size_t j = prejoin_it->dimension_; + ofs << " <h5>" << j << "-dimensional ambient simplices</h5>\n"; + ofs << " <ul>\n"; + for (; prejoin_it->dimension_ == j; ++prejoin_it) { + ofs << " <li>" << simplex_format(prejoin_it->simplex_, true) + << " join = " << simplex_format(prejoin_it->join_, true) << " boundary:\n" + << " <ul>\n"; + for (const auto& face : prejoin_it->faces_) ofs << " <li>" << simplex_format(face) << "</li>"; + ofs << " </ul>\n"; + switch (prejoin_it->status_) { + case (CC_prejoin_info::Result_type::join_single): + ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, true) + << " as it has a single face.</p>"; + break; + case (CC_prejoin_info::Result_type::join_is_face): + ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, true) + << " as its join " << simplex_format(prejoin_it->join_, true) << " is one of the faces.</p>"; + break; + case (CC_prejoin_info::Result_type::join_different): + ofs << " <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, true) + << " and replaced by its join " << simplex_format(prejoin_it->join_, true) << ".</p>"; + break; + case (CC_prejoin_info::Result_type::join_same): + ofs << " <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, true) + << ".</p>"; + } + ofs << " </li>"; + } + ofs << " </ul>\n"; + } + } + if (i < cc_boundary_detail_lists.size()) { + ofs << " <h4> Details for boundary simplices</h4>\n" + << " <ul>\n"; + for (const CC_detail_info& cc_info : cc_boundary_detail_lists[i]) { + if (cc_info.status_ == CC_detail_info::Result_type::join_single) { + ofs << " <li style=\"color:magenta\" id = \"" << id_from_simplex("B" + cc_info.simplex_) + << "\"> Simplex " << simplex_format(cc_info.simplex_, true) << " has only one face (" + << simplex_format(cc_info.trigger_, true) << ") and is deleted."; + continue; + } + if (cc_info.status_ == CC_detail_info::Result_type::join_single) { + ofs << " <li style=\"color:darkmagenta\" id = \"" << id_from_simplex("B" + cc_info.simplex_) + << "\"> The join of the simplex " << simplex_format(cc_info.simplex_, true) << " is one of its faces (" + << simplex_format(cc_info.trigger_, true) << "), hence it is is deleted."; + continue; + } + ofs << " <li> Insert_simplex called on " << simplex_format(cc_info.simplex_, true); + ofs << " <ul>\n"; + // for (const std::string& cof: cc_info.faces_) + // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, true) + // << " is a face of " << simplex_format(cof, true) << "\n"; + ofs << " </ul>\n"; + ofs << " <ul>\n"; + if (cc_info.status_ == CC_detail_info::Result_type::self) { + ofs << " <p><span style=\"color:blue\">The simplex " << simplex_format(cc_info.simplex_, true) + << " already exists in the cell complex!</span></p>\n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::face) { + ofs << " <p><span style=\"color:red\">The simplex " << simplex_format(cc_info.simplex_, true) + << " is a face of the simplex " << simplex_format(cc_info.trigger_, true) << "!</span><br>\n"; + ofs << " <ul>\n"; + for (const std::string post_face : cc_info.post_faces_) + ofs << " <li id=\"" << id_from_simplex("B" + post_face) << "\">Post deleting " + << simplex_format(post_face, true) << "</li>\n"; + ofs << " </ul>\n"; + ofs << " </p>\n"; + ofs << " <p id=\"" << id_from_simplex(cc_info.trigger_) << "\">Deleting " + << simplex_format(cc_info.trigger_, true) << "</p>\n"; + } + // for (const std::string& fac: cc_info.cofaces_) + // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, true) + // << " is a coface of " << simplex_format(fac, true) << "\n"; + ofs << " </ul>\n"; + ofs << " </li>\n"; + if (cc_info.status_ == CC_detail_info::Result_type::coface) { + ofs << " <p><span style=\"color:darkorange\">The simplex " + << simplex_format(cc_info.simplex_, true) << " is a coface of the simplex " + << simplex_format(cc_info.trigger_, true) << "!</span><p>\n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::inserted) { + ofs << " <p><span style=\"color:green\">Successfully inserted " + << simplex_format(cc_info.simplex_, true) << "!</span><p>\n"; + } + } + ofs << " </ul>\n"; + } + } + ofs << " </div>\n"; +} + +void write_visu(std::ofstream& ofs) { + ofs << " <div id=\"visu\">\n" + << " <h2> Visualization details debug trace </h2>\n"; + // std::vector<std::map<std::string, std::string> > vs_maps(cc_interior_summary_lists.size()); + std::map<std::string, std::string> vs_map; + for (const auto& sv_pair : simplex_vlist_map) vs_map.emplace(sv_pair.second, sv_pair.first); + ofs << " <ul>\n"; + for (const auto& vs_pair : vs_map) { + std::string w_simplex = vs_pair.second.substr(1); + bool is_boundary = vs_pair.second[0] == 'B'; + ofs << " <li><b>" << vs_pair.first << "</b>: " << simplex_format(w_simplex, is_boundary) << "</li>\n"; + } + ofs << " </ul>\n"; + ofs << " </div>\n"; +} + +void write_to_html(std::string file_name) { + std::ofstream ofs(file_name + ".html", std::ofstream::out); + ofs << "<!DOCTYPE html>\n" + << "<html>\n"; + write_head(ofs); + ofs << " <body>\n"; + write_nav(ofs); + ofs << " <h1> Debug traces for " << file_name << " </h1>\n"; + write_mt(ofs); + write_cc(ofs); + write_visu(ofs); + ofs << " </body>\n"; + ofs << "</html>\n"; + + ofs.close(); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif // DEBUG_TRACES +#endif // IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_ diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h new file mode 100644 index 00000000..f69d8b29 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h @@ -0,0 +1,154 @@ +/* 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 IO_OUTPUT_MESHES_TO_MEDIT_H_ +#define IO_OUTPUT_MESHES_TO_MEDIT_H_ + +#include <gudhi/IO/Mesh_medit.h> + +#include <Eigen/Dense> + +#include <cstdlib> // for std::size_t +#include <fstream> // for std::ofstream +#include <vector> +#include <type_traits> // for std::enable_if +#include <tuple> // for std::get +#include <utility> // for std::make_pair + +namespace Gudhi { + +namespace coxeter_triangulation { + +using Vertex_points = Mesh_medit::Vertex_points; +using Mesh_elements = Mesh_medit::Mesh_elements; +using Scalar_field_range = Mesh_medit::Scalar_field_range; + +template <std::size_t I = 0, typename... Meshes> +typename std::enable_if<I == sizeof...(Meshes), void>::type fill_meshes(Vertex_points& vertex_points, + Mesh_elements& edges, Mesh_elements& triangles, + Mesh_elements& tetrahedra, + Scalar_field_range& triangles_scalar_range, + Scalar_field_range& tetrahedra_scalar_range, + std::size_t index, const Meshes&... meshes) {} + +template <std::size_t I = 0, typename... Meshes> +typename std::enable_if<I != sizeof...(Meshes), void>::type fill_meshes(Vertex_points& vertex_points, + Mesh_elements& edges, Mesh_elements& triangles, + Mesh_elements& tetrahedra, + Scalar_field_range& triangles_scalar_range, + Scalar_field_range& tetrahedra_scalar_range, + std::size_t index, const Meshes&... meshes) { + auto mesh = std::get<I>(std::forward_as_tuple(meshes...)); + for (const auto& v : mesh.vertex_points) vertex_points.push_back(v); + for (const auto& e : mesh.edges) { + std::vector<std::size_t> edge; + for (const auto& v_i : e.first) edge.push_back(v_i + index); + edges.emplace_back(edge, e.second); + } + for (const auto& t : mesh.triangles) { + std::vector<std::size_t> triangle; + for (const auto& v_i : t.first) triangle.push_back(v_i + index); + triangles.emplace_back(triangle, t.second); + } + for (const auto& t : mesh.tetrahedra) { + std::vector<std::size_t> tetrahedron; + for (const auto& v_i : t.first) tetrahedron.push_back(v_i + index); + tetrahedra.emplace_back(tetrahedron, t.second); + } + for (const auto& b : mesh.triangles_scalar_range) triangles_scalar_range.push_back(b); + for (const auto& b : mesh.tetrahedra_scalar_range) tetrahedra_scalar_range.push_back(b); + fill_meshes<I + 1, Meshes...>(vertex_points, edges, triangles, tetrahedra, triangles_scalar_range, + tetrahedra_scalar_range, index + mesh.vertex_points.size(), meshes...); +} + +/** \brief Outputs a text file with specified meshes that can be visualized in + * <a target="_blank" href="https://www.ljll.math.upmc.fr/frey/software.html">Medit</a>. + * + * @param[in] amb_d Ambient dimension. Can be 2 or 3. + * @param[in] file_name The name of the output file. + * @param[in] meshes A pack of meshes to be specified separated by commas. + * + * @ingroup coxeter_triangulation + */ +template <typename... Meshes> +void output_meshes_to_medit(std::size_t amb_d, std::string file_name, const Meshes&... meshes) { + Vertex_points vertex_points; + Mesh_elements edges, triangles, tetrahedra; + Scalar_field_range triangles_scalar_range, tetrahedra_scalar_range; + fill_meshes(vertex_points, edges, triangles, tetrahedra, triangles_scalar_range, tetrahedra_scalar_range, 0, + meshes...); + + std::ofstream ofs(file_name + ".mesh", std::ofstream::out); + std::ofstream ofs_bb(file_name + ".bb", std::ofstream::out); + + if (amb_d == 2) { + ofs << "MeshVersionFormatted 1\nDimension 2\n"; + ofs_bb << "2 1 "; + ofs << "Vertices\n" << vertex_points.size() << "\n"; + for (auto p : vertex_points) { + ofs << p[0] << " " << p[1] << " 2\n"; + } + ofs << "Edges " << edges.size() << "\n"; + for (auto e : edges) { + for (auto v : e.first) ofs << v << " "; + ofs << e.second << std::endl; + } + ofs << "Triangles " << triangles.size() << "\n"; + for (auto s : triangles) { + for (auto v : s.first) { + ofs << v << " "; + } + ofs << s.second << std::endl; + } + + ofs_bb << triangles_scalar_range.size() << " 1\n"; + for (auto& b : triangles_scalar_range) ofs_bb << b << "\n"; + + } else { + ofs << "MeshVersionFormatted 1\nDimension 3\n"; + ofs_bb << "3 1 "; + ofs << "Vertices\n" << vertex_points.size() << "\n"; + for (auto p : vertex_points) { + ofs << p[0] << " " << p[1] << " " << p[2] << " 2\n"; + } + ofs << "Edges " << edges.size() << "\n"; + for (auto e : edges) { + for (auto v : e.first) ofs << v << " "; + ofs << e.second << std::endl; + } + ofs << "Triangles " << triangles.size() << "\n"; + for (auto s : triangles) { + for (auto v : s.first) { + ofs << v << " "; + } + ofs << s.second << std::endl; + } + ofs << "Tetrahedra " << tetrahedra.size() << "\n"; + for (auto s : tetrahedra) { + for (auto v : s.first) { + ofs << v << " "; + } + ofs << s.second << std::endl; + } + + ofs_bb << triangles_scalar_range.size() + tetrahedra_scalar_range.size() << " 1\n"; + for (auto& b : triangles_scalar_range) ofs_bb << b << "\n"; + for (auto& b : tetrahedra_scalar_range) ofs_bb << b << "\n"; + } + + ofs.close(); + ofs_bb.close(); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h new file mode 100644 index 00000000..277f8b6c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h @@ -0,0 +1,261 @@ +/* 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 IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_ +#define IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_ + +#include <Eigen/Dense> + +#include <gudhi/Permutahedral_representation/face_from_indices.h> +#include <gudhi/Functions/Constant_function.h> +#include <gudhi/Functions/PL_approximation.h> +#include <gudhi/Coxeter_triangulation/Query_result.h> +#include <gudhi/Debug_utils.h> // for GUDHI_CHECK + +#include <vector> +#include <limits> // for std::numeric_limits<> +#include <cmath> // for std::fabs + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Implicit_manifold_intersection_oracle + * \brief An oracle that supports the intersection query on an implicit manifold. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * \tparam Domain_function_ The domain function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_, class Domain_function_ = Constant_function> +class Implicit_manifold_intersection_oracle { + /* Computes the affine coordinates of the intersection point of the implicit manifold + * and the affine hull of the simplex. */ + template <class Simplex_handle, class Triangulation> + Eigen::VectorXd compute_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const { + std::size_t cod_d = this->cod_d(); + Eigen::MatrixXd matrix(cod_d + 1, cod_d + 1); + for (std::size_t i = 0; i < cod_d + 1; ++i) matrix(0, i) = 1; + std::size_t j = 0; + for (auto v : simplex.vertex_range()) { + Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v)); + for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1); + j++; + } + Eigen::VectorXd z(cod_d + 1); + z(0) = 1; + for (std::size_t i = 1; i < cod_d + 1; ++i) z(i) = 0; + Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z); + if (!z.isApprox(matrix*lambda)) { + // NaN non valid results + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) = + std::numeric_limits<double>::quiet_NaN(); + } + return lambda; + } + + /* Computes the affine coordinates of the intersection point of the boundary + * of the implicit manifold and the affine hull of the simplex. */ + template <class Simplex_handle, class Triangulation> + Eigen::VectorXd compute_boundary_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const { + std::size_t cod_d = this->cod_d(); + Eigen::MatrixXd matrix(cod_d + 2, cod_d + 2); + for (std::size_t i = 0; i < cod_d + 2; ++i) matrix(0, i) = 1; + std::size_t j = 0; + for (auto v : simplex.vertex_range()) { + Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v)); + for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1); + Eigen::VectorXd bv_coords = domain_fun_(triangulation.cartesian_coordinates(v)); + matrix(cod_d + 1, j) = bv_coords(0); + j++; + } + Eigen::VectorXd z(cod_d + 2); + z(0) = 1; + for (std::size_t i = 1; i < cod_d + 2; ++i) z(i) = 0; + Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z); + if (!z.isApprox(matrix*lambda)) { + // NaN non valid results + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) = + std::numeric_limits<double>::quiet_NaN(); + } + return lambda; + } + + /* Computes the intersection result for a given simplex in a triangulation. */ + template <class Simplex_handle, class Triangulation> + Query_result<Simplex_handle> intersection_result(const Eigen::VectorXd& lambda, const Simplex_handle& simplex, + const Triangulation& triangulation) const { + using QR = Query_result<Simplex_handle>; + std::size_t amb_d = triangulation.dimension(); + std::size_t cod_d = simplex.dimension(); + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) { + if (std::isnan(lambda(i))) return QR({Eigen::VectorXd(), false}); + GUDHI_CHECK((std::fabs(lambda(i) - 1.) > std::numeric_limits<double>::epsilon() && + std::fabs(lambda(i) - 0.) > std::numeric_limits<double>::epsilon()), + std::invalid_argument("A vertex of the triangulation lies exactly on the manifold")); + if (lambda(i) < 0. || lambda(i) > 1.) return QR({Eigen::VectorXd(), false}); + } + Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d); + auto v_range = simplex.vertex_range(); + auto v_it = v_range.begin(); + for (std::size_t i = 0; i < cod_d + 1 && v_it != v_range.end(); ++v_it, ++i) { + Eigen::VectorXd v_coords = triangulation.cartesian_coordinates(*v_it); + for (std::size_t j = 0; j < amb_d; ++j) vertex_matrix(i, j) = v_coords(j); + } + Eigen::VectorXd intersection = lambda.transpose() * vertex_matrix; + return QR({intersection, true}); + } + + public: + /** \brief Ambient dimension of the implicit manifold. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Codimension of the implicit manifold. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief The seed point of the implicit manifold. */ + Eigen::VectorXd seed() const { return fun_.seed(); } + + /** \brief Intersection query with the relative interior of the manifold. + * + * \details The returned structure Query_result contains the boolean value + * that is true only if the intersection point of the query simplex and + * the relative interior of the manifold exists, the intersection point + * and the face of the query simplex that contains + * the intersection point. + * + * \tparam Simplex_handle The class of the query simplex. + * Needs to be a model of the concept SimplexInCoxeterTriangulation. + * \tparam Triangulation The class of the triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * + * @param[in] simplex The query simplex. The dimension of the simplex + * should be the same as the codimension of the manifold + * (the codomain dimension of the function). + * @param[in] triangulation The ambient triangulation. The dimension of + * the triangulation should be the same as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + template <class Simplex_handle, class Triangulation> + Query_result<Simplex_handle> intersects(const Simplex_handle& simplex, const Triangulation& triangulation) const { + Eigen::VectorXd lambda = compute_lambda(simplex, triangulation); + return intersection_result(lambda, simplex, triangulation); + } + + /** \brief Intersection query with the boundary of the manifold. + * + * \details The returned structure Query_result contains the boolean value + * that is true only if the intersection point of the query simplex and + * the boundary of the manifold exists, the intersection point + * and the face of the query simplex that contains + * the intersection point. + * + * \tparam Simplex_handle The class of the query simplex. + * Needs to be a model of the concept SimplexInCoxeterTriangulation. + * \tparam Triangulation The class of the triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * + * @param[in] simplex The query simplex. The dimension of the simplex + * should be the same as the codimension of the boundary of the manifold + * (the codomain dimension of the function + 1). + * @param[in] triangulation The ambient triangulation. The dimension of + * the triangulation should be the same as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + template <class Simplex_handle, class Triangulation> + Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex, + const Triangulation& triangulation) const { + //std::cout << "intersects_boundary" << std::endl; + Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation); + return intersection_result(lambda, simplex, triangulation); + } + + /** \brief Returns true if the input point lies inside the piecewise-linear + * domain induced by the given ambient triangulation that defines the relative + * interior of the piecewise-linear approximation of the manifold. + * + * @param p The input point. Needs to have the same dimension as the ambient + * dimension of the manifold (the domain dimension of the function). + * @param triangulation The ambient triangulation. Needs to have the same + * dimension as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + template <class Triangulation> + bool lies_in_domain(const Eigen::VectorXd& p, const Triangulation& triangulation) const { + Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p); + return pl_p(0) < 0; + } + + /** \brief Returns the function that defines the interior of the manifold */ + const Function_& function() const { return fun_; } + + /** \brief Constructs an intersection oracle for an implicit manifold potentially + * with boundary from given function and domain. + * + * @param function The input function that represents the implicit manifold + * before the restriction with the domain. + * @param domain_function The input domain function that can be used to define an implicit + * manifold with boundary. + */ + Implicit_manifold_intersection_oracle(const Function_& function, const Domain_function_& domain_function) + : fun_(function), domain_fun_(domain_function) {} + + /** \brief Constructs an intersection oracle for an implicit manifold + * without boundary from a given function. + * + * \details To use this constructor, the template Domain_function_ needs to be left + * at its default value (Gudhi::coxeter_triangulation::Constant_function). + * + * @param function The input function that represents the implicit manifold + * without boundary. + */ + Implicit_manifold_intersection_oracle(const Function_& function) + : fun_(function), domain_fun_(function.amb_d(), 1, Eigen::VectorXd::Constant(1, -1)) {} + + private: + Function_ fun_; + Domain_function_ domain_fun_; +}; + +/** \brief Static constructor of an intersection oracle from a function with a domain. + * + * @param function The input function that represents the implicit manifold + * before the restriction with the domain. + * @param domain_function The input domain function that can be used to define an implicit + * manifold with boundary. + * + * \ingroup coxeter_triangulation + */ +template <class Function_, class Domain_function_> +Implicit_manifold_intersection_oracle<Function_, Domain_function_> make_oracle( + const Function_& function, const Domain_function_& domain_function) { + return Implicit_manifold_intersection_oracle<Function_, Domain_function_>(function, domain_function); +} + +/** \brief Static constructor of an intersection oracle from a function without a domain. + * + * @param function The input function that represents the implicit manifold + * without boundary. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Implicit_manifold_intersection_oracle<Function_> make_oracle(const Function_& function) { + return Implicit_manifold_intersection_oracle<Function_>(function); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h new file mode 100644 index 00000000..d61bbed7 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h @@ -0,0 +1,270 @@ +/* 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 MANIFOLD_TRACING_H_ +#define MANIFOLD_TRACING_H_ + +#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES +#include <gudhi/Coxeter_triangulation/Query_result.h> + +#include <boost/functional/hash.hpp> + +#include <Eigen/Dense> + +#include <queue> +#include <unordered_map> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \ingroup coxeter_triangulation + */ + +/** \class Manifold_tracing + * \brief A class that assembles methods for manifold tracing algorithm. + * + * \tparam Triangulation_ The type of the ambient triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + */ +template <class Triangulation_> +class Manifold_tracing { + public: + using Simplex_handle = typename Triangulation_::Simplex_handle; + + struct Simplex_hash { + typedef Simplex_handle argument_type; + typedef std::size_t result_type; + result_type operator()(const argument_type& s) const noexcept { + return boost::hash<typename Simplex_handle::Vertex>()(s.vertex()); + } + }; + + public: + /** \brief Type of the output simplex map with keys of type Triangulation_::Simplex_handle + * and values of type Eigen::VectorXd. + * This type should be used for the output in the method manifold_tracing_algorithm. + */ + typedef std::unordered_map<Simplex_handle, Eigen::VectorXd, Simplex_hash> Out_simplex_map; + + /** + * \brief Computes the set of k-simplices that intersect + * a boundaryless implicit manifold given by an intersection oracle, where k + * is the codimension of the manifold. + * The computation is based on the seed propagation --- it starts at the + * given seed points and then propagates along the manifold. + * + * \tparam Point_range Range of points of type Eigen::VectorXd. + * \tparam Intersection_oracle Intersection oracle that represents the manifold. + * Needs to be a model of the concept IntersectionOracle. + * + * \param[in] seed_points The range of points on the manifold from which + * the computation begins. + * \param[in] triangulation The ambient triangulation. + * \param[in] oracle The intersection oracle for the manifold. + * The ambient dimension needs to match the dimension of the + * triangulation. + * \param[out] out_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the input manifold and the mapped values + * are the intersection points. + */ + template <class Point_range, class Intersection_oracle> + void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation, + const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) { + std::size_t cod_d = oracle.cod_d(); + std::queue<Simplex_handle> queue; + + for (const auto& p : seed_points) { + Simplex_handle full_simplex = triangulation.locate_point(p); + for (Simplex_handle face : full_simplex.face_range(cod_d)) { + Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation); + if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) { +#ifdef DEBUG_TRACES + mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false)); +#endif + queue.emplace(face); + break; + } + } + } + + while (!queue.empty()) { + Simplex_handle s = queue.front(); + queue.pop(); + for (auto cof : s.coface_range(cod_d + 1)) { + for (auto face : cof.face_range(cod_d)) { + Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation); + if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face); + } + } + } + } + + /** + * \brief Computes the set of k-simplices that intersect + * the dimensional manifold given by an intersection oracle, where k + * is the codimension of the manifold. + * The computation is based on the seed propagation --- it starts at the + * given seed points and then propagates along the manifold. + * + * \tparam Point_range Range of points of type Eigen::VectorXd. + * \tparam Intersection_oracle Intersection oracle that represents the manifold. + * Needs to be a model of the concept IntersectionOracle. + * + * \param[in] seed_points The range of points on the manifold from which + * the computation begins. + * \param[in] triangulation The ambient triangulation. + * \param[in] oracle The intersection oracle for the manifold. + * The ambient dimension needs to match the dimension of the + * triangulation. + * \param[out] interior_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the relative interior of the input manifold + * and the mapped values are the intersection points. + * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the boundary of the input manifold + * and the mapped values are the intersection points. + */ + template <class Point_range, class Intersection_oracle> + void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation, + const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map, + Out_simplex_map& boundary_simplex_map) { + std::size_t cod_d = oracle.cod_d(); + std::queue<Simplex_handle> queue; + + for (const auto& p : seed_points) { + Simplex_handle full_simplex = triangulation.locate_point(p); + for (Simplex_handle face : full_simplex.face_range(cod_d)) { + auto qr = oracle.intersects(face, triangulation); +#ifdef DEBUG_TRACES + mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false)); +#endif + if (qr.success) { + if (oracle.lies_in_domain(qr.intersection, triangulation)) { + if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face); + } else { + for (Simplex_handle cof : face.coface_range(cod_d + 1)) { + auto qrb = oracle.intersects_boundary(cof, triangulation); +#ifdef DEBUG_TRACES + mt_seed_inserted_list.push_back(MT_inserted_info(qrb, cof, true)); +#endif + if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection); + } + } + // break; + } + } + } + + while (!queue.empty()) { + Simplex_handle s = queue.front(); + queue.pop(); + for (auto cof : s.coface_range(cod_d + 1)) { + for (auto face : cof.face_range(cod_d)) { + auto qr = oracle.intersects(face, triangulation); +#ifdef DEBUG_TRACES + mt_inserted_list.push_back(MT_inserted_info(qr, face, false)); +#endif + if (qr.success) { + if (oracle.lies_in_domain(qr.intersection, triangulation)) { + if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face); + } else { + auto qrb = oracle.intersects_boundary(cof, triangulation); +#ifdef DEBUG_TRACES + mt_inserted_list.push_back(MT_inserted_info(qrb, cof, true)); +#endif + if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection); + } + } + } + } + } + } + + /** \brief Empty constructor */ + Manifold_tracing() {} +}; + +/** + * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm + * that computes the set of k-simplices that intersect + * a boundaryless implicit manifold given by an intersection oracle, where k + * is the codimension of the manifold. + * The computation is based on the seed propagation --- it starts at the + * given seed points and then propagates along the manifold. + * + * \tparam Point_range Range of points of type Eigen::VectorXd. + * \tparam Triangulation_ The type of the ambient triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * \tparam Intersection_oracle Intersection oracle that represents the manifold. + * Needs to be a model of the concept IntersectionOracle. + * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map. + * + * \param[in] seed_points The range of points on the manifold from which + * the computation begins. + * \param[in] triangulation The ambient triangulation. + * \param[in] oracle The intersection oracle for the manifold. + * The ambient dimension needs to match the dimension of the + * triangulation. + * \param[out] out_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the input manifold and the mapped values + * are the intersection points. + * + * \ingroup coxeter_triangulation + */ +template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map> +void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation, + const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) { + Manifold_tracing<Triangulation> mt; + mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map); +} + +/** + * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm + * the dimensional manifold given by an intersection oracle, where k + * is the codimension of the manifold. + * The computation is based on the seed propagation --- it starts at the + * given seed points and then propagates along the manifold. + * + * \tparam Point_range Range of points of type Eigen::VectorXd. + * \tparam Triangulation_ The type of the ambient triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * \tparam Intersection_oracle Intersection oracle that represents the manifold. + * Needs to be a model of the concept IntersectionOracle. + * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map. + * + * \param[in] seed_points The range of points on the manifold from which + * the computation begins. + * \param[in] triangulation The ambient triangulation. + * \param[in] oracle The intersection oracle for the manifold. + * The ambient dimension needs to match the dimension of the + * triangulation. + * \param[out] interior_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the relative interior of the input manifold + * and the mapped values are the intersection points. + * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in + * the input triangulation that intersect the boundary of the input manifold + * and the mapped values are the intersection points. + * + * \ingroup coxeter_triangulation + */ +template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map> +void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation, + const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map, + Out_simplex_map& boundary_simplex_map) { + Manifold_tracing<Triangulation> mt; + mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h new file mode 100644 index 00000000..76438c91 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h @@ -0,0 +1,216 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_H_ +#define PERMUTAHEDRAL_REPRESENTATION_H_ + +#include <gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h> + +#include <utility> // for std::make_pair + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Permutahedral_representation + * \brief A class that stores the permutahedral representation of a simplex + * in a Coxeter triangulation or a Freudenthal-Kuhn triangulation. + * + * \ingroup coxeter_triangulation + * + * \details The data structure is a record consisting of a range that + * represents the vertex and a range that represents the ordered set + * partition, both of which identify the simplex in the triangulation. + * + * \tparam Vertex_ needs to be a random-access range. + * \tparam Ordered_set_partition_ needs to be a a random-access range that consists of + * random-access ranges. + */ +template <class Vertex_, class Ordered_set_partition_> +class Permutahedral_representation { + typedef Permutahedral_representation<Vertex_, Ordered_set_partition_> Self; + + public: + /** \brief Type of the vertex. */ + typedef Vertex_ Vertex; + + /** \brief Type of the ordered partition. */ + typedef Ordered_set_partition_ OrderedSetPartition; + + /** \brief Permutahedral_representation constructor from a vertex and an ordered set partition. + * + * @param[in] vertex Vertex. + * @param[in] partition Ordered set partition. + * + * \details If the size of vertex is d, the ranges in partition must consist + * of the integers 0,...,d without repetition or collision between the ranges. + */ + Permutahedral_representation(const Vertex& vertex, const OrderedSetPartition& partition) + : vertex_(vertex), partition_(partition) {} + + /** \brief Constructor for an empty permutahedral representation that does not correspond + * to any simplex. + */ + Permutahedral_representation() {} + + /** \brief Dimension of the simplex. */ + std::size_t dimension() const { return partition_.size() - 1; } + + /** \brief Lexicographically-minimal vertex. */ + Vertex& vertex() { return vertex_; } + + /** \brief Lexicographically-minimal vertex. */ + const Vertex& vertex() const { return vertex_; } + + /** \brief Ordered set partition. */ + OrderedSetPartition& partition() { return partition_; } + + /** \brief Identifying vertex. */ + const OrderedSetPartition& partition() const { return partition_; } + + /** \brief Equality operator. + * Returns true if an only if both vertex and the ordered set partition coincide. + */ + bool operator==(const Permutahedral_representation& other) const { + if (dimension() != other.dimension()) return false; + if (vertex_ != other.vertex_) return false; + for (std::size_t k = 0; k < partition_.size(); ++k) + if (partition_[k] != other.partition_[k]) return false; + return true; + } + + /** \brief Inequality operator. + * Returns true if an only if either vertex or the ordered set partition are different. + */ + bool operator!=(const Permutahedral_representation& other) const { return !(*this == other); } + + typedef Gudhi::coxeter_triangulation::Vertex_iterator<Self> Vertex_iterator; + typedef boost::iterator_range<Vertex_iterator> Vertex_range; + /** \brief Returns a range of vertices of the simplex. + * The type of vertices is Vertex. + */ + Vertex_range vertex_range() const { return Vertex_range(Vertex_iterator(*this), Vertex_iterator()); } + + typedef Gudhi::coxeter_triangulation::Face_iterator<Self> Face_iterator; + typedef boost::iterator_range<Face_iterator> Face_range; + /** \brief Returns a range of permutahedral representations of faces of the simplex. + * @param[in] value_dim The dimension of the faces. Must be between 0 and the dimension of the simplex. + */ + Face_range face_range(std::size_t value_dim) const { + return Face_range(Face_iterator(*this, value_dim), Face_iterator()); + } + + /** \brief Returns a range of permutahedral representations of facets of the simplex. + * The dimension of the simplex must be strictly positive. + */ + Face_range facet_range() const { return Face_range(Face_iterator(*this, dimension() - 1), Face_iterator()); } + + typedef Gudhi::coxeter_triangulation::Coface_iterator<Self> Coface_iterator; + typedef boost::iterator_range<Coface_iterator> Coface_range; + /** \brief Returns a range of permutahedral representations of cofaces of the simplex. + * @param[in] value_dim The dimension of the cofaces. Must be between the dimension of the simplex and the ambient + * dimension (the size of the vertex). + */ + Coface_range coface_range(std::size_t value_dim) const { + return Coface_range(Coface_iterator(*this, value_dim), Coface_iterator()); + } + + /** \brief Returns a range of permutahedral representations of cofacets of the simplex. + * The dimension of the simplex must be strictly different from the ambient dimension (the size of the vertex). + */ + Coface_range cofacet_range() const { + return Coface_range(Coface_iterator(*this, dimension() + 1), Coface_iterator()); + } + + /** \brief Returns true, if the simplex is a face of other simplex. + * + * @param[in] other A simplex that is potential a coface of the current simplex. + */ + bool is_face_of(const Permutahedral_representation& other) const { + using Part = typename OrderedSetPartition::value_type; + + if (other.dimension() < dimension()) return false; + if (other.vertex_.size() != vertex_.size()) + std::cerr << "Error: Permutahedral_representation::is_face_of: incompatible ambient dimensions.\n"; + + Vertex v_self = vertex_, v_other = other.vertex_; + auto self_partition_it = partition_.begin(); + auto other_partition_it = other.partition_.begin(); + while (self_partition_it != partition_.end()) { + while (other_partition_it != other.partition_.end() && v_self != v_other) { + const Part& other_part = *other_partition_it++; + if (other_partition_it == other.partition_.end()) return false; + for (const auto& k : other_part) v_other[k]++; + } + if (other_partition_it == other.partition_.end()) return false; + const Part& self_part = *self_partition_it++; + if (self_partition_it == partition_.end()) return true; + for (const auto& k : self_part) v_self[k]++; + } + return true; + } + + private: + Vertex vertex_; + OrderedSetPartition partition_; +}; + +/** \brief Print a permutahedral representation to a stream. + * \ingroup coxeter_triangulation + * + * @param[in] os The output stream. + * @param[in] simplex A simplex represented by its permutahedral representation. + */ +template <class Vertex, class OrderedSetPartition> +std::ostream& operator<<(std::ostream& os, const Permutahedral_representation<Vertex, OrderedSetPartition>& simplex) { + // vertex part + os << "("; + if (simplex.vertex().empty()) { + os << ")"; + return os; + } + auto v_it = simplex.vertex().begin(); + os << *v_it++; + for (; v_it != simplex.vertex().end(); ++v_it) os << ", " << *v_it; + os << ")"; + + // ordered partition part + using Part = typename OrderedSetPartition::value_type; + auto print_part = [&os](const Part& p) { + os << "{"; + if (p.empty()) { + os << "}"; + } + auto p_it = p.begin(); + os << *p_it++; + for (; p_it != p.end(); ++p_it) os << ", " << *p_it; + os << "}"; + }; + os << " ["; + if (simplex.partition().empty()) { + os << "]"; + return os; + } + auto o_it = simplex.partition().begin(); + print_part(*o_it++); + for (; o_it != simplex.partition().end(); ++o_it) { + os << ", "; + print_part(*o_it); + } + os << "]"; + return os; +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif // PERMUTAHEDRAL_REPRESENTATION_H_ diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h new file mode 100644 index 00000000..5f382e31 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h @@ -0,0 +1,83 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_ + +#include <vector> +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +typedef unsigned uint; + +/** \brief Class that allows the user to generate combinations of + * k elements in a set of n elements. + * Based on the algorithm by Mifsud. + */ +class Combination_iterator + : public boost::iterator_facade<Combination_iterator, std::vector<uint> const, boost::forward_traversal_tag> { + typedef std::vector<uint> value_t; + + protected: + friend class boost::iterator_core_access; + + bool equal(Combination_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void increment() { + if (value_[0] == n_ - k_) { + is_end_ = true; + return; + } + uint j = k_ - 1; + if (value_[j] < n_ - 1) { + value_[j]++; + return; + } + for (; j > 0; --j) + if (value_[j - 1] < n_ - k_ + j - 1) { + value_[j - 1]++; + for (uint s = j; s < k_; s++) value_[s] = value_[j - 1] + s - (j - 1); + return; + } + } + + public: + Combination_iterator(const uint& n, const uint& k) : value_(k), is_end_(n == 0), n_(n), k_(k) { + for (uint i = 0; i < k; ++i) value_[i] = i; + } + + // Used for the creating an end iterator + Combination_iterator() : is_end_(true), n_(0), k_(0) {} + + void reinitialize() { + if (n_ > 0) { + is_end_ = false; + for (uint i = 0; i < n_; ++i) value_[i] = i; + } + } + + private: + value_t value_; // the dereference value + bool is_end_; // is true when the current permutation is the final one + + uint n_; + uint k_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h new file mode 100644 index 00000000..594b6fbf --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h @@ -0,0 +1,113 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_ + +#include <vector> +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +typedef unsigned uint; + +/** \brief Class that allows the user to generate combinations of + * k elements in a set of n elements. + * Based on the algorithm by Mifsud. + */ +class Integer_combination_iterator + : public boost::iterator_facade<Integer_combination_iterator, std::vector<uint> const, + boost::forward_traversal_tag> { + using value_t = std::vector<uint>; + + private: + friend class boost::iterator_core_access; + + bool equal(Integer_combination_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void increment() { + uint j1 = 0; + uint s = 0; + while (value_[j1] == 0 && j1 < k_) j1++; + uint j2 = j1 + 1; + while (value_[j2] == bounds_[j2]) { + if (bounds_[j2] != 0) { + s += value_[j1]; + value_[j1] = 0; + j1 = j2; + } + j2++; + } + if (j2 >= k_) { + is_end_ = true; + return; + } + s += value_[j1] - 1; + value_[j1] = 0; + value_[j2]++; + uint i = 0; + while (s >= bounds_[i]) { + value_[i] = bounds_[i]; + s -= bounds_[i]; + i++; + } + value_[i++] = s; + } + + public: + template <class Bound_range> + Integer_combination_iterator(const uint& n, const uint& k, const Bound_range& bounds) + : value_(k + 2), is_end_(n == 0 || k == 0), k_(k) { + bounds_.reserve(k + 2); + uint sum_radices = 0; + for (auto b : bounds) { + bounds_.push_back(b); + sum_radices += b; + } + bounds_.push_back(2); + bounds_.push_back(1); + if (n > sum_radices) { + is_end_ = true; + return; + } + uint i = 0; + uint s = n; + while (s >= bounds_[i]) { + value_[i] = bounds_[i]; + s -= bounds_[i]; + i++; + } + value_[i++] = s; + + while (i < k_) value_[i++] = 0; + value_[k] = 1; + value_[k + 1] = 0; + } + + // Used for the creating an end iterator + Integer_combination_iterator() : is_end_(true), k_(0) {} + + private: + value_t value_; // the dereference value + bool is_end_; // is true when the current integer combination is the final one + + uint k_; + std::vector<uint> bounds_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h new file mode 100644 index 00000000..866079fa --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h @@ -0,0 +1,93 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_ + +#include <vector> +#include <limits> + +#include <gudhi/Permutahedral_representation/Permutation_iterator.h> +#include <gudhi/Permutahedral_representation/Set_partition_iterator.h> + +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +typedef unsigned uint; + +/** \brief Class that represents an ordered set partition of a set {0,...,n-1} in k parts as + * a pair of an unordered set partition given in lexicographic order and + * a permutation of the parts. + */ +struct Ordered_set_partition { + Set_partition_iterator s_it_; + Permutation_iterator p_it_; + + // Ordered_set_partition(const Set_partition_iterator& s_it, const Permutation_iterator& p_it) + // : s_it_(s_it), p_it_(p_it) {} + + const std::vector<uint> operator[](const uint& i) const { return (*s_it_)[(*p_it_)[i]]; } + + std::size_t size() const { return s_it_->size(); } +}; + +/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts. + * + */ +class Ordered_set_partition_iterator + : public boost::iterator_facade<Ordered_set_partition_iterator, Ordered_set_partition const, + boost::forward_traversal_tag> { + using value_t = Ordered_set_partition; + + private: + friend class boost::iterator_core_access; + + bool equal(Ordered_set_partition_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void increment() { + if (++value_.p_it_ == p_end_) { + if (++value_.s_it_ == s_end_) { + is_end_ = true; + return; + } else + value_.p_it_.reinitialize(); + } + } + + public: + Ordered_set_partition_iterator(const uint& n, const uint& k) + : value_({Set_partition_iterator(n, k), Permutation_iterator(k)}), is_end_(n == 0) {} + + // Used for the creating an end iterator + Ordered_set_partition_iterator() : is_end_(true) {} + + void reinitialize() { + is_end_ = false; + value_.p_it_.reinitialize(); + value_.s_it_.reinitialize(); + } + + private: + Set_partition_iterator s_end_; // Set partition iterator and the corresponding end iterator + Permutation_iterator p_end_; // Permutation iterator and the corresponding end iterator + value_t value_; // the dereference value + bool is_end_; // is true when the current permutation is the final one +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h new file mode 100644 index 00000000..1a63d2f7 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h @@ -0,0 +1,256 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_ +#define PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_ + +#include <gudhi/Permutahedral_representation/Size_range.h> +#include <gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h> +#include <gudhi/Permutahedral_representation/Integer_combination_iterator.h> +#include <gudhi/Permutahedral_representation/Combination_iterator.h> +#include <gudhi/Permutahedral_representation/face_from_indices.h> +#include <boost/iterator/iterator_facade.hpp> + +#include <vector> +#include <iostream> +#include <algorithm> // for std::find + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \addtogroup coxeter_triangulation + * Iterator types for Permutahedral_representation + * @{ + */ + +/** \brief Iterator over the vertices of a simplex + * represented by its permutahedral representation. + * + * Forward iterator, 'value_type' is Permutahedral_representation::Vertex.*/ +template <class Permutahedral_representation> +class Vertex_iterator + : public boost::iterator_facade<Vertex_iterator<Permutahedral_representation>, + typename Permutahedral_representation::Vertex const, boost::forward_traversal_tag> { + private: + friend class boost::iterator_core_access; + + using Vertex = typename Permutahedral_representation::Vertex; + using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition; + + using value_t = Vertex; + + bool equal(Vertex_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void update_value() { + std::size_t d = value_.size(); + for (auto i : *o_it_) + if (i != d) + value_[i]++; + else + for (std::size_t j = 0; j < d; j++) value_[j]--; + } + + void increment() { + if (is_end_) return; + update_value(); + if (++o_it_ == o_end_) is_end_ = true; + } + + public: + Vertex_iterator(const Permutahedral_representation& simplex) + : o_it_(simplex.partition().begin()), + o_end_(simplex.partition().end()), + value_(simplex.vertex()), + is_end_(o_it_ == o_end_) {} + + Vertex_iterator() : is_end_(true) {} + + private: + typename Ordered_partition::const_iterator o_it_, o_end_; + value_t value_; + bool is_end_; + +}; // Vertex_iterator + +/*---------------------------------------------------------------------------*/ +/** \brief Iterator over the k-faces of a simplex + * given by its permutahedral representation. + * + * Forward iterator, value_type is Permutahedral_representation. */ +template <class Permutahedral_representation> +class Face_iterator : public boost::iterator_facade<Face_iterator<Permutahedral_representation>, + Permutahedral_representation const, boost::forward_traversal_tag> { + using value_t = Permutahedral_representation; + + private: + friend class boost::iterator_core_access; + + using Vertex = typename Permutahedral_representation::Vertex; + using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition; + + bool equal(Face_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void increment() { + if (++c_it_ == c_end_) { + is_end_ = true; + return; + } + update_value(); + } + + void update_value() { + // Combination *c_it_ is supposed to be sorted in increasing order + value_ = face_from_indices<Permutahedral_representation>(simplex_, *c_it_); + } + + public: + Face_iterator(const Permutahedral_representation& simplex, const uint& k) + : simplex_(simplex), + k_(k), + l_(simplex.dimension()), + c_it_(l_ + 1, k_ + 1), + is_end_(k_ > l_), + value_({Vertex(simplex.vertex().size()), Ordered_partition(k + 1)}) { + update_value(); + } + + // Used for the creating an end iterator + Face_iterator() : is_end_(true) {} + + private: + Permutahedral_representation simplex_; // Input simplex + uint k_; + uint l_; // Dimension of the input simplex + Combination_iterator c_it_, c_end_; // indicates the vertices in the current face + + bool is_end_; // is true when the current permutation is the final one + value_t value_; // the dereference value + +}; // Face_iterator + +/*---------------------------------------------------------------------------*/ +/** \brief Iterator over the k-cofaces of a simplex + * given by its permutahedral representation. + * + * Forward iterator, value_type is Permutahedral_representation. */ +template <class Permutahedral_representation> +class Coface_iterator + : public boost::iterator_facade<Coface_iterator<Permutahedral_representation>, Permutahedral_representation const, + boost::forward_traversal_tag> { + using value_t = Permutahedral_representation; + + private: + friend class boost::iterator_core_access; + + using Vertex = typename Permutahedral_representation::Vertex; + using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition; + + bool equal(Coface_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void increment() { + uint i = 0; + for (; i < k_ + 1; i++) { + if (++(o_its_[i]) != o_end_) break; + } + if (i == k_ + 1) { + if (++i_it_ == i_end_) { + is_end_ = true; + return; + } + o_its_.clear(); + for (uint j = 0; j < k_ + 1; j++) + o_its_.emplace_back(simplex_.partition()[j].size(), (*i_it_)[j] + 1); + } else + for (uint j = 0; j < i; j++) o_its_[j].reinitialize(); + update_value(); + } + + void update_value() { + value_.vertex() = simplex_.vertex(); + for (auto& p : value_.partition()) p.clear(); + uint u_ = 0; // the part in o_its_[k_] that contains t_ + for (; u_ <= (*i_it_)[k_]; u_++) { + auto range = (*o_its_[k_])[u_]; + if (std::find(range.begin(), range.end(), t_) != range.end()) break; + } + uint i = 0; + for (uint j = u_ + 1; j <= (*i_it_)[k_]; j++, i++) + for (uint b : (*o_its_[k_])[j]) { + uint c = simplex_.partition()[k_][b]; + value_.partition()[i].push_back(c); + value_.vertex()[c]--; + } + for (uint h = 0; h < k_; h++) + for (uint j = 0; j <= (*i_it_)[h]; j++, i++) { + for (uint b : (*o_its_[h])[j]) value_.partition()[i].push_back(simplex_.partition()[h][b]); + } + for (uint j = 0; j <= u_; j++, i++) + for (uint b : (*o_its_[k_])[j]) value_.partition()[i].push_back(simplex_.partition()[k_][b]); + // sort the values in each part (probably not needed) + for (auto& part : value_.partition()) std::sort(part.begin(), part.end()); + } + + public: + Coface_iterator(const Permutahedral_representation& simplex, const uint& l) + : simplex_(simplex), + d_(simplex.vertex().size()), + l_(l), + k_(simplex.dimension()), + i_it_(l_ - k_, k_ + 1, Size_range<Ordered_partition>(simplex.partition())), + is_end_(k_ > l_), + value_({Vertex(d_), Ordered_partition(l_ + 1)}) { + uint j = 0; + for (; j < simplex_.partition()[k_].size(); j++) + if (simplex_.partition()[k_][j] == d_) { + t_ = j; + break; + } + if (j == simplex_.partition()[k_].size()) { + std::cerr << "Coface iterator: the argument simplex is not a permutahedral representation\n"; + is_end_ = true; + return; + } + for (uint i = 0; i < k_ + 1; i++) + o_its_.emplace_back(simplex_.partition()[i].size(), (*i_it_)[i] + 1); + update_value(); + } + + // Used for the creating an end iterator + Coface_iterator() : is_end_(true) {} + + private: + Permutahedral_representation simplex_; // Input simplex + uint d_; // Ambient dimension + uint l_; // Dimension of the coface + uint k_; // Dimension of the input simplex + uint t_; // The position of d in simplex_.partition()[k_] + Integer_combination_iterator i_it_, i_end_; // indicates in how many parts each simplex_[i] is subdivided + std::vector<Ordered_set_partition_iterator> o_its_; // indicates subdivision for each simplex_[i] + Ordered_set_partition_iterator o_end_; // one end for all o_its_ + + bool is_end_; // is true when the current permutation is the final one + value_t value_; // the dereference value + +}; // Coface_iterator + +/** @} */ + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h new file mode 100644 index 00000000..0f91d41c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h @@ -0,0 +1,120 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_ + +#include <cstdlib> // for std::size_t +#include <vector> + +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +typedef unsigned uint; + +/** \brief Class that allows the user to generate permutations. + * Based on the optimization of the Heap's algorithm by Sedgewick. + */ +class Permutation_iterator + : public boost::iterator_facade<Permutation_iterator, std::vector<uint> const, boost::forward_traversal_tag> { + using value_t = std::vector<uint>; + + private: + friend class boost::iterator_core_access; + + bool equal(Permutation_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void swap_two_indices(std::size_t i, std::size_t j) { + uint t = value_[i]; + value_[i] = value_[j]; + value_[j] = t; + } + + void elementary_increment() { + uint j = 0; + while (d_[j] == j + 1) { + d_[j] = 0; + ++j; + } + if (j == n_ - 1) { + is_end_ = true; + return; + } + uint k = j + 1; + uint x = (k % 2 ? d_[j] : 0); + swap_two_indices(k, x); + ++d_[j]; + } + + void elementary_increment_optim_3() { + if (ct_ != 0) { + --ct_; + swap_two_indices(1 + (ct_ % 2), 0); + } else { + ct_ = 5; + uint j = 2; + while (d_[j] == j + 1) { + d_[j] = 0; + ++j; + } + if (j == n_ - 1) { + is_end_ = true; + return; + } + uint k = j + 1; + uint x = (k % 2 ? d_[j] : 0); + swap_two_indices(k, x); + ++d_[j]; + } + } + + void increment() { + if (optim_3_) + elementary_increment_optim_3(); + else + elementary_increment(); + } + + public: + Permutation_iterator(const uint& n) : value_(n), is_end_(n == 0), optim_3_(n >= 3), n_(n), d_(n), ct_(5) { + for (uint i = 0; i < n; ++i) { + value_[i] = i; + d_[i] = 0; + } + if (n > 0) d_[n - 1] = -1; + } + + // Used for the creating an end iterator + Permutation_iterator() : is_end_(true), n_(0) {} + + void reinitialize() { + if (n_ > 0) is_end_ = false; + } + + private: + value_t value_; // the dereference value + bool is_end_; // is true when the current permutation is the final one + bool optim_3_; // true if n>=3. for n >= 3, the algorithm is optimized + + uint n_; + std::vector<uint> d_; // mix radix digits with radix [2 3 4 ... n-1 (sentinel=-1)] + uint ct_; // counter with values in {0,...,5} used in the n>=3 optimization. +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h new file mode 100644 index 00000000..94ac10c2 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h @@ -0,0 +1,111 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_ + +#include <vector> +#include <limits> +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +typedef unsigned uint; + +/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts. + * + */ +class Set_partition_iterator + : public boost::iterator_facade<Set_partition_iterator, std::vector<std::vector<uint>> const, + boost::forward_traversal_tag> { + using value_t = std::vector<std::vector<uint>>; + + private: + friend class boost::iterator_core_access; + + bool equal(Set_partition_iterator const& other) const { return (is_end_ && other.is_end_); } + + value_t const& dereference() const { return value_; } + + void update_value() { + for (uint i = 0; i < k_; i++) value_[i].clear(); + for (uint i = 0; i < n_; i++) value_[rgs_[i]].push_back(i); + } + + void increment() { + if (k_ <= 1) { + is_end_ = true; + return; + } + uint i = n_ - 1; + while (rgs_[i] + 1 > max_[i] || rgs_[i] + 1 >= k_) i--; + if (i == 0) { + is_end_ = true; + return; + } + rgs_[i]++; + uint mm = max_[i]; + mm += (rgs_[i] >= mm); + max_[i + 1] = mm; + while (++i < n_) { + rgs_[i] = 0; + max_[i + 1] = mm; + } + uint p = k_; + if (mm < p) do { + max_[i] = p; + --i; + --p; + rgs_[i] = p; + } while (max_[i] < p); + update_value(); + } + + public: + Set_partition_iterator(const uint& n, const uint& k) + : value_(k), rgs_(n, 0), max_(n + 1), is_end_(n == 0), n_(n), k_(k) { + max_[0] = std::numeric_limits<uint>::max(); + for (uint i = 0; i <= n - k; ++i) value_[0].push_back(i); + for (uint i = n - k + 1, j = 1; i < n; ++i, ++j) { + rgs_[i] = j; + value_[j].push_back(i); + } + for (uint i = 1; i <= n; i++) max_[i] = rgs_[i - 1] + 1; + update_value(); + } + + // Used for creating an end iterator + Set_partition_iterator() : is_end_(true), n_(0), k_(0) {} + + void reinitialize() { + if (n_ > 0) is_end_ = false; + for (uint i = 0; i <= n_ - k_; ++i) rgs_[i] = 0; + for (uint i = n_ - k_ + 1, j = 1; i < n_; ++i, ++j) rgs_[i] = j; + for (uint i = 1; i <= n_; i++) max_[i] = rgs_[i - 1] + 1; + update_value(); + } + + private: + value_t value_; // the dereference value + std::vector<uint> rgs_; // restricted growth string + std::vector<uint> max_; // max_[i] = max(rgs_[0],...,rgs[i-1]) + 1 + bool is_end_; // is true when the current permutation is the final one + + uint n_; + uint k_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h new file mode 100644 index 00000000..905d68d5 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h @@ -0,0 +1,54 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_ +#define PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_ + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \class Simplex_comparator + * \brief A comparator class for Permutahedral_representation. + * The comparison is in lexicographic order first on + * vertices and then on ordered partitions with sorted parts. + * The lexicographic order forces that any face is larger than + * a coface. + * + * \tparam Permutahdral_representation_ Needs to be + * Permutahedral_representation<Vertex_, Ordered_set_partition_> + * + * \ingroup coxeter_triangulation + */ +template <class Permutahedral_representation_> +struct Simplex_comparator { + /** \brief Comparison between two permutahedral representations. + * Both permutahedral representations need to be valid and + * the vertices of both permutahedral representations need to be of the same size. + */ + bool operator()(const Permutahedral_representation_& lhs, const Permutahedral_representation_& rhs) const { + if (lhs.vertex() < rhs.vertex()) return true; + if (lhs.vertex() > rhs.vertex()) return false; + + if (lhs.partition().size() > rhs.partition().size()) return true; + if (lhs.partition().size() < rhs.partition().size()) return false; + + if (lhs.partition() < rhs.partition()) return true; + if (lhs.partition() > rhs.partition()) return false; + + return false; + } +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h new file mode 100644 index 00000000..6b137744 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h @@ -0,0 +1,73 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_ +#define PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_ + +#include <cstdlib> // for std::size_t + +#include <boost/range/iterator_range.hpp> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \brief Auxiliary iterator class for sizes of parts in an ordered set partition. + */ +template <class T_it> +class Size_iterator + : public boost::iterator_facade<Size_iterator<T_it>, std::size_t const, boost::forward_traversal_tag> { + friend class boost::iterator_core_access; + + private: + bool equal(Size_iterator const& other) const { return (is_end_ && other.is_end_); } + + std::size_t const& dereference() const { return value_; } + + void increment() { + if (++t_it_ == t_end_) { + is_end_ = true; + return; + } + value_ = t_it_->size() - 1; + } + + public: + Size_iterator(const T_it& t_begin, const T_it& t_end) : t_it_(t_begin), t_end_(t_end), is_end_(t_begin == t_end) { + if (!is_end_) value_ = t_it_->size() - 1; + } + + private: + T_it t_it_, t_end_; + bool is_end_; + std::size_t value_; +}; + +template <class T> +class Size_range { + const T& t_; + + public: + typedef Size_iterator<typename T::const_iterator> iterator; + + Size_range(const T& t) : t_(t) {} + + std::size_t operator[](std::size_t i) const { return t_[i].size() - 1; } + + iterator begin() const { return iterator(t_.begin(), t_.end()); } + + iterator end() const { return iterator(t_.end(), t_.end()); } +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h new file mode 100644 index 00000000..47120689 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h @@ -0,0 +1,66 @@ +/* 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 PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_ +#define PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_ + +#include <cstdlib> // for std::size_t +#include <algorithm> + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \brief Computes the permutahedral representation of a face of a given simplex + * and a range of the vertex indices that compose the face. + * + * \tparam Permutahedral_representation has to be Permutahedral_representation + * \tparam Index_range is a range of unsigned integers taking values in 0,...,k, + * where k is the dimension of the simplex simplex. + * + * @param[in] simplex Input simplex. + * @param[in] indices Input range of indices. + */ +template <class Permutahedral_representation, class Index_range> +Permutahedral_representation face_from_indices(const Permutahedral_representation& simplex, + const Index_range& indices) { + using range_index = typename Index_range::value_type; + using Ordered_set_partition = typename Permutahedral_representation::OrderedSetPartition; + using Part = typename Ordered_set_partition::value_type; + using part_index = typename Part::value_type; + Permutahedral_representation value; + std::size_t d = simplex.vertex().size(); + value.vertex() = simplex.vertex(); + std::size_t k = indices.size() - 1; + value.partition().resize(k + 1); + std::size_t l = simplex.partition().size() - 1; + for (std::size_t h = 1; h < k + 1; h++) + for (range_index i = indices[h - 1]; i < indices[h]; i++) + for (part_index j : simplex.partition()[i]) value.partition()[h - 1].push_back(j); + for (range_index i = indices[k]; i < l + 1; i++) + for (part_index j : simplex.partition()[i]) value.partition()[k].push_back(j); + for (range_index i = 0; i < indices[0]; i++) + for (part_index j : simplex.partition()[i]) { + if (j != d) + value.vertex()[j]++; + else + for (std::size_t l = 0; l < d; l++) value.vertex()[l]--; + value.partition()[k].push_back(j); + } + // sort the values in each part (probably not needed) + for (auto& part : value.partition()) std::sort(part.begin(), part.end()); + return value; +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif |