From ec9953f0dbe0f69074f25cc95442ea0012db7d98 Mon Sep 17 00:00:00 2001 From: Siargey Kachanovich Date: Thu, 17 Oct 2019 21:40:20 +0200 Subject: Added the files from the coxeter branch --- .../include/gudhi/Cell_complex.h | 376 +++++++++++++ .../gudhi/Cell_complex/Hasse_diagram_cell.h | 310 +++++++++++ .../include/gudhi/Coxeter_triangulation.h | 94 ++++ .../include/gudhi/Freudenthal_triangulation.h | 239 +++++++++ .../include/gudhi/Functions/Cartesian_product.h | 171 ++++++ .../include/gudhi/Functions/Constant_function.h | 71 +++ .../include/gudhi/Functions/Embed_in_Rd.h | 106 ++++ .../include/gudhi/Functions/Function.h | 54 ++ .../include/gudhi/Functions/Function_Sm_in_Rd.h | 131 +++++ .../gudhi/Functions/Function_affine_plane_in_Rd.h | 101 ++++ .../include/gudhi/Functions/Function_chair_in_R3.h | 87 +++ .../include/gudhi/Functions/Function_iron_in_R3.h | 73 +++ .../Function_lemniscate_revolution_in_R3.h | 90 ++++ .../gudhi/Functions/Function_moment_curve_in_Rd.h | 89 ++++ .../include/gudhi/Functions/Function_torus_in_R3.h | 77 +++ .../Functions/Function_whitney_umbrella_in_R3.h | 82 +++ .../gudhi/Functions/Linear_transformation.h | 96 ++++ .../include/gudhi/Functions/Negation.h | 92 ++++ .../include/gudhi/Functions/PL_approximation.h | 125 +++++ .../include/gudhi/Functions/Translate.h | 96 ++++ .../gudhi/Functions/random_orthogonal_matrix.h | 70 +++ .../include/gudhi/IO/Mesh_medit.h | 57 ++ .../gudhi/IO/build_mesh_from_cell_complex.h | 174 ++++++ .../include/gudhi/IO/output_debug_traces_to_html.h | 589 +++++++++++++++++++++ .../include/gudhi/IO/output_meshes_to_medit.h | 181 +++++++ .../gudhi/Implicit_manifold_intersection_oracle.h | 279 ++++++++++ .../include/gudhi/Manifold_tracing.h | 307 +++++++++++ .../include/gudhi/Permutahedral_representation.h | 257 +++++++++ .../Combination_iterator.h | 101 ++++ .../Integer_combination_iterator.h | 128 +++++ .../Ordered_set_partition_iterator.h | 108 ++++ .../Permutahedral_representation_iterators.h | 288 ++++++++++ .../Permutation_iterator.h | 137 +++++ .../Set_partition_iterator.h | 137 +++++ .../Simplex_comparator.h | 64 +++ .../Permutahedral_representation/Size_range.h | 89 ++++ .../face_from_indices.h | 71 +++ .../include/gudhi/Query_result.h | 42 ++ 38 files changed, 5639 insertions(+) create mode 100644 src/Coxeter_triangulation/include/gudhi/Cell_complex.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Negation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/Translate.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h create mode 100644 src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h create mode 100644 src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h create mode 100644 src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h create mode 100644 src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h create mode 100644 src/Coxeter_triangulation/include/gudhi/Query_result.h (limited to 'src/Coxeter_triangulation/include/gudhi') diff --git a/src/Coxeter_triangulation/include/gudhi/Cell_complex.h b/src/Coxeter_triangulation/include/gudhi/Cell_complex.h new file mode 100644 index 00000000..51109bcb --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Cell_complex.h @@ -0,0 +1,376 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include +#include // for Hasse_diagram_persistence.h + +#include // for Hasse_cell + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \ingroup coxeter_triangulation + */ + +/** \class Cell_complex + * \brief A class that constructs the cell complex from the output provided by the class + * Gudhi::Manifold_tracing. + * + * \tparam Out_simplex_map_ The type of a map from a simplex type that is a + * model of SimplexInCoxeterTriangulation to Eigen::VectorXd. + * + * \ingroup coxeter_triangulation + */ +template +class Cell_complex { +public: + + /** \brief Type of a simplex in the ambient triangulation. + * Is a model of the concept SimplexInCoxeterTriangulation. + */ + typedef typename Out_simplex_map_::key_type Simplex_handle; + /** \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. + */ + typedef Gudhi::Hasse_diagram::Hasse_diagram_cell Hasse_cell; + /** \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. + */ + typedef std::map > Simplex_cell_map; + /** \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. + */ + typedef std::vector Simplex_cell_maps; + + /** \brief Type of a map from cells in the cell complex to the permutahedral representations + * of the corresponding simplices in the ambient triangulation. + */ + typedef std::map Cell_simplex_map; + + /** \brief Type of a map from vertex cells in the cell complex to the permutahedral representations + * of their Cartesian coordinates. + */ + typedef std::map Cell_point_map; + +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 GUDHI_COX_OUTPUT_TO_HTML + 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(CC_detail_info(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(std::make_pair(simplex, new_cell)); + cell_simplex_map_.emplace(std::make_pair(new_cell, simplex)); +#ifdef GUDHI_COX_OUTPUT_TO_HTML + cc_detail_list.back().status_ = CC_detail_info::Result_type::inserted; +#endif + return new_cell; + } +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(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(std::make_pair(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(std::make_pair(cell, 1)); + } + } + } + } + + void construct_complex_(const Out_simplex_map_& out_simplex_map) { +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(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 GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(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(std::make_pair(new_cell, point)); + } +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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 GUDHI_COX_OUTPUT_TO_HTML + 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 Conxtructor 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_cells_; +}; + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h b/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h new file mode 100644 index 00000000..bd730368 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h @@ -0,0 +1,310 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Pawel Dlotko + * + * Copyright (C) 2017 Swansea University UK + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include +#include +#include +#include +#include +#include + + + + +#ifndef HASSE_DIAGRAM_CELL_H +#define HASSE_DIAGRAM_CELL_H + + +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 + * This is a data structure to store a cell in a general Hasse diagram data structure. It store 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). + * + * Please refer to \ref Hasse_diagram for examples. + * + * 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 +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 >; + + /** + * 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){} + + /** + * Construcor 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 retrive 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 + friend std::vector 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 ( 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 ( 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& c ) + { + //cout << "position : " << c.position << ", dimension : " << c.dimension << ", filtration: " << c.filtration << ", size of boudary : " << c.boundary.size() << "\n"; + out << c.position << " " << c.dimension << " " << c.filtration << std::endl; + for ( 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; + size_t size_of_boundary = this->boundary.size(); + result.reserve( size_of_boundary ); + for ( 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; + size_t size_of_boundary = this->boundary.size(); + result.reserve( size_of_boundary ); + for ( 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::value ) + { + result += " Additional_information: "; + result += std::to_string(this->additional_info); + } + result += " boundary "; + for ( 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 ( 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.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h new file mode 100644 index 00000000..bc59ac12 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h @@ -0,0 +1,94 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include //iota + +#include +#include +#include + +#include +#include +#include + +#include +#include + +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 , std::vector > > > +class Coxeter_triangulation : public Freudenthal_triangulation { + + typedef Eigen::MatrixXd Matrix; + + Matrix root_matrix(unsigned d) { + Matrix cartan(d,d); + for (unsigned i = 0; i < d; i++) { + cartan(i,i) = 1.0; + } + for (unsigned i = 1; i < d; i++) { + cartan(i-1,i) = -0.5; + cartan(i,i-1) = -0.5; + } + for (unsigned i = 0; i < d; i++) + for (unsigned j = 0; j < d; j++) + if (j+1 < i || j > i+1) + cartan(i,j) = 0; + Eigen::SelfAdjointEigenSolver saes(cartan); + Eigen::VectorXd sqrt_diag(d); + for (unsigned i = 0; i < d; ++i) + sqrt_diag(i) = std::sqrt(saes.eigenvalues()[i]); + + Matrix lower(d,d); + for (unsigned i = 0; i < d; i++) + for (unsigned j = 0; j < d; j++) + if (i < j) + lower(i,j) = 0; + else + lower(i,j) = 1; + 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(dimension, root_matrix(dimension)) {} +}; + +} // 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..be5d275c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h @@ -0,0 +1,239 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include //iota + +#include +#include +#include + +#include + +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 , std::vector > > > +class Freudenthal_triangulation { + + typedef Eigen::MatrixXd Matrix; + typedef Eigen::VectorXd Vector; + typedef Eigen::SparseMatrix SparseMatrix; + typedef Eigen::Triplet Triplet; + +public: + + /** \brief Type of the simplices in the triangulation. */ + typedef Permutahedral_representation_ Simplex_handle; + + /** \brief Type of the vertices in the triangulation. */ + typedef typename Permutahedral_representation_::Vertex Vertex_handle; + + + /** \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) + : dimension_(dimension), + matrix_(Matrix::Identity(dimension, dimension)), + offset_(Vector::Zero(dimension)), + colpivhouseholderqr_(matrix_.colPivHouseholderQr()), + 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) + : dimension_(dimension), + matrix_(matrix), + offset_(Vector::Zero(dimension)), + colpivhouseholderqr_(matrix_.colPivHouseholderQr()), + is_freudenthal(false) { + } + + /** \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. + */ + Freudenthal_triangulation(unsigned dimension, const Matrix& matrix, const Vector& offset) + : dimension_(dimension), + matrix_(matrix), + offset_(offset), + colpivhouseholderqr_(matrix_.colPivHouseholderQr()), + is_freudenthal(false) { + } + + + /** \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. + */ + template + 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(); + assert(d == dimension_); + double error = 1e-9; + Simplex_handle output; + std::vector 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_); + for (size_t i = 0; i < dimension_; ++i) + res_vector(i) = 0; + 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 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..43198b89 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h @@ -0,0 +1,171 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include + +#include +#include + +namespace Gudhi { + +namespace coxeter_triangulation { + +/* Get the domain dimension of the tuple of functions. + */ +template +inline typename std::enable_if::type +get_amb_d (const std::tuple& tuple) { + return 0; +} + +template +inline typename std::enable_if::type +get_amb_d (const std::tuple& tuple) { + return std::get(tuple).amb_d() + get_amb_d(tuple); +} + +/* Get the codomain dimension of the tuple of functions. + */ +template +inline typename std::enable_if::type +get_cod_d (const std::tuple& tuple) { + return 0; +} + +template +inline typename std::enable_if::type +get_cod_d (const std::tuple& tuple) { + return std::get(tuple).cod_d() + get_cod_d(tuple); +} + +/* Get the seed of the tuple of functions. + */ +template +inline typename std::enable_if::type +get_seed (const std::tuple& tuple, Eigen::VectorXd& point, std::size_t i = 0) { +} + +template +inline typename std::enable_if::type +get_seed (const std::tuple& tuple, Eigen::VectorXd& point, std::size_t i = 0) { + const auto& f = std::get(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(tuple, point, i+n); +} + +/* Get the seed of the tuple of functions. + */ +template +inline typename std::enable_if::type +get_value (const std::tuple& tuple, const Eigen::VectorXd& x, Eigen::VectorXd& point, std::size_t i = 0, std::size_t j = 0) { +} + +template +inline typename std::enable_if::type +get_value (const std::tuple& tuple, const Eigen::VectorXd& x, Eigen::VectorXd& point, std::size_t i = 0, std::size_t j = 0) { + const auto& f = std::get(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(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. + * + * \ingroup coxeter_triangulation + */ +template +struct Cartesian_product : public Function { + + /** + * \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 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 +Cartesian_product make_product_function(const Functions&... functions) { + return Cartesian_product(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..ea354fac --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h @@ -0,0 +1,71 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Constant_function : public 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 { + Eigen::VectorXd result = value_; + 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 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) {} + + 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..9476944b --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h @@ -0,0 +1,106 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +#include +#include + + +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. + * + * \ingroup coxeter_triangulation + */ +template +struct Embed_in_Rd : public Function { + + /** + * \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) { + } + 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 +Embed_in_Rd make_embedding(const Function_& function, std::size_t d) { + return Embed_in_Rd(function, d); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function.h new file mode 100644 index 00000000..dbcb0142 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function.h @@ -0,0 +1,54 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2019 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef FUNCTIONS_FUNCTION_H_ +#define FUNCTIONS_FUNCTION_H_ + +#include + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** + * \class Function + * \brief The parent class for all functions implemented in the module. + * Contains virtual methods needed to be a model of the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +struct Function { + + /** \brief Virtual method for the value of the function at a specified point. + * @param[in] p The input point. + */ + virtual Eigen::VectorXd evaluate(const Eigen::VectorXd& p) const {return Eigen::VectorXd(0);} + + /** \brief Virtual method for the domain dimension. */ + virtual std::size_t amb_d() const {return 0;}; + + /** \brief Virtual method for the codomain dimension. */ + virtual std::size_t cod_d() const {return 0;}; + + /** \brief Virtual method for the seed point. */ + virtual Eigen::VectorXd seed() const {return Eigen::VectorXd(0);} + + /** \brief Virtual destructor. */ + virtual ~Function() {} +}; + +} // 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..a7d4a965 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h @@ -0,0 +1,131 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_Sm_in_Rd: public Function { + +/** \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_) {} + + +protected: + 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..cb3af848 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h @@ -0,0 +1,101 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_affine_plane_in_Rd : public Function { + + /** + * \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 corespond 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()), + m_(d_ - k_), + 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 corespond 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()), + m_(d_ - k_), + off_(Eigen::VectorXd::Zero(d_)) { + normal_matrix_.colwise().normalize(); + } + +protected: + Eigen::MatrixXd normal_matrix_; + std::size_t d_, k_, m_; + 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..2a76e631 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h @@ -0,0 +1,87 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_chair_in_R3 : public Function { + + /** + * \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..af323c94 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h @@ -0,0 +1,73 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_iron_in_R3 : public Function { + + /** + * \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) {} + +protected: + 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..cd03a0a5 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h @@ -0,0 +1,90 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_lemniscate_revolution_in_R3 : public Function { + + /** + * \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(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(-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) {} + +protected: + 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..91eb016b --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h @@ -0,0 +1,89 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_moment_curve_in_Rd : public Function { + +/** \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 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) + : m_(1), 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) + : m_(1), k_(d-1), d_(d), r_(r), off_(offset) {} + +protected: + std::size_t m_, 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..3cda2518 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h @@ -0,0 +1,77 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_torus_in_R3 : public Function { + + /** + * \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) {} + +protected: + 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..3d825f8f --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h @@ -0,0 +1,82 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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. + * + * \ingroup coxeter_triangulation + */ +struct Function_whitney_umbrella_in_R3 : public Function { + + /** + * \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) {} + 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..4cb5ca84 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h @@ -0,0 +1,96 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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. + * + * \ingroup coxeter_triangulation + */ +template +struct Linear_transformation : public Function { + + /** + * \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. + */ +template +Linear_transformation make_linear_transformation(const Function_& function, + const Eigen::MatrixXd& matrix) { + return Linear_transformation(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..6b7feff5 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h @@ -0,0 +1,92 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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. + * + * \ingroup coxeter_triangulation + */ +template +struct Negation : public Function { + + /** + * \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) { + } + Function_ fun_; +}; + + +/** + * \brief Static constructor of the negative 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 +Negation negation(const Function_& function) { + return Negation(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..25c8a139 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h @@ -0,0 +1,125 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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. + * + * \ingroup coxeter_triangulation + */ +template +struct PL_approximation : public Function { + + /** + * \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 +PL_approximation +make_pl_approximation(const Function_& function, + const Triangulation_& triangulation) { + return PL_approximation(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..6edac28c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h @@ -0,0 +1,96 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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. + * + * \ingroup coxeter_triangulation + */ +template +struct Translate : public Function { + + /** + * \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) { + } + 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 +Translate translate(const Function_& function, Eigen::VectorXd off) { + return Translate(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..8979669e --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h @@ -0,0 +1,70 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +#include +#include +#include + +#include +#include + +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 Kernel; + typedef typename Kernel::Point_d Point_d; + if (d == 1) + return Eigen::VectorXd::Constant(1, 1.0); + if (d == 2) { + double X = 2 * 3.14159265358; + double alpha = static_cast (rand()) / (static_cast (RAND_MAX/X)); + Eigen::Matrix2d rot; + rot << cos(alpha), -sin(alpha), 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(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..b1cc0fe5 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h @@ -0,0 +1,57 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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_ + +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 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::size_t> Mesh_element; + /** \brief Type of a range of mesh elements. */ + typedef std::vector Mesh_elements; + /** \brief Type of a range of scalar field . */ + typedef std::vector 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..5044660e --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h @@ -0,0 +1,174 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +namespace Gudhi { + +namespace coxeter_triangulation { + +struct Configuration { + bool toggle_edges = true, + toggle_triangles = true, + toggle_tetrahedra = true; + std::size_t ref_edges = 1, + ref_triangles = 1, + ref_tetrahedra = 1; + + 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) {} +}; + + +template +void populate_mesh(Mesh_medit& output, + Simplex_cell_map& sc_map, + Configuration configuration, + std::size_t amb_d, + std::map vi_map) { + using Mesh_element_vertices = Mesh_medit::Mesh_elements::value_type::first_type; + std::map 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 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(std::make_pair(cell, index++)); + output.vertex_points.emplace_back((1./vertex_indices.size()) * barycenter); +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(to_string(cell), vlist)); +#endif + } + + if (configuration.toggle_edges && sc_map.size() >= 2) + for (const auto& sc_map: sc_map[1]) { + Hasse_cell* edge_cell = sc_map.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(std::make_pair(edge, configuration.ref_edges)); +#ifdef GUDHI_COX_OUTPUT_TO_HTML + std::string vlist; + for (const std::size_t& v: edge) + vlist += " " + std::to_string(v); + cell_vlist_map.emplace(std::make_pair(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(std::make_pair(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 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 GUDHI_COX_OUTPUT_TO_HTML + std::string vlist = " (" + std::to_string(index) + ")"; + for (const std::size_t& v: vertex_indices) + vlist += " " + std::to_string(v); + cell_vlist_map.emplace(std::make_pair(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(std::make_pair(tetrahedron, configuration.ref_tetrahedra)); + } + index++; + } +} + +template +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 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 GUDHI_COX_OUTPUT_TO_HTML + std::string vlist; + vlist += " " + std::to_string(index); + cell_vlist_map.emplace(std::make_pair(to_string(cp_pair.first), vlist)); +#endif + vi_map.emplace(std::make_pair(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 GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(simplex, vlist)); + } +#endif + populate_mesh(output, cell_complex.boundary_simplex_cell_maps(), b_configuration, amb_d, vi_map); +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(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..43cca967 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h @@ -0,0 +1,589 @@ +#ifdef GUDHI_DEBUG +#define GUDHI_COX_OUTPUT_TO_HTML + +#include +#include +#include +#include +#include +#include + +#include + +namespace Gudhi { + +namespace coxeter_triangulation { + +template +std::ostream& operator<<(std::ostream& os, const std::vector& 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 +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 + 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_seed_inserted_list, mt_inserted_list; + +struct CC_summary_info { + std::string face_, cell_; + template + 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; +std::vector 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 faces_, post_faces_, cofaces_; + template + CC_detail_info(const Simplex_handle& simplex) + : simplex_(to_string(simplex)) {} +}; +using CC_detail_list = std::list; +std::vector cc_interior_detail_lists, cc_boundary_detail_lists; +std::vector 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 faces_; + std::size_t dimension_; + Result_type status_; + template + CC_prejoin_info(const Simplex_handle& simplex) + : simplex_(to_string(simplex)), dimension_(simplex.dimension()) {} +}; +using CC_prejoin_list = std::list; +std::vector 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 boundary_faces_; + std::list faces_, post_faces_, cofaces_; + template + CC_join_info(const Simplex_handle& simplex) + : simplex_(to_string(simplex)) {} +}; +bool join_switch = false; +std::vector cc_interior_join_detail_lists, cc_boundary_join_detail_lists; + +std::map cell_vlist_map; +std::map simplex_vlist_map; + +std::ostringstream mt_ostream, vis_ostream; +std::vector 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)"" + b_simplex + + "" + tooltiptext + ""; +} + +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)"" + b_simplex + + "" + tooltiptext + ""; +} + + +void write_head(std::ofstream& ofs) { + ofs << " \n" + << " Cell complex debug trace\n" + << " \n" + << " \n"; +} + +void write_nav(std::ofstream& ofs) { + ofs << "
\n" + << " \n" + << "
\n"; +} + +void write_mt(std::ofstream& ofs) { + ofs << "
\n"; + ofs << "

Manifold debug trace

\n"; + ofs << "

Simplices inserted during the seed phase

\n"; + ofs << "
    \n"; + for (const MT_inserted_info& mt_info: mt_seed_inserted_list) { + if (mt_info.qr_success_) { + ofs << "
  • 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_ << "
    • \n"; + } + else + ofs << "
  • Failed to insert " << mt_info.init_face_ + << "
    • \n"; + } + ofs << "
\n"; + ofs << "

Simplices inserted during the while loop phase

\n"; + ofs << "
    \n"; + for (const MT_inserted_info& mt_info: mt_inserted_list) { + if (mt_info.qr_success_) { + ofs << "
  • 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_ << "
    • \n"; + } + else + ofs << "
  • Failed to insert " << mt_info.init_face_ + << ")
    • \n"; + } + ofs << "
\n"; + ofs << "
\n"; +} + +void write_cc(std::ofstream& ofs) { + ofs << "
\n" + << "

Cell complex debug trace

\n" + << "

Go to:

\n" + << "
    \n"; + for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) { + ofs << "
  • Dimension " << i << "
    • \n"; + } + ofs << "
\n"; + for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) { + ofs << "

Dimension " << i << "

\n"; + ofs << "

Summary for interior simplices

\n"; + if (i < cc_boundary_summary_lists.size()) + ofs << "

Go to boundary

\n"; + ofs << "
    \n"; + for (const CC_summary_info& cc_info: cc_interior_summary_lists[i]) + ofs << "
  • " + << simplex_format(cc_info.face_, false) + << " cell =" << cc_info.cell_ << "
    • \n"; + ofs << "
\n"; + ofs << "

Prejoin state of the interior cells of dimension " << i << "

\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 << "
" << j << "-dimensional ambient simplices
\n"; + ofs << "
    \n"; + for (; prejoin_it->dimension_ == j; ++prejoin_it) { + ofs << "
  • " << simplex_format(prejoin_it->simplex_, false) + << " join = " << simplex_format(prejoin_it->join_, false) + << " boundary:\n" + << "
      \n"; + for (const auto& face: prejoin_it->faces_) + ofs << "
    • " << simplex_format(face) << "
      • "; + ofs << "
    \n"; + switch (prejoin_it->status_) { + case (CC_prejoin_info::Result_type::join_single): + ofs << "

    Deleted " + << simplex_format(prejoin_it->simplex_, false) + << " as it has a single face.

    "; + break; + case (CC_prejoin_info::Result_type::join_is_face): + ofs << "

    Deleted " + << simplex_format(prejoin_it->simplex_, false) + << " as its join " << simplex_format(prejoin_it->join_, false) + << " is one of the faces.

    "; + break; + case (CC_prejoin_info::Result_type::join_different): + ofs << "

    Deleted " << simplex_format(prejoin_it->simplex_, false) + << " and replaced by its join " << simplex_format(prejoin_it->join_, false) + << ".

    "; + break; + case (CC_prejoin_info::Result_type::join_same): + ofs << "

    Kept " << simplex_format(prejoin_it->simplex_, false) + << ".

    "; + } + ofs << "
    • "; + } + ofs << "
\n"; + } + ofs << "

Details for interior simplices

\n"; + ofs << "
    \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 << "
  • 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 << "
  • 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 << "
  • Insert_cell called for " << simplex_format(cc_info.simplex_, false) + << "\n"; + ofs << "
      \n"; + // for (const std::string& cof: cc_info.faces_) + // ofs << "
    • Checking if " << simplex_format(cc_info.simplex_, false) + // << " is a face of " << simplex_format(cof, false) << "\n"; + ofs << "
    \n"; + ofs << "
      \n"; + if (cc_info.status_ == CC_detail_info::Result_type::self) { + ofs << "

      The simplex " + << simplex_format(cc_info.simplex_, false) + << " already exists in the cell complex!

      \n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::face) { + ofs << "

      The simplex " + << simplex_format(cc_info.simplex_, false) << " is a face of the simplex " + << simplex_format(cc_info.trigger_, false) << "!
      \n"; + ofs << "

        \n"; + for (const std::string post_face: cc_info.post_faces_) + ofs << "
      • " + << "Post deleting " << simplex_format(post_face, false) << "
        • \n"; + ofs << "
      \n"; + ofs << "

      \n"; + ofs << "

      " + << "Deleting " << simplex_format(cc_info.trigger_, false) << "

      \n"; + } + // for (const std::string& fac: cc_info.cofaces_) + // ofs << "
    • 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 << "

      The simplex " + << simplex_format(cc_info.simplex_, false) << " is a coface of the simplex " + << simplex_format(cc_info.trigger_, false) << "!

      \n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::inserted) { + ofs << "

      Successfully inserted " + << simplex_format(cc_info.simplex_, false) << "!

      \n"; + } + ofs << "

    \n"; + ofs << "
    • \n"; + } + ofs << "
\n"; + + if (i < cc_boundary_summary_lists.size()) { + ofs << "

Summary for boundary simplices

\n"; + ofs << "

Go to interior

\n"; + ofs << "
    \n"; + for (const CC_summary_info& cc_info: cc_boundary_summary_lists[i]) + ofs << "
  • " + << simplex_format(cc_info.face_, true) + << " cell =" << cc_info.cell_ << "
    • \n"; + ofs << "
\n"; + ofs << "

Prejoin state of the boundary cells of dimension " << i << "

\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 << "
" << j << "-dimensional ambient simplices
\n"; + ofs << "
    \n"; + for (; prejoin_it->dimension_ == j; ++prejoin_it) { + ofs << "
  • " << simplex_format(prejoin_it->simplex_, true) + << " join = " << simplex_format(prejoin_it->join_, true) + << " boundary:\n" + << "
      \n"; + for (const auto& face: prejoin_it->faces_) + ofs << "
    • " << simplex_format(face) << "
      • "; + ofs << "
    \n"; + switch (prejoin_it->status_) { + case (CC_prejoin_info::Result_type::join_single): + ofs << "

    Deleted " + << simplex_format(prejoin_it->simplex_, true) + << " as it has a single face.

    "; + break; + case (CC_prejoin_info::Result_type::join_is_face): + ofs << "

    Deleted " + << simplex_format(prejoin_it->simplex_, true) + << " as its join " << simplex_format(prejoin_it->join_, true) + << " is one of the faces.

    "; + break; + case (CC_prejoin_info::Result_type::join_different): + ofs << "

    Deleted " << simplex_format(prejoin_it->simplex_, true) + << " and replaced by its join " << simplex_format(prejoin_it->join_, true) + << ".

    "; + break; + case (CC_prejoin_info::Result_type::join_same): + ofs << "

    Kept " << simplex_format(prejoin_it->simplex_, true) + << ".

    "; + } + ofs << "
    • "; + } + ofs << "
\n"; + } + } + if (i < cc_boundary_detail_lists.size()) { + ofs << "

Details for boundary simplices

\n" + << "
    \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 << "
  • 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 << "
  • 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 << "
  • Insert_simplex called on " << simplex_format(cc_info.simplex_, true); + ofs << "
      \n"; + // for (const std::string& cof: cc_info.faces_) + // ofs << "
    • Checking if " << simplex_format(cc_info.simplex_, true) + // << " is a face of " << simplex_format(cof, true) << "\n"; + ofs << "
    \n"; + ofs << "
      \n"; + if (cc_info.status_ == CC_detail_info::Result_type::self) { + ofs << "

      The simplex " + << simplex_format(cc_info.simplex_, true) + << " already exists in the cell complex!

      \n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::face) { + ofs << "

      The simplex " + << simplex_format(cc_info.simplex_, true) << " is a face of the simplex " + << simplex_format(cc_info.trigger_, true) << "!
      \n"; + ofs << "

        \n"; + for (const std::string post_face: cc_info.post_faces_) + ofs << "
      • Post deleting " << simplex_format(post_face, true) + << "
        • \n"; + ofs << "
      \n"; + ofs << "

      \n"; + ofs << "

      Deleting " << simplex_format(cc_info.trigger_, true) << "

      \n"; + } + // for (const std::string& fac: cc_info.cofaces_) + // ofs << "
    • Checking if " << simplex_format(cc_info.simplex_, true) + // << " is a coface of " << simplex_format(fac, true) << "\n"; + ofs << "
    \n"; + ofs << "
    • \n"; + if (cc_info.status_ == CC_detail_info::Result_type::coface) { + ofs << "

    The simplex " + << simplex_format(cc_info.simplex_, true) << " is a coface of the simplex " + << simplex_format(cc_info.trigger_, true) << "!

    \n"; + } + if (cc_info.status_ == CC_detail_info::Result_type::inserted) { + ofs << "

    Successfully inserted " + << simplex_format(cc_info.simplex_, true) << "!

    \n"; + } + } + ofs << "

\n"; + } + } + ofs << "
\n"; +} + +void write_visu(std::ofstream& ofs) { + ofs << "
\n" + << "

Visualization details debug trace

\n"; + // std::vector > vs_maps(cc_interior_summary_lists.size()); + std::map vs_map; + for (const auto& sv_pair: simplex_vlist_map) + vs_map.emplace(std::make_pair(sv_pair.second, sv_pair.first)); + ofs << "
    \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 << "
  • " << vs_pair.first << ": " + << simplex_format(w_simplex, is_boundary) << "
    • \n"; + } + ofs << "
\n"; + ofs << "
\n"; +} + +void write_to_html(std::string file_name) { + std::ofstream ofs(file_name + ".html", std::ofstream::out); + ofs << "\n" + << "\n"; + write_head(ofs); + ofs << " \n"; + write_nav(ofs); + ofs << "

Debug traces for " << file_name << "

\n"; + write_mt(ofs); + write_cc(ofs); + write_visu(ofs); + ofs << " \n"; + ofs << "\n"; + + ofs.close(); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif + 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..fe379242 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h @@ -0,0 +1,181 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include + +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 +typename std::enable_if::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 +typename std::enable_if::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(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 edge; + for (const auto& v_i: e.first) + edge.push_back(v_i + index); + edges.emplace_back(std::make_pair(edge, e.second)); + } + for (const auto& t: mesh.triangles) { + std::vector triangle; + for (const auto& v_i: t.first) + triangle.push_back(v_i + index); + triangles.emplace_back(std::make_pair(triangle, t.second)); + } + for (const auto& t: mesh.tetrahedra) { + std::vector tetrahedron; + for (const auto& v_i: t.first) + tetrahedron.push_back(v_i + index); + tetrahedra.emplace_back(std::make_pair(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(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 Medit. + * + * @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. + */ +template +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..a522f691 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h @@ -0,0 +1,279 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include +#include +#include + +#include + +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 Implicit_manifold_intersection_oracle { + + /* Computes the affine coordinates of the intersection point of the implicit manifold + * and the affine hull of the simplex. */ + template + 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); + 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 + 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); + return lambda; + } + + /* Computes the intersection result for a given simplex in a triangulation. */ + template + Query_result intersection_result(const Eigen::VectorXd& lambda, + const Simplex_handle& simplex, + const Triangulation& triangulation) const { + using QR = Query_result; + 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 (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 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 + Query_result 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 + Query_result intersects_boundary(const Simplex_handle& simplex, + const Triangulation& triangulation) const { + 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 + 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 +Implicit_manifold_intersection_oracle +make_oracle(const Function_& function, + const Domain_function_& domain_function){ + return Implicit_manifold_intersection_oracle(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 +Implicit_manifold_intersection_oracle make_oracle(const Function_& function){ + return Implicit_manifold_intersection_oracle(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..fdb4f630 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h @@ -0,0 +1,307 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include + +#include + +#include +#include + +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 Manifold_tracing { + + typedef typename Triangulation_::Simplex_handle 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()(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 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 + 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 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 qr = oracle.intersects(face, triangulation); + if (qr.success && + out_simplex_map.emplace(std::make_pair(face, qr.intersection)).second) { +#ifdef GUDHI_COX_OUTPUT_TO_HTML + 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 qr = oracle.intersects(face, triangulation); + if (qr.success && + out_simplex_map.emplace(std::make_pair(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 + 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 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 GUDHI_COX_OUTPUT_TO_HTML + 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(std::make_pair(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 GUDHI_COX_OUTPUT_TO_HTML + 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 GUDHI_COX_OUTPUT_TO_HTML + 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 GUDHI_COX_OUTPUT_TO_HTML + mt_inserted_list.push_back(MT_inserted_info(qrb, cof, true)); +#endif + // assert (qrb.success); // always a success + if (qrb.success) + boundary_simplex_map.emplace(cof, qrb.intersection); + } + } + } + } + } + } + + /** \brief Empty constructor */ + Manifold_tracing() {} + +}; + +/** + * \brief Static method for Manifold_tracing::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::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 +void manifold_tracing_algorithm(const Point_range& seed_points, + const Triangulation& triangulation, + const Intersection_oracle& oracle, + Out_simplex_map& out_simplex_map) { + Manifold_tracing mt; + mt.manifold_tracing_algorithm(seed_points, + triangulation, + oracle, + out_simplex_map); +} + +/** + * \brief Static method for Manifold_tracing::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::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 +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 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..a3c2b2c7 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h @@ -0,0 +1,257 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + +#include +#include + +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 Permutahedral_representation { + + typedef Permutahedral_representation 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 Vertex_iterator; + typedef boost::iterator_range 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 Face_iterator; + typedef boost::iterator_range 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 Coface_iterator; + typedef boost::iterator_range 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 +std::ostream& operator<<(std::ostream& os, + const Permutahedral_representation& 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..20a4a8ff --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h @@ -0,0 +1,101 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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 const, + boost::forward_traversal_tag> { + typedef std::vector 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; + } + } + +protected: + 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..47eb8b98 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h @@ -0,0 +1,128 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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 const, + boost::forward_traversal_tag> { + typedef std::vector value_t; + +protected: + 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 + Integer_combination_iterator(const uint& n, const uint& k, const Bound_range& bounds) + : + value_(k+2), + is_end_(n == 0 || k == 0), + n_(n), + 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), n_(0), k_(0) {} + +protected: + value_t value_; // the dereference value + bool is_end_; // is true when the current integer combination is the final one + + uint n_; + uint k_; + std::vector 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..55c32664 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h @@ -0,0 +1,108 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include +#include +#include + +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 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> { + typedef Ordered_set_partition value_t; + +protected: + 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(); + } + +protected: + 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..4528ef20 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h @@ -0,0 +1,288 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include +#include +#include +#include + +#include +#include +#include + +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 Vertex_iterator : public boost::iterator_facade< Vertex_iterator, + typename Permutahedral_representation::Vertex const, + boost::forward_traversal_tag> { +protected: + friend class boost::iterator_core_access; + + using Vertex = typename Permutahedral_representation::Vertex; + using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition; + + typedef Vertex value_t; + + + 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) {} + + +protected: + 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 Face_iterator : public boost::iterator_facade< Face_iterator, + Permutahedral_representation const, + boost::forward_traversal_tag> { + typedef Permutahedral_representation value_t; + +protected: + 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(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) {} + +protected: + 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 Coface_iterator : public boost::iterator_facade< Coface_iterator, + Permutahedral_representation const, + boost::forward_traversal_tag> { + typedef Permutahedral_representation value_t; + +protected: + 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(Ordered_set_partition_iterator(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(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(Ordered_set_partition_iterator(simplex_.partition()[i].size(), (*i_it_)[i]+1)); + update_value(); + } + + // Used for the creating an end iterator + Coface_iterator() : is_end_(true) {} + +protected: + 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 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..7cf6158b --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h @@ -0,0 +1,137 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include + +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 const, + boost::forward_traversal_tag> { + typedef std::vector value_t; + +protected: + 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; + } + +protected: + 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 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..46f67752 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h @@ -0,0 +1,137 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include +#include + +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 > const, + boost::forward_traversal_tag> { + typedef std::vector > value_t; + +protected: + 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::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(); + } + +protected: + value_t value_; // the dereference value + std::vector rgs_; // restricted growth string + std::vector 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..5b3c29e9 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h @@ -0,0 +1,64 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 + * + * \ingroup coxeter_triangulation + */ +template +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..dd9c20dc --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h @@ -0,0 +1,89 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +#include // std::size_t + +namespace Gudhi { + +namespace coxeter_triangulation { + +/** \brief Auxillary iterator class for sizes of parts in an ordered set partition. + */ +template +class Size_iterator : public boost::iterator_facade< Size_iterator, + std::size_t const, + boost::forward_traversal_tag> { + friend class boost::iterator_core_access; + +protected: + 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; + } + +protected: + T_it t_it_, t_end_; + bool is_end_; + std::size_t value_; +}; + +template +class Size_range { + const T& t_; + +public: + typedef Size_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..461a01e7 --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h @@ -0,0 +1,71 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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_ + +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 +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 diff --git a/src/Coxeter_triangulation/include/gudhi/Query_result.h b/src/Coxeter_triangulation/include/gudhi/Query_result.h new file mode 100644 index 00000000..c7384c0b --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Query_result.h @@ -0,0 +1,42 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * 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 +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 -- cgit v1.2.3