summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorVincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com>2021-11-02 13:26:16 +0100
committerGitHub <noreply@github.com>2021-11-02 13:26:16 +0100
commit07333d6e217eca028c41886e1c17dd78477adcf4 (patch)
tree0530c8987c824da79ecddd22485f143831a4c71f /src
parent8adb46d8a54f1a0dd71ea686473cc4ca9f5d2f67 (diff)
parent003abb66b7c657356328bc24306e2c5aef02b1d4 (diff)
Merge pull request #397 from VincentRouvreau/coxeter_integration
Coxeter integration
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt1
-rw-r--r--src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h46
-rw-r--r--src/Coxeter_triangulation/concept/IntersectionOracle.h104
-rw-r--r--src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h81
-rw-r--r--src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h56
-rw-r--r--src/Coxeter_triangulation/doc/custom_function.pngbin0 -> 256301 bytes
-rw-r--r--src/Coxeter_triangulation/doc/flat_torus_with_boundary.pngbin0 -> 222900 bytes
-rw-r--r--src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h240
-rw-r--r--src/Coxeter_triangulation/doc/manifold_tracing_on_custom_function_example.pngbin0 -> 589120 bytes
-rw-r--r--src/Coxeter_triangulation/doc/two_triangulations.pngbin0 -> 39507 bytes
-rw-r--r--src/Coxeter_triangulation/example/CMakeLists.txt19
-rw-r--r--src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp55
-rw-r--r--src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt26
-rw-r--r--src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp87
-rw-r--r--src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp72
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h77
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h340
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h285
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h40
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h219
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h157
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h64
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h93
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h110
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h91
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h80
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h69
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h85
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h79
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h71
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h78
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h88
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Negation.h84
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h111
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/Translate.h89
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h72
-rw-r--r--src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h60
-rw-r--r--src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h171
-rw-r--r--src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h550
-rw-r--r--src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h154
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h261
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h270
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h216
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h83
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h114
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h93
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h254
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h120
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h111
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h54
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h73
-rw-r--r--src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h66
-rw-r--r--src/Coxeter_triangulation/test/CMakeLists.txt30
-rw-r--r--src/Coxeter_triangulation/test/cell_complex_test.cpp59
-rw-r--r--src/Coxeter_triangulation/test/freud_triang_test.cpp114
-rw-r--r--src/Coxeter_triangulation/test/function_test.cpp158
-rw-r--r--src/Coxeter_triangulation/test/manifold_tracing_test.cpp62
-rw-r--r--src/Coxeter_triangulation/test/oracle_test.cpp56
-rw-r--r--src/Coxeter_triangulation/test/perm_rep_test.cpp61
-rw-r--r--src/Coxeter_triangulation/test/random_orthogonal_matrix_function_test.cpp36
-rw-r--r--src/common/doc/examples.h188
-rw-r--r--src/common/doc/installation.h8
-rw-r--r--src/common/doc/main_page.md28
63 files changed, 6526 insertions, 93 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 79ec42c1..8f6a1ccc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,6 +27,7 @@ add_gudhi_module(Bottleneck_distance)
add_gudhi_module(Cech_complex)
add_gudhi_module(Contraction)
add_gudhi_module(Collapse)
+add_gudhi_module(Coxeter_triangulation)
add_gudhi_module(Hasse_complex)
add_gudhi_module(Persistence_representations)
add_gudhi_module(Persistent_cohomology)
diff --git a/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h b/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h
new file mode 100644
index 00000000..210d804e
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h
@@ -0,0 +1,46 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_FUNCTION_FOR_IMPLICIT_MANIFOLD_H_
+#define CONCEPT_COXETER_TRIANGULATION_FUNCTION_FOR_IMPLICIT_MANIFOLD_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept FunctionForImplicitManifold describes the requirements
+ * for a type to implement an implicit function class used for example in Manifold_tracing.
+ */
+struct FunctionForImplicitManifold {
+ /** \brief Value of the function at a specified point 'p'.
+ * @param[in] p The input point given by its Cartesian coordinates.
+ * Its size needs to be equal to amb_d().
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const;
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const;
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const;
+
+ /** \brief Returns a point on the zero-set of the function. */
+ Eigen::VectorXd seed() const;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/IntersectionOracle.h b/src/Coxeter_triangulation/concept/IntersectionOracle.h
new file mode 100644
index 00000000..e4e397fa
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/IntersectionOracle.h
@@ -0,0 +1,104 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_INTERSECTION_ORACLE_H_
+#define CONCEPT_COXETER_TRIANGULATION_INTERSECTION_ORACLE_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept IntersectionOracle describes the requirements
+ * for a type to implement an intersection oracle class used for example in Manifold_tracing.
+ *
+ */
+struct IntersectionOracle {
+ /** \brief Returns the domain (ambient) dimension of the underlying manifold. */
+ std::size_t amb_d() const;
+
+ /** \brief Returns the codomain dimension of the underlying manifold. */
+ std::size_t cod_d() const;
+
+ /** \brief Intersection query with the relative interior of the manifold.
+ *
+ * \details The returned structure Query_result contains the boolean value
+ * that is true only if the intersection point of the query simplex and
+ * the relative interior of the manifold exists, the intersection point
+ * and the face of the query simplex that contains
+ * the intersection point.
+ *
+ * \tparam Simplex_handle The class of the query simplex.
+ * Needs to be a model of the concept SimplexInCoxeterTriangulation.
+ * \tparam Triangulation The class of the triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ *
+ * @param[in] simplex The query simplex. The dimension of the simplex
+ * should be the same as the codimension of the manifold
+ * (the codomain dimension of the function).
+ * @param[in] triangulation The ambient triangulation. The dimension of
+ * the triangulation should be the same as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Simplex_handle, class Triangulation>
+ Query_result<Simplex_handle> intersects(const Simplex_handle& simplex, const Triangulation& triangulation) const;
+
+ /** \brief Intersection query with the boundary of the manifold.
+ *
+ * \details The returned structure Query_result contains the boolean value
+ * that is true only if the intersection point of the query simplex and
+ * the boundary of the manifold exists, the intersection point
+ * and the face of the query simplex that contains
+ * the intersection point.
+ *
+ * \tparam Simplex_handle The class of the query simplex.
+ * Needs to be a model of the concept SimplexInCoxeterTriangulation.
+ * \tparam Triangulation The class of the triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ *
+ * @param[in] simplex The query simplex. The dimension of the simplex
+ * should be the same as the codimension of the boundary of the manifold
+ * (the codomain dimension of the function + 1).
+ * @param[in] triangulation The ambient triangulation. The dimension of
+ * the triangulation should be the same as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Simplex_handle, class Triangulation>
+ Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex,
+ const Triangulation& triangulation) const;
+
+ /** \brief Returns true if the input point lies inside the piecewise-linear
+ * domain induced by the given ambient triangulation that defines the relative
+ * interior of the piecewise-linear approximation of the manifold.
+ *
+ * @param p The input point. Needs to have the same dimension as the ambient
+ * dimension of the manifold (the domain dimension of the function).
+ * @param triangulation The ambient triangulation. Needs to have the same
+ * dimension as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Triangulation>
+ bool lies_in_domain(const Eigen::VectorXd& p, const Triangulation& triangulation) const {
+ Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p);
+ return pl_p(0) < 0;
+ }
+
+ /** \brief Returns the function that defines the interior of the manifold */
+ const Function_& function() const;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h b/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h
new file mode 100644
index 00000000..dac8e66d
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h
@@ -0,0 +1,81 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_SIMPLEX_IN_COXETER_TRIANGULATION_H_
+#define CONCEPT_COXETER_TRIANGULATION_SIMPLEX_IN_COXETER_TRIANGULATION_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <gudhi/Permutahedral_representation.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept SimplexInCoxeterTriangulation describes the requirements
+ * for a type to implement a representation of simplices in Freudenthal_triangulation
+ * or in Coxeter_triangulation.
+ */
+struct SimplexInCoxeterTriangulation {
+ /** \brief Type of the vertex. */
+ typedef Vertex_ Vertex;
+
+ /** \brief Type of the ordered partition. */
+ typedef Ordered_set_partition_ OrderedSetPartition;
+
+ /** \brief Dimension of the simplex. */
+ unsigned dimension() const;
+
+ /** \brief Type of a range of vertices, each of type Vertex. */
+ typedef Vertex_range;
+
+ /** \brief Returns a range of vertices of the simplex.
+ */
+ Vertex_range vertex_range() const;
+
+ /** \brief Type of a range of faces, each of type that
+ * is a model of the concept SimplexInCoxeterTriangulation.
+ */
+ typedef Face_range;
+
+ /** \brief Returns a range of permutahedral representations of k-dimensional faces
+ * of the simplex for some given integer parameter 'k'.
+ */
+ Face_range face_range(std::size_t k) const;
+
+ /** \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;
+
+ /** \brief Type of a range of cofaces, each of type that
+ * is a model of the concept SimplexInCoxeterTriangulation.
+ */
+ typedef Coface_range;
+
+ /** \brief Returns a range of permutahedral representations of k-dimensional cofaces
+ * of the simplex for some given integer parameter 'k'.
+ */
+ Coface_range coface_range(std::size_t k) const;
+
+ /** \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.
+ */
+ Coface_range cofacet_range() const;
+
+ /** \brief Returns true, if the simplex is a face of other simplex. */
+ bool is_face_of(const Permutahedral_representation& other) const;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h b/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h
new file mode 100644
index 00000000..2b5d568c
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h
@@ -0,0 +1,56 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_TRIANGULATION_FOR_MANIFOLD_TRACING_H_
+#define CONCEPT_COXETER_TRIANGULATION_TRIANGULATION_FOR_MANIFOLD_TRACING_H_
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept TriangulationForManifoldTracing describes the requirements
+ * for a type to implement a triangulation class used for example in Manifold_tracing.
+ */
+struct TriangulationForManifoldTracing {
+ /** \brief Type of the simplices in the triangulation.
+ * Needs to be a model of the concept SimplexInCoxeterTriangulation. */
+ typedef Simplex_handle;
+
+ /** \brief Type of the vertices in the triangulation.
+ * Needs to be a random-access range of integer values. */
+ typedef Vertex_handle;
+
+ /** \brief Returns the permutahedral representation of the simplex in the
+ * triangulation that contains a given query point 'p'.
+ * \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.
+ */
+ template <class Point_d>
+ Simplex_handle locate_point(const Point_d& point) const;
+
+ /** \brief Returns the Cartesian coordinates of the given vertex 'v'.
+ * @param[in] v The input vertex.
+ */
+ Eigen::VectorXd cartesian_coordinates(const Vertex_handle& v) const;
+
+ /** \brief Returns the Cartesian coordinates of the barycenter of a given simplex 's'.
+ * @param[in] s The input simplex given by permutahedral representation.
+ */
+ Eigen::VectorXd barycenter(const Simplex_handle& s) const;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/doc/custom_function.png b/src/Coxeter_triangulation/doc/custom_function.png
new file mode 100644
index 00000000..8bb8ba9a
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/custom_function.png
Binary files differ
diff --git a/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png
new file mode 100644
index 00000000..338b39fe
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png
Binary files differ
diff --git a/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h b/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h
new file mode 100644
index 00000000..395996c9
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h
@@ -0,0 +1,240 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef DOC_COXETER_TRIANGULATION_INTRO_COXETER_TRIANGULATION_H_
+#define DOC_COXETER_TRIANGULATION_INTRO_COXETER_TRIANGULATION_H_
+
+// needs namespaces for Doxygen to link on classes
+namespace Gudhi {
+namespace coxeter_triangulation {
+
+/** \defgroup coxeter_triangulation Coxeter triangulation
+
+\author Siargey Kachanovich
+
+@{
+
+\section overview Module overview
+
+Coxeter triangulation module is designed to provide tools for constructing a piecewise-linear approximation of an
+\f$m\f$-dimensional smooth manifold embedded in \f$ \mathbb{R}^d \f$ using an ambient triangulation.
+For a more detailed description of the module see \cite KachanovichThesis.
+
+\section manifoldtracing Manifold tracing algorithm
+The central piece of the module is the manifold tracing algorithm represented by the class
+\ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing".
+The manifold tracing algorithm takes as input a manifold of some dimension \f$m\f$ embedded in \f$\mathbb{R}^d\f$
+represented by an intersection oracle (see Section \ref intersectionoracle "Intersection oracle"), a point on the
+manifold and an ambient triangulation (see Section \ref ambienttriangulations "Ambient triangulations").
+The output consists of one map (or two maps in the case of manifolds with boundary) from the \f$(d-m)\f$-dimensional
+(and \f$(d-m+1)\f$-dimensional in the case of manifolds with boundary) simplices in the ambient triangulation that
+intersect the manifold to their intersection points.
+From this output, it is possible to construct the cell complex of the piecewise-linear approximation of the input
+manifold.
+
+There are two methods that execute the manifold tracing algorithm: the method
+\ref Gudhi::coxeter_triangulation::Manifold_tracing::manifold_tracing_algorithm() "Manifold_tracing::manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map)"
+for manifolds without boundary and
+\ref Gudhi::coxeter_triangulation::Manifold_tracing::manifold_tracing_algorithm() "Manifold_tracing::manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map,boundary_simplex_map)"
+for manifolds with boundary. The algorithm functions as follows. It starts at the specified seed points and inserts a
+\f$(d-m)\f$-dimensional simplices nearby each seed point that intersect the manifold into the output. Starting from
+this simplex, the algorithm propagates the search for other \f$(d-m)\f$-dimensional simplices that intersect the
+manifold by marching from a simplex to neighbouring simplices via their common cofaces.
+
+This class \ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing" has one template parameter
+`Triangulation_` which specifies the ambient triangulation which is used by the algorithm.
+The template type `Triangulation_` has to be a model of the concept
+\ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing".
+
+The module also provides two static methods:
+\ref Gudhi::coxeter_triangulation::manifold_tracing_algorithm() "manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map)"
+for manifolds without boundary and
+\ref manifold_tracing_algorithm() "manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map)"
+for manifolds with boundary. For these static methods it is not necessary to specify any template arguments.
+
+\section ambienttriangulations Ambient triangulations
+
+The ambient triangulations supported by the manifold tracing algorithm have to be models of the concept
+\ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing".
+This module offers two such models: the class
+\ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" and the derived class
+\ref Gudhi::coxeter_triangulation::Coxeter_triangulation "Coxeter_triangulation".
+
+Both these classes encode affine transformations of the so-called Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$.
+The Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$ is defined as the simplicial subdivision of the unit cubic
+partition of \f$\mathbb{R}^d\f$.
+Each simplex is encoded using the permutahedral representation, which consists of an integer-valued vector \f$y\f$ that
+positions the simplex in a specific cube in the cubical partition and an ordered partition \f$\omega\f$ of the set
+\f$\{1,\ldots,d+1\}\f$, which positions the simplex in the simplicial subdivision of the cube.
+The default constructor
+\ref Gudhi::coxeter_triangulation::Freudenthal_triangulation::Freudenthal_triangulation(std::size_t)
+"Freudenthal_triangulation(d)" the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$. The class
+\ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" can also encode any affine
+transformation of the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$ using an invertible matrix \f$\Lambda\f$ and
+an offset vector \f$b\f$ that can be specified in the constructor and which can be changed using the methods
+change_matrix and change_offset. The class
+\ref Gudhi::coxeter_triangulation::Coxeter_triangulation "Coxeter_triangulation" is derived from
+\ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" and its default constructor
+\ref Gudhi::coxeter_triangulation::Coxeter_triangulation::Coxeter_triangulation(std::size_t) "Coxeter_triangulation(d)"
+builds a Coxeter triangulation of type \f$\tilde{A}_d\f$, which has the best simplex quality of all linear
+transformations of the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$.
+
+\image html two_triangulations.png "Coxeter (on the left) and Freudenthal-Kuhn triangulation (on the right)"
+
+
+\section intersectionoracle Intersection oracle
+
+The input \f$m\f$-dimensional manifold in \f$\mathbb{R}^d\f$ needs to be given via the intersection oracle that answers
+the following query: given a \f$(d-m)\f$-dimensional simplex, does it intersect the manifold?
+The concept \ref Gudhi::coxeter_triangulation::IntersectionOracle "IntersectionOracle" describes all requirements for
+an intersection oracle class to be compatible with the class
+\ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing".
+This module offers one model of the concept
+\ref Gudhi::coxeter_triangulation::IntersectionOracle "IntersectionOracle", which is the class
+\ref Gudhi::coxeter_triangulation::Implicit_manifold_intersection_oracle "Implicit_manifold_intersection_oracle".
+This class represents a manifold given as the zero-set of a specified function
+\f$F: \mathbb{R}^d \rightarrow \mathbb{R}^{d-m}\f$.
+The function \f$F\f$ is given by a class which is a model of the concept
+\ref Gudhi::coxeter_triangulation::FunctionForImplicitManifold "FunctionForImplicitManifold".
+There are multiple function classes that are already implemented in this module.
+
+\li \ref Gudhi::coxeter_triangulation::Constant_function(std::size_t, std::size_t, Eigen::VectorXd)
+"Constant_function(d,k,v)" defines a constant function \f$F\f$ such that for all \f$x \in \mathbb{R}^d\f$, we have
+ \f$F(x) = v \in \mathbb{R}^k\f$.
+ The class Constant_function does not define an implicit manifold, but is useful as the domain function when defining
+ boundaryless implicit manifolds.
+\li \ref Gudhi::coxeter_triangulation::Function_affine_plane_in_Rd(N,b) "Function_affine_plane_in_Rd(N,b)" defines an
+ \f$m\f$-dimensional implicit affine plane in the \f$d\f$-dimensional Euclidean space given by a normal matrix \f$N\f$
+ and an offset vector \f$b\f$.
+\li \ref Gudhi::coxeter_triangulation::Function_Sm_in_Rd(r,m,d,center) "Function_Sm_in_Rd(r,m,d,center)" defines an
+ \f$m\f$-dimensional implicit sphere embedded in the \f$d\f$-dimensional Euclidean space of radius \f$r\f$ centered at
+ the point 'center'.
+\li \ref Gudhi::coxeter_triangulation::Function_moment_curve_in_Rd(r,d) "Function_moment_curve(r,d)" defines the moment
+ curve in the \f$d\f$-dimensional Euclidean space of radius \f$r\f$ given as the parameterized curve (but implemented
+ as an implicit curve):
+ \f[ (r, rt, \ldots, rt^{d-1}) \in \mathbb{R}^d,\text{ for $t \in \mathbb{R}$.} \f]
+\li \ref Gudhi::coxeter_triangulation::Function_torus_in_R3(R, r) "Function_torus_in_R3(R, r)" defines a torus in
+ \f$\mathbb{R}^3\f$ with the outer radius \f$R\f$ and the inner radius, given by the equation:
+ \f[ z^2 + (\sqrt{x^2 + y^2} - r)^2 - R^2 = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_chair_in_R3(a, b, k) "Function_chair_in_R3(a, b, k)" defines the
+ \"Chair\" surface in \f$\mathbb{R}^3\f$ defined by the equation:
+ \f[ (x^2 + y^2 + z^2 - ak^2)^2 - b((z-k)^2 - 2x^2)((z+k)^2 - 2y^2) = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_iron_in_R3() "Function_iron_in_R3()" defines the \"Iron\" surface in
+ \f$\mathbb{R}^3\f$ defined by the equation:
+ \f[ \frac{-x^6-y^6-z^6}{300} + \frac{xy^2z}{2.1} + y^2 + (z-2)^2 = 1. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_lemniscate_revolution_in_R3(a) "Function_lemniscate_revolution_in_R3(a)"
+ defines a revolution surface in \f$\mathbb{R}^3\f$ obtained from the lemniscate of Bernoulli defined by the equation:
+ \f[ (x^2 + y^2 + z^2)^2 - 2a^2(x^2 - y^2 - z^2) = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_whitney_umbrella_in_R3() "Function_whitney_umbrella_in_R3()" defines
+ the Whitney umbrella surface in \f$\mathbb{R}^3\f$ defined by the equation:
+ \f[ x^2 - y^2z = 0. \f]
+
+The base function classes above can be composed or modified into new functions using the following classes and methods:
+
+\li \ref Gudhi::coxeter_triangulation::Cartesian_product "Cartesian_product(functions...)" expresses the Cartesian
+ product \f$F_1^{-1}(0) \times \ldots \times F_k^{-1}(0)\f$ of multiple implicit manifolds as an implicit manifold.
+ For convenience, a static function
+ \ref Gudhi::coxeter_triangulation::make_product_function() "make_product_function(functions...)" is provided that
+ takes a pack of function-typed objects as the argument.
+\li \ref Gudhi::coxeter_triangulation::Embed_in_Rd "Embed_in_Rd(F, d)" expresses an implicit manifold given as the
+ zero-set of a function \f$F\f$ embedded in a higher-dimensional Euclidean space \f$\mathbb{R}^d\f$.
+ For convenience, a static function \ref Gudhi::coxeter_triangulation::make_embedding() "make_embedding(F, d)" is
+ provided.
+\li \ref Gudhi::coxeter_triangulation::Linear_transformation "Linear_transformation(F, M)" applies a linear
+ transformation given by a matrix \f$M\f$ on an implicit manifold given as the zero-set of the function \f$F\f$.
+ For convenience, a static function
+ \ref Gudhi::coxeter_triangulation::make_linear_transformation() "make_linear_transformation(F, M)" is provided.
+\li \ref Gudhi::coxeter_triangulation::Translate "Translate(F, v)" translates an implicit manifold given as the
+ zero-set of ththe function \f$F\f$ by a vector \f$v\f$.
+ For convenience, a static function \ref Gudhi::coxeter_triangulation::translate() "translate(F, v)" is provided.
+\li \ref Gudhi::coxeter_triangulation::Negation() "Negation(F)" defines the negative of the given function \f$F\f$.
+ This class is useful to define the complementary of a given domain, when defining a manifold with boundary.
+ For convenience, a static function \ref Gudhi::coxeter_triangulation::negation() "negation(F)" is provided.
+\li \ref Gudhi::coxeter_triangulation::PL_approximation "PL_approximation(F, T)" defines a piecewise-linear
+ approximation of a given function \f$F\f$ induced by an ambient triangulation \f$T\f$.
+ The purpose of this class is to define a piecewise-linear function that is compatible with the requirements for the
+ domain function \f$D\f$ when defining a manifold with boundary.
+ For convenience, a static function
+ \ref Gudhi::coxeter_triangulation::make_pl_approximation() "make_pl_approximation(F, T)" is provided.
+ The type of \f$T\f$ is required to be a model of the concept
+ \ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing".
+
+It is also possible to implement your own function as detailed in this \ref exampleswithcustomfunction.
+
+\section cellcomplex Cell complex construction
+
+The output of the manifold tracing algorithm can be transformed into the Hasse diagram of a cell complex that
+approximates the input manifold using the class \ref Gudhi::coxeter_triangulation::Cell_complex "Cell_complex".
+The type of the cells in the Hasse diagram is
+\ref Gudhi::Hasse_diagram::Hasse_diagram_cell "Hasse_cell<int, double, bool>" provided by the module Hasse diagram.
+The cells in the cell complex given by an object of the class
+\ref Gudhi::coxeter_triangulation::Cell_complex "Cell_complex" are accessed through several maps that are accessed
+through the following methods.
+
+\li The method
+\ref Gudhi::coxeter_triangulation::Cell_complex::interior_simplex_cell_maps() "interior_simplex_cell_maps()"
+returns a vector of maps from the cells of various dimensions in the interior of the cell complex to the permutahedral
+representations of the corresponding simplices in the ambient triangulation.
+Each individual map for cells of a specific dimension \f$l\f$ can be accessed using the method
+\ref Gudhi::coxeter_triangulation::Cell_complex::interior_simplex_cell_map() "interior_simplex_cell_map(l)".
+\li The method
+\ref Gudhi::coxeter_triangulation::Cell_complex::boundary_simplex_cell_maps() "boundary_simplex_cell_maps()"
+returns a vector of maps from the cells of various dimensions on the boundary of the cell complex to the permutahedral
+representations of the corresponding simplices in the ambient triangulation.
+Each individual map for cells of a specific dimension \f$l\f$ can be accessed using the method
+\ref Gudhi::coxeter_triangulation::Cell_complex::boundary_simplex_cell_map() "boundary_simplex_cell_map(l)".
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::cell_simplex_map() "cell_simplex_map()" returns a map
+from the cells in the cell complex to the permutahedral representations of the corresponding simplices in the ambient
+triangulation.
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::cell_point_map() "cell_point_map()" returns a map from
+the vertex cells in the cell complex to their Cartesian coordinates.
+
+The use and interfaces of this \ref Gudhi::coxeter_triangulation::Cell_complex "Cell_complex" is limited to the
+Coxeter_triangulation implementation.
+
+\section example Examples
+
+\subsection examplewithoutboundaries Basic example without boundaries
+\include cell_complex_from_basic_circle_manifold.cpp
+
+The program output is:
+
+\include cell_complex_from_basic_circle_manifold_for_doc.txt
+
+\subsection exampleswithboundaries Example with boundaries
+
+Here is an example of constructing a piecewise-linear approximation of a flat torus embedded in \f$\mathbb{R}^4\f$,
+rotated by a random rotation in \f$\mathbb{R}^4\f$ and cut by a hyperplane.
+
+\include manifold_tracing_flat_torus_with_boundary.cpp
+
+The output in <a target="_blank" href="https://www.ljll.math.upmc.fr/frey/software.html">medit</a> is:
+
+\image html "flat_torus_with_boundary.png" "Output from the example of a flat torus with boundary"
+
+\subsection exampleswithcustomfunction Example with a custom function
+
+In the following more complex example, we define a custom function for the implicit manifold.
+
+\include manifold_tracing_custom_function.cpp
+
+The output in <a target="_blank" href="https://www.ljll.math.upmc.fr/frey/software.html">medit</a> looks as follows:
+
+\image html "custom_function.png" "Output from the example with a custom function"
+
+
+ */
+/** @} */ // end defgroup coxeter_triangulation
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif // DOC_COXETER_TRIANGULATION_INTRO_COXETER_TRIANGULATION_H_
diff --git a/src/Coxeter_triangulation/doc/manifold_tracing_on_custom_function_example.png b/src/Coxeter_triangulation/doc/manifold_tracing_on_custom_function_example.png
new file mode 100644
index 00000000..04912729
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/manifold_tracing_on_custom_function_example.png
Binary files differ
diff --git a/src/Coxeter_triangulation/doc/two_triangulations.png b/src/Coxeter_triangulation/doc/two_triangulations.png
new file mode 100644
index 00000000..055d93e7
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/two_triangulations.png
Binary files differ
diff --git a/src/Coxeter_triangulation/example/CMakeLists.txt b/src/Coxeter_triangulation/example/CMakeLists.txt
new file mode 100644
index 00000000..7f81c599
--- /dev/null
+++ b/src/Coxeter_triangulation/example/CMakeLists.txt
@@ -0,0 +1,19 @@
+project(Coxeter_triangulation_example)
+
+if (NOT EIGEN3_VERSION VERSION_LESS 3.1.0)
+ # because of random_orthogonal_matrix inclusion
+ if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+ add_executable ( Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example manifold_tracing_flat_torus_with_boundary.cpp )
+ target_link_libraries(Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example ${CGAL_LIBRARY})
+ add_test(NAME Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example
+ COMMAND $<TARGET_FILE:Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example>)
+ endif()
+
+ add_executable ( Coxeter_triangulation_manifold_tracing_custom_function_example manifold_tracing_custom_function.cpp )
+ add_test(NAME Coxeter_triangulation_manifold_tracing_custom_function_example
+ COMMAND $<TARGET_FILE:Coxeter_triangulation_manifold_tracing_custom_function_example>)
+
+ add_executable ( Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example cell_complex_from_basic_circle_manifold.cpp )
+ add_test(NAME Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example
+ COMMAND $<TARGET_FILE:Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example>)
+endif() \ No newline at end of file
diff --git a/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp
new file mode 100644
index 00000000..dfaaffa8
--- /dev/null
+++ b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp
@@ -0,0 +1,55 @@
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h> // for Gudhi::coxeter_triangulation::make_oracle
+#include <gudhi/Manifold_tracing.h>
+#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+int main(int argc, char** argv) {
+ // Oracle is a circle of radius 1
+ double radius = 1.;
+ auto oracle = make_oracle(Function_Sm_in_Rd(radius, 1));
+
+ // Define a Coxeter triangulation.
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ // Theory forbids that a vertex of the triangulation lies exactly on the circle.
+ // Add some offset to avoid algorithm degeneracies.
+ cox_tr.change_offset(-Eigen::VectorXd::Random(oracle.amb_d()));
+ // For a better manifold approximation, one can change the circle radius value or change the linear transformation
+ // matrix.
+ // The number of points and edges will increase with a better resolution.
+ //cox_tr.change_matrix(0.5 * cox_tr.matrix());
+
+ // Manifold tracing algorithm
+ using Out_simplex_map = typename Manifold_tracing<Coxeter_triangulation<> >::Out_simplex_map;
+
+ std::vector<Eigen::VectorXd> seed_points(1, oracle.seed());
+ Out_simplex_map interior_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map);
+
+ // Constructing the cell complex
+ std::size_t intr_d = oracle.amb_d() - oracle.cod_d();
+ Cell_complex<Out_simplex_map> cell_complex(intr_d);
+ cell_complex.construct_complex(interior_simplex_map);
+
+ // List of Hasse_cell pointers to retrieve vertices values from edges
+ std::map<Cell_complex<Out_simplex_map>::Hasse_cell*, std::size_t> vi_map;
+ std::size_t index = 0;
+
+ std::clog << "Vertices:" << std::endl;
+ for (const auto& cp_pair : cell_complex.cell_point_map()) {
+ std::clog << index << " : (" << cp_pair.second(0) << ", " << cp_pair.second(1) << ")" << std::endl;
+ vi_map.emplace(cp_pair.first, index++);
+ }
+
+ std::clog << "Edges:" << std::endl;
+ for (const auto& sc_pair : cell_complex.interior_simplex_cell_map(1)) {
+ Cell_complex<Out_simplex_map>::Hasse_cell* edge_cell = sc_pair.second;
+ for (const auto& vi_pair : edge_cell->get_boundary()) std::clog << vi_map[vi_pair.first] << " ";
+ std::clog << std::endl;
+ }
+ return 0;
+}
diff --git a/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt
new file mode 100644
index 00000000..b323cca3
--- /dev/null
+++ b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt
@@ -0,0 +1,26 @@
+Vertices:
+0 : (-0.680375, 0.523483)
+1 : (0.147642, 0.887879)
+2 : (-0.847996, 0.30801)
+3 : (-0.881369, 0.0951903)
+4 : (0.638494, -0.550215)
+5 : (0.415344, 0.843848)
+6 : (0.812453, -0.0815816)
+7 : (0.319625, -0.7709)
+8 : (0.319625, 0.889605)
+9 : (0.579487, 0.638553)
+10 : (-0.680375, -0.461325)
+11 : (-0.364269, -0.760962)
+Edges:
+3 2
+3 10
+10 11
+11 7
+7 4
+2 0
+0 1
+6 9
+6 4
+1 8
+8 5
+5 9
diff --git a/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp b/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp
new file mode 100644
index 00000000..fe2051bb
--- /dev/null
+++ b/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp
@@ -0,0 +1,87 @@
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+#include <gudhi/Manifold_tracing.h>
+#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h>
+#include <gudhi/Functions/Linear_transformation.h>
+
+#include <gudhi/IO/build_mesh_from_cell_complex.h>
+#include <gudhi/IO/output_meshes_to_medit.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+/* A definition of a function that defines a 2d surface embedded in R^4, but that normally
+ * lives on a complex projective plane.
+ * In terms of harmonic coordinates [x:y:z] of points on the complex projective plane,
+ * the equation of the manifold is x^3*y + y^3*z + z^3*x = 0.
+ * The embedding consists of restricting the manifold to the affine subspace z = 1.
+ */
+struct Function_surface_on_CP2_in_R4 {
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ // The real and imaginary parts of the variables x and y
+ double xr = p(0), xi = p(1), yr = p(2), yi = p(3);
+ Eigen::VectorXd result(cod_d());
+
+ // Squares and cubes of real and imaginary parts used in the computations
+ double xr2 = xr * xr, xi2 = xi * xi, yr2 = yr * yr, yi2 = yi * yi, xr3 = xr2 * xr, xi3 = xi2 * xi, yr3 = yr2 * yr,
+ yi3 = yi2 * yi;
+
+ // The first coordinate of the output is Re(x^3*y + y^3 + x)
+ result(0) = xr3 * yr - 3 * xr * xi2 * yr - 3 * xr2 * xi * yi + xi3 * yi + yr3 - 3 * yr * yi2 + xr;
+ // The second coordinate of the output is Im(x^3*y + y^3 + x)
+ result(1) = 3 * xr2 * xi * yr + xr3 * yi - 3 * xr * xi2 * yi - xi3 * yr + 3 * yr2 * yi - yi3 + xi;
+ return result;
+ }
+
+ std::size_t amb_d() const { return 4; };
+ std::size_t cod_d() const { return 2; };
+
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = Eigen::VectorXd::Zero(4);
+ return result;
+ }
+
+ Function_surface_on_CP2_in_R4() {}
+};
+
+int main(int argc, char** argv) {
+ // The function for the (non-compact) manifold
+ Function_surface_on_CP2_in_R4 fun;
+
+ // Seed of the function
+ Eigen::VectorXd seed = fun.seed();
+
+ // Creating the function that defines the boundary of a compact region on the manifold
+ double radius = 3.0;
+ Function_Sm_in_Rd fun_sph(radius, 3, seed);
+
+ // Defining the intersection oracle
+ auto oracle = make_oracle(fun, fun_sph);
+
+ // Define a Coxeter triangulation scaled by a factor lambda.
+ // The triangulation is translated by a random vector to avoid violating the genericity hypothesis.
+ double lambda = 0.2;
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+ cox_tr.change_matrix(lambda * cox_tr.matrix());
+
+ // Manifold tracing algorithm
+ using MT = Manifold_tracing<Coxeter_triangulation<> >;
+ using Out_simplex_map = typename MT::Out_simplex_map;
+ std::vector<Eigen::VectorXd> seed_points(1, seed);
+ Out_simplex_map interior_simplex_map, boundary_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map, boundary_simplex_map);
+
+ // Constructing the cell complex
+ std::size_t intr_d = oracle.amb_d() - oracle.cod_d();
+ Cell_complex<Out_simplex_map> cell_complex(intr_d);
+ cell_complex.construct_complex(interior_simplex_map, boundary_simplex_map);
+
+ // Output the cell complex to a file readable by medit
+ output_meshes_to_medit(3, "manifold_on_CP2_with_boundary",
+ build_mesh_from_cell_complex(cell_complex, Configuration(true, true, true, 1, 5, 3),
+ Configuration(true, true, true, 2, 13, 14)));
+ return 0;
+}
diff --git a/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp b/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp
new file mode 100644
index 00000000..59fe2e2b
--- /dev/null
+++ b/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp
@@ -0,0 +1,72 @@
+// workaround for the annoying boost message in boost 1.69
+#define BOOST_PENDING_INTEGER_LOG2_HPP
+#include <boost/integer/integer_log2.hpp>
+// end workaround
+
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Functions/Function_affine_plane_in_Rd.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Cartesian_product.h>
+#include <gudhi/Functions/Linear_transformation.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+#include <gudhi/Manifold_tracing.h>
+#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h>
+#include <gudhi/Functions/random_orthogonal_matrix.h> // requires CGAL
+
+#include <gudhi/IO/build_mesh_from_cell_complex.h>
+#include <gudhi/IO/output_meshes_to_medit.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+int main(int argc, char** argv) {
+ // Creating a circle S1 in R2 of specified radius
+ double radius = 1.0;
+ Function_Sm_in_Rd fun_circle(radius, 1);
+
+ // Creating a flat torus S1xS1 in R4 from two circle functions
+ auto fun_flat_torus = make_product_function(fun_circle, fun_circle);
+
+ // Apply a random rotation in R4
+ auto matrix = random_orthogonal_matrix(4);
+ auto fun_flat_torus_rotated = make_linear_transformation(fun_flat_torus, matrix);
+
+ // Computing the seed of the function fun_flat_torus
+ Eigen::VectorXd seed = fun_flat_torus_rotated.seed();
+
+ // Defining a domain function that defines the boundary, which is a hyperplane passing by the origin and orthogonal to
+ // x.
+ Eigen::MatrixXd normal_matrix = Eigen::MatrixXd::Zero(4, 1);
+ for (std::size_t i = 0; i < 4; i++) normal_matrix(i, 0) = -seed(i);
+ Function_affine_plane_in_Rd fun_bound(normal_matrix, -seed / 2);
+
+ // Defining the intersection oracle
+ auto oracle = make_oracle(fun_flat_torus_rotated, fun_bound);
+
+ // Define a Coxeter triangulation scaled by a factor lambda.
+ // The triangulation is translated by a random vector to avoid violating the genericity hypothesis.
+ double lambda = 0.2;
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+ cox_tr.change_matrix(lambda * cox_tr.matrix());
+
+ // Manifold tracing algorithm
+ using MT = Manifold_tracing<Coxeter_triangulation<> >;
+ using Out_simplex_map = typename MT::Out_simplex_map;
+ std::vector<Eigen::VectorXd> seed_points(1, seed);
+ Out_simplex_map interior_simplex_map, boundary_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map, boundary_simplex_map);
+
+ // Constructing the cell complex
+ std::size_t intr_d = oracle.amb_d() - oracle.cod_d();
+ Cell_complex<Out_simplex_map> cell_complex(intr_d);
+ cell_complex.construct_complex(interior_simplex_map, boundary_simplex_map);
+
+ // Output the cell complex to a file readable by medit
+ output_meshes_to_medit(3, "flat_torus_with_boundary",
+ build_mesh_from_cell_complex(cell_complex, Configuration(true, true, true, 1, 5, 3),
+ Configuration(true, true, true, 2, 13, 14)));
+
+ return 0;
+}
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h
new file mode 100644
index 00000000..de68acb6
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h
@@ -0,0 +1,77 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef COXETER_TRIANGULATION_H_
+#define COXETER_TRIANGULATION_H_
+
+#include <vector>
+#include <cmath> // for std::sqrt
+
+#include <boost/range/iterator_range.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+
+#include <Eigen/Eigenvalues>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+
+#include <gudhi/Freudenthal_triangulation.h>
+#include <gudhi/Permutahedral_representation.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Coxeter_triangulation
+ * \brief A class that stores Coxeter triangulation of type \f$\tilde{A}_d\f$.
+ * This triangulation has the greatest simplex quality out of all linear transformations
+ * of the Freudenthal-Kuhn triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation.
+ * Needs to be a model of SimplexInCoxeterTriangulation.
+ */
+template <class Permutahedral_representation_ =
+ Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
+class Coxeter_triangulation : public Freudenthal_triangulation<Permutahedral_representation_> {
+ using Matrix = Eigen::MatrixXd;
+
+ Matrix root_matrix(unsigned d) {
+ Matrix cartan(Matrix::Identity(d, d));
+ for (unsigned i = 1; i < d; i++) {
+ cartan(i - 1, i) = -0.5;
+ cartan(i, i - 1) = -0.5;
+ }
+ Eigen::SelfAdjointEigenSolver<Matrix> saes(cartan);
+ Eigen::VectorXd sqrt_diag(d);
+ for (unsigned i = 0; i < d; ++i) sqrt_diag(i) = std::sqrt(saes.eigenvalues()[i]);
+
+ Matrix lower(Matrix::Ones(d, d));
+ lower = lower.triangularView<Eigen::Lower>();
+
+ Matrix result = (lower * saes.eigenvectors() * sqrt_diag.asDiagonal()).inverse();
+ return result;
+ }
+
+ public:
+ /** \brief Constructor of Coxeter triangulation of a given dimension.
+ * @param[in] dimension The dimension of the triangulation.
+ */
+ Coxeter_triangulation(std::size_t dimension)
+ : Freudenthal_triangulation<Permutahedral_representation_>(dimension, root_matrix(dimension)) {}
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h
new file mode 100644
index 00000000..de342ecc
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h
@@ -0,0 +1,340 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CELL_COMPLEX_H_
+#define CELL_COMPLEX_H_
+
+#include <Eigen/Dense>
+
+#include <vector>
+#include <map>
+#include <utility> // for std::make_pair
+
+#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES
+#include <gudhi/Permutahedral_representation/Simplex_comparator.h>
+#include <gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h> // for Hasse_cell
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Cell_complex
+ * \brief A class that constructs the cell complex from the output provided by the class
+ * \ref Gudhi::coxeter_triangulation::Manifold_tracing.
+ *
+ * The use and interfaces of this cell complex is limited to the \ref coxeter_triangulation implementation.
+ *
+ * \tparam Out_simplex_map_ The type of a map from a simplex type that is a
+ * model of SimplexInCoxeterTriangulation to Eigen::VectorXd.
+ */
+template <class Out_simplex_map_>
+class Cell_complex {
+ public:
+ /** \brief Type of a simplex in the ambient triangulation.
+ * Is a model of the concept SimplexInCoxeterTriangulation.
+ */
+ using Simplex_handle = typename Out_simplex_map_::key_type;
+ /** \brief Type of a cell in the cell complex.
+ * Always is Gudhi::Hasse_cell from the Hasse diagram module.
+ * The additional information is the boolean that is true if and only if the cell lies
+ * on the boundary.
+ */
+ using Hasse_cell = Gudhi::Hasse_diagram::Hasse_diagram_cell<int, double, bool>;
+ /** \brief Type of a map from permutahedral representations of simplices in the
+ * ambient triangulation to the corresponding cells in the cell complex of some
+ * specific dimension.
+ */
+ using Simplex_cell_map = std::map<Simplex_handle, Hasse_cell*, Simplex_comparator<Simplex_handle> >;
+ /** \brief Type of a vector of maps from permutahedral representations of simplices in the
+ * ambient triangulation to the corresponding cells in the cell complex of various dimensions.
+ */
+ using Simplex_cell_maps = std::vector<Simplex_cell_map>;
+
+ /** \brief Type of a map from cells in the cell complex to the permutahedral representations
+ * of the corresponding simplices in the ambient triangulation.
+ */
+ using Cell_simplex_map = std::map<Hasse_cell*, Simplex_handle>;
+
+ /** \brief Type of a map from vertex cells in the cell complex to the permutahedral representations
+ * of their Cartesian coordinates.
+ */
+ using Cell_point_map = std::map<Hasse_cell*, Eigen::VectorXd>;
+
+ private:
+ Hasse_cell* insert_cell(const Simplex_handle& simplex, std::size_t cell_d, bool is_boundary) {
+ Simplex_cell_maps& simplex_cell_maps = (is_boundary ? boundary_simplex_cell_maps_ : interior_simplex_cell_maps_);
+#ifdef DEBUG_TRACES
+ CC_detail_list& cc_detail_list =
+ (is_boundary ? cc_boundary_detail_lists[cell_d] : cc_interior_detail_lists[cell_d]);
+ cc_detail_list.emplace_back(simplex);
+#endif
+ Simplex_cell_map& simplex_cell_map = simplex_cell_maps[cell_d];
+ auto map_it = simplex_cell_map.find(simplex);
+ if (map_it == simplex_cell_map.end()) {
+ hasse_cells_.push_back(new Hasse_cell(is_boundary, cell_d));
+ Hasse_cell* new_cell = hasse_cells_.back();
+ simplex_cell_map.emplace(simplex, new_cell);
+ cell_simplex_map_.emplace(new_cell, simplex);
+#ifdef DEBUG_TRACES
+ cc_detail_list.back().status_ = CC_detail_info::Result_type::inserted;
+#endif
+ return new_cell;
+ }
+#ifdef DEBUG_TRACES
+ CC_detail_info& cc_info = cc_detail_list.back();
+ cc_info.trigger_ = to_string(map_it->first);
+ cc_info.status_ = CC_detail_info::Result_type::self;
+#endif
+ return map_it->second;
+ }
+
+ void expand_level(std::size_t cell_d) {
+ bool is_manifold_with_boundary = boundary_simplex_cell_maps_.size() > 0;
+ for (auto& sc_pair : interior_simplex_cell_maps_[cell_d - 1]) {
+ const Simplex_handle& simplex = sc_pair.first;
+ Hasse_cell* cell = sc_pair.second;
+ for (Simplex_handle coface : simplex.coface_range(cod_d_ + cell_d)) {
+ Hasse_cell* new_cell = insert_cell(coface, cell_d, false);
+ new_cell->get_boundary().emplace_back(cell, 1);
+ }
+ }
+
+ if (is_manifold_with_boundary) {
+ for (auto& sc_pair : boundary_simplex_cell_maps_[cell_d - 1]) {
+ const Simplex_handle& simplex = sc_pair.first;
+ Hasse_cell* cell = sc_pair.second;
+ if (cell_d != intr_d_)
+ for (Simplex_handle coface : simplex.coface_range(cod_d_ + cell_d + 1)) {
+ Hasse_cell* new_cell = insert_cell(coface, cell_d, true);
+ new_cell->get_boundary().emplace_back(cell, 1);
+ }
+ auto map_it = interior_simplex_cell_maps_[cell_d].find(simplex);
+ if (map_it == interior_simplex_cell_maps_[cell_d].end())
+ std::cerr << "Cell_complex::expand_level error: A boundary cell does not have an interior counterpart.\n";
+ else {
+ Hasse_cell* i_cell = map_it->second;
+ i_cell->get_boundary().emplace_back(cell, 1);
+ }
+ }
+ }
+ }
+
+ void construct_complex_(const Out_simplex_map_& out_simplex_map) {
+#ifdef DEBUG_TRACES
+ cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size());
+ cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size());
+ cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size());
+#endif
+ for (auto& os_pair : out_simplex_map) {
+ const Simplex_handle& simplex = os_pair.first;
+ const Eigen::VectorXd& point = os_pair.second;
+ Hasse_cell* new_cell = insert_cell(simplex, 0, false);
+ cell_point_map_.emplace(new_cell, point);
+ }
+ for (std::size_t cell_d = 1;
+ cell_d < interior_simplex_cell_maps_.size() && !interior_simplex_cell_maps_[cell_d - 1].empty(); ++cell_d) {
+ expand_level(cell_d);
+ }
+ }
+
+ void construct_complex_(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map) {
+#ifdef DEBUG_TRACES
+ cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size());
+ cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size());
+ cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size());
+ cc_boundary_summary_lists.resize(boundary_simplex_cell_maps_.size());
+ cc_boundary_prejoin_lists.resize(boundary_simplex_cell_maps_.size());
+ cc_boundary_detail_lists.resize(boundary_simplex_cell_maps_.size());
+#endif
+ for (auto& os_pair : boundary_simplex_map) {
+ const Simplex_handle& simplex = os_pair.first;
+ const Eigen::VectorXd& point = os_pair.second;
+ Hasse_cell* new_cell = insert_cell(simplex, 0, true);
+ cell_point_map_.emplace(new_cell, point);
+ }
+ for (auto& os_pair : interior_simplex_map) {
+ const Simplex_handle& simplex = os_pair.first;
+ const Eigen::VectorXd& point = os_pair.second;
+ Hasse_cell* new_cell = insert_cell(simplex, 0, false);
+ cell_point_map_.emplace(new_cell, point);
+ }
+#ifdef DEBUG_TRACES
+ for (const auto& sc_pair : interior_simplex_cell_maps_[0])
+ cc_interior_summary_lists[0].push_back(CC_summary_info(sc_pair));
+ for (const auto& sc_pair : boundary_simplex_cell_maps_[0])
+ cc_boundary_summary_lists[0].push_back(CC_summary_info(sc_pair));
+#endif
+
+ for (std::size_t cell_d = 1;
+ cell_d < interior_simplex_cell_maps_.size() && !interior_simplex_cell_maps_[cell_d - 1].empty(); ++cell_d) {
+ expand_level(cell_d);
+
+#ifdef DEBUG_TRACES
+ for (const auto& sc_pair : interior_simplex_cell_maps_[cell_d])
+ cc_interior_summary_lists[cell_d].push_back(CC_summary_info(sc_pair));
+ if (cell_d < boundary_simplex_cell_maps_.size())
+ for (const auto& sc_pair : boundary_simplex_cell_maps_[cell_d])
+ cc_boundary_summary_lists[cell_d].push_back(CC_summary_info(sc_pair));
+#endif
+ }
+ }
+
+ public:
+ /**
+ * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold
+ * without boundary embedded in the \f$ d \f$-dimensional Euclidean space
+ * from the output of the class Gudhi::Manifold_tracing.
+ *
+ * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$
+ * in the ambient triangulation that intersect the relative interior of the manifold
+ * to the intersection points.
+ */
+ void construct_complex(const Out_simplex_map_& out_simplex_map) {
+ interior_simplex_cell_maps_.resize(intr_d_ + 1);
+ if (!out_simplex_map.empty()) cod_d_ = out_simplex_map.begin()->first.dimension();
+ construct_complex_(out_simplex_map);
+ }
+
+ /**
+ * \brief Constructs the skeleton of the cell complex that approximates
+ * an \f$m\f$-dimensional manifold without boundary embedded
+ * in the \f$d\f$-dimensional Euclidean space
+ * up to a limit dimension from the output of the class Gudhi::Manifold_tracing.
+ *
+ * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$
+ * in the ambient triangulation that intersect the relative interior of the manifold
+ * to the intersection points.
+ * \param[in] limit_dimension The dimension of the constructed skeleton.
+ */
+ void construct_complex(const Out_simplex_map_& out_simplex_map, std::size_t limit_dimension) {
+ interior_simplex_cell_maps_.resize(limit_dimension + 1);
+ if (!out_simplex_map.empty()) cod_d_ = out_simplex_map.begin()->first.dimension();
+ construct_complex_(out_simplex_map);
+ }
+
+ /**
+ * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold
+ * with boundary embedded in the \f$ d \f$-dimensional Euclidean space
+ * from the output of the class Gudhi::Manifold_tracing.
+ *
+ * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$
+ * in the ambient triangulation that intersect the relative interior of the manifold
+ * to the intersection points.
+ * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$
+ * in the ambient triangulation that intersect the boundary of the manifold
+ * to the intersection points.
+ */
+ void construct_complex(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map) {
+ interior_simplex_cell_maps_.resize(intr_d_ + 1);
+ boundary_simplex_cell_maps_.resize(intr_d_);
+ if (!interior_simplex_map.empty()) cod_d_ = interior_simplex_map.begin()->first.dimension();
+ construct_complex_(interior_simplex_map, boundary_simplex_map);
+ }
+
+ /**
+ * \brief Constructs the skeleton of the cell complex that approximates
+ * an \f$m\f$-dimensional manifold with boundary embedded
+ * in the \f$d\f$-dimensional Euclidean space
+ * up to a limit dimension from the output of the class Gudhi::Manifold_tracing.
+ *
+ * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$
+ * in the ambient triangulation that intersect the relative interior of the manifold
+ * to the intersection points.
+ * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$
+ * in the ambient triangulation that intersect the boundary of the manifold
+ * to the intersection points.
+ * \param[in] limit_dimension The dimension of the constructed skeleton.
+ */
+ void construct_complex(const Out_simplex_map_& interior_simplex_map, const Out_simplex_map_& boundary_simplex_map,
+ std::size_t limit_dimension) {
+ interior_simplex_cell_maps_.resize(limit_dimension + 1);
+ boundary_simplex_cell_maps_.resize(limit_dimension);
+ if (!interior_simplex_map.empty()) cod_d_ = interior_simplex_map.begin()->first.dimension();
+ construct_complex_(interior_simplex_map, boundary_simplex_map);
+ }
+
+ /**
+ * \brief Returns the dimension of the cell complex.
+ */
+ std::size_t intrinsic_dimension() const { return intr_d_; }
+
+ /**
+ * \brief Returns a vector of maps from the cells of various dimensions in the interior
+ * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations
+ * of the corresponding simplices in the ambient triangulation.
+ */
+ const Simplex_cell_maps& interior_simplex_cell_maps() const { return interior_simplex_cell_maps_; }
+
+ /**
+ * \brief Returns a vector of maps from the cells of various dimensions on the boundary
+ * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations
+ * of the corresponding simplices in the ambient triangulation.
+ */
+ const Simplex_cell_maps& boundary_simplex_cell_maps() const { return boundary_simplex_cell_maps_; }
+
+ /**
+ * \brief Returns a map from the cells of a given dimension in the interior
+ * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations
+ * of the corresponding simplices in the ambient triangulation.
+ *
+ * \param[in] cell_d The dimension of the cells.
+ */
+ const Simplex_cell_map& interior_simplex_cell_map(std::size_t cell_d) const {
+ return interior_simplex_cell_maps_[cell_d];
+ }
+
+ /**
+ * \brief Returns a map from the cells of a given dimension on the boundary
+ * of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations
+ * of the corresponding simplices in the ambient triangulation.
+ *
+ * \param[in] cell_d The dimension of the cells.
+ */
+ const Simplex_cell_map& boundary_simplex_cell_map(std::size_t cell_d) const {
+ return boundary_simplex_cell_maps_[cell_d];
+ }
+
+ /**
+ * \brief Returns a map from the cells in the cell complex of type Gudhi::Hasse_cell
+ * to the permutahedral representations of the corresponding simplices in the
+ * ambient triangulation.
+ */
+ const Cell_simplex_map& cell_simplex_map() const { return cell_simplex_map_; }
+
+ /**
+ * \brief Returns a map from the vertex cells in the cell complex of type Gudhi::Hasse_cell
+ * to their Cartesian coordinates.
+ */
+ const Cell_point_map& cell_point_map() const { return cell_point_map_; }
+
+ /**
+ * \brief Constructor for the class Cell_complex.
+ *
+ * \param[in] intrinsic_dimension The dimension of the cell complex.
+ */
+ Cell_complex(std::size_t intrinsic_dimension) : intr_d_(intrinsic_dimension) {}
+
+ ~Cell_complex() {
+ for (Hasse_cell* hs_ptr : hasse_cells_) delete hs_ptr;
+ }
+
+ private:
+ std::size_t intr_d_, cod_d_;
+ Simplex_cell_maps interior_simplex_cell_maps_, boundary_simplex_cell_maps_;
+ Cell_simplex_map cell_simplex_map_;
+ Cell_point_map cell_point_map_;
+ std::vector<Hasse_cell*> hasse_cells_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h
new file mode 100644
index 00000000..59e9a350
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Cell_complex/Hasse_diagram_cell.h
@@ -0,0 +1,285 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Pawel Dlotko
+ *
+ * Copyright (C) 2017 Swansea University UK
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef HASSE_DIAGRAM_CELL_H
+#define HASSE_DIAGRAM_CELL_H
+
+#include <vector>
+#include <utility> // for std::pair
+#include <ostream>
+#include <string>
+#include <type_traits> // for std::is_same
+#include <cstdlib> // for std::size_t
+
+namespace Gudhi {
+namespace Hasse_diagram {
+
+template <typename Cell_type>
+class Hasse_diagram;
+
+/**
+ * \class Hasse_diagram_cell
+ * \brief Data structure to store a cell in a Hasse diagram.
+ *
+ * \ingroup Hasse_diagram
+ *
+ * \details
+ * The use and interfaces of this Hasse diagram cell is limited to the \ref coxeter_triangulation implementation.
+ *
+ * This is a data structure to store a cell in a general Hasse diagram data structure. It stores the following
+ * information about the cell: References to boundary and coBoundary elements, dimension of a cell and its filtration.
+ * It also allow to store any additional information of a type Additional_information which is a template parameter of
+ * the class (set by default to void).
+ *
+ * The complex is a template class requiring the following parameters:
+ * Incidence_type_ - determine the type of incidence coefficients. Use integers in most general case.
+ * Filtration_type_ - type of filtration of cells.
+ * Additional_information_ (set by default to void) - allows to store any
+ * additional information in the cells of Hasse diagrams.
+ *
+ */
+template <typename Incidence_type_, typename Filtration_type_, typename Additional_information_ = void>
+class Hasse_diagram_cell {
+ public:
+ typedef Incidence_type_ Incidence_type;
+ typedef Filtration_type_ Filtration_type;
+ typedef Additional_information_ Additional_information;
+ using Cell_range = std::vector<std::pair<Hasse_diagram_cell*, Incidence_type> >;
+
+ /**
+ * Default constructor.
+ **/
+ Hasse_diagram_cell() : dimension(0), position(0), deleted_(false) {}
+
+ /**
+ * Constructor of a cell of dimension dim.
+ **/
+ Hasse_diagram_cell(int dim) : dimension(dim), position(0), deleted_(false) {}
+
+ /**
+ * Constructor of a cell of dimension dim.
+ **/
+ Hasse_diagram_cell(int dim, Filtration_type filt_)
+ : dimension(dim), position(0), deleted_(false), filtration(filt_) {}
+
+ /**
+ * Constructor of a cell of dimension dim with a given boundary.
+ **/
+ Hasse_diagram_cell(const Cell_range& boundary_, int dim)
+ : dimension(dim), boundary(boundary_), position(0), deleted_(false) {}
+
+ /**
+ * Constructor of a cell of dimension dim with a given boundary and coboundary.
+ **/
+ Hasse_diagram_cell(const Cell_range& boundary_, const Cell_range& coboundary_, int dim)
+ : dimension(dim), boundary(boundary_), coBoundary(coboundary_), position(0), deleted_(false) {}
+
+ /**
+ * Constructor of a cell of dimension dim with a given boundary, coboundary and
+ * additional information.
+ **/
+ Hasse_diagram_cell(const Cell_range& boundary_, const Cell_range& coboundary_, const Additional_information& ai,
+ int dim)
+ : dimension(dim),
+ boundary(boundary_),
+ coBoundary(coboundary_),
+ additional_info(ai),
+ position(0),
+ deleted_(false) {}
+
+ /**
+ * 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 <typename Complex_type, typename Cell_type>
+ friend std::vector<Cell_type*> convert_to_vector_of_Cell_type(Complex_type& cmplx);
+
+ /**
+ * Procedure to remove deleted boundary and coboundary elements from the
+ * vectors of boundary and coboundary elements of this cell.
+ **/
+ void remove_deleted_elements_from_boundary_and_coboundary() {
+ Cell_range new_boundary;
+ new_boundary.reserve(this->boundary.size());
+ for (std::size_t bd = 0; bd != this->boundary.size(); ++bd) {
+ if (!this->boundary[bd].first->deleted()) {
+ new_boundary.push_back(this->boundary[bd]);
+ }
+ }
+ this->boundary.swap(new_boundary);
+
+ Cell_range new_coBoundary;
+ new_coBoundary.reserve(this->coBoundary.size());
+ for (std::size_t cbd = 0; cbd != this->coBoundary.size(); ++cbd) {
+ if (!this->coBoundary[cbd].first->deleted()) {
+ new_coBoundary.push_back(this->coBoundary[cbd]);
+ }
+ }
+ this->coBoundary.swap(new_coBoundary);
+ }
+
+ /**
+ * Writing to a stream operator.
+ **/
+ friend std::ostream& operator<<(
+ std::ostream& out, const Hasse_diagram_cell<Incidence_type, Filtration_type, Additional_information>& c) {
+ // cout << "position : " << c.position << ", dimension : " << c.dimension << ", filtration: " << c.filtration << ",
+ // size of boudary : " << c.boundary.size() << "\n";
+ out << c.position << " " << c.dimension << " " << c.filtration << std::endl;
+ for (std::size_t bd = 0; bd != c.boundary.size(); ++bd) {
+ // do not write out the cells that has been deleted
+ if (c.boundary[bd].first->deleted()) continue;
+ out << c.boundary[bd].first->position << " " << c.boundary[bd].second << " ";
+ }
+ out << std::endl;
+ return out;
+ }
+
+ /**
+ * Procedure that return vector of pointers to boundary elements of a given cell.
+ **/
+ inline std::vector<Hasse_diagram_cell*> get_list_of_boundary_elements() {
+ std::vector<Hasse_diagram_cell*> result;
+ std::size_t size_of_boundary = this->boundary.size();
+ result.reserve(size_of_boundary);
+ for (std::size_t bd = 0; bd != size_of_boundary; ++bd) {
+ result.push_back(this->boundary[bd].first);
+ }
+ return result;
+ }
+
+ /**
+ * Procedure that return vector of positios of boundary elements of a given cell.
+ **/
+ inline std::vector<unsigned> get_list_of_positions_of_boundary_elements() {
+ std::vector<unsigned> result;
+ std::size_t size_of_boundary = this->boundary.size();
+ result.reserve(size_of_boundary);
+ for (std::size_t bd = 0; bd != size_of_boundary; ++bd) {
+ result.push_back(this->boundary[bd].first->position);
+ }
+ return result;
+ }
+
+ /**
+ * Function that display a string being a signature of a structure.
+ * Used mainly for debugging purposes.
+ **/
+ std::string full_signature_of_the_structure() {
+ std::string result;
+ result += "dimension: ";
+ result += std::to_string(this->dimension);
+ result += " filtration: ";
+ result += std::to_string(this->filtration);
+ result += " position: ";
+ result += std::to_string(this->position);
+ result += " deleted_: ";
+ result += std::to_string(this->deleted_);
+
+ // if the Additional_information is not void, add them to
+ // the signature as well.
+ if (std::is_same<Additional_information, void>::value) {
+ result += " Additional_information: ";
+ result += std::to_string(this->additional_info);
+ }
+ result += " boundary ";
+ for (std::size_t bd = 0; bd != this->boundary.size(); ++bd) {
+ result += "( " + std::to_string(this->boundary[bd].first->position);
+ result += " " + std::to_string(this->boundary[bd].second);
+ result += ") ";
+ }
+
+ result += " coBoundary ";
+ for (std::size_t cbd = 0; cbd != this->coBoundary.size(); ++cbd) {
+ result += "( " + std::to_string(this->coBoundary[cbd].first->position);
+ result += " " + std::to_string(this->coBoundary[cbd].second);
+ result += ") ";
+ }
+
+ return result;
+ }
+
+ protected:
+ Cell_range boundary;
+ Cell_range coBoundary;
+ int dimension;
+ Additional_information additional_info;
+ unsigned position;
+ bool deleted_;
+ Filtration_type filtration;
+
+ /**
+ * A procedure to delete a cell. It is a private function of the Hasse_diagram_cell
+ * class, since in the Hasse_diagram class I want to have a control
+ * of removal of cells. Therefore, to remove cell please use
+ * remove_cell in the Hasse_diagram structure.
+ **/
+ void delete_cell() { this->deleted_ = true; }
+}; // Hasse_diagram_cell
+
+} // namespace Hasse_diagram
+} // namespace Gudhi
+
+#endif // CELL_H
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h
new file mode 100644
index 00000000..5543c2fb
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation/Query_result.h
@@ -0,0 +1,40 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef QUERY_RESULT_H_
+#define QUERY_RESULT_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Query_result
+ * \brief The result of a query by an oracle such as Implicit_manifold_intersection_oracle.
+ *
+ * \tparam Simplex_handle The class of the query simplex.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Simplex_handle>
+struct Query_result {
+ /** \brief The potentially lower-dimensional face of the query simplex
+ * that contains the intersection point. OBSOLETE: as the snapping is removed. */
+ // Simplex_handle face;
+ /** \brief The intersection point. */
+ Eigen::VectorXd intersection;
+ /** \brief True if the query simplex intersects the manifold. */
+ bool success;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h
new file mode 100644
index 00000000..873c5c9b
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h
@@ -0,0 +1,219 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FREUDENTHAL_TRIANGULATION_H_
+#define FREUDENTHAL_TRIANGULATION_H_
+
+#include <vector>
+#include <algorithm> // for std::sort
+#include <cmath> // for std::floor
+#include <numeric> // for std::iota
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Eigenvalues>
+#include <Eigen/SVD>
+
+#include <gudhi/Permutahedral_representation.h>
+#include <gudhi/Debug_utils.h> // for GUDHI_CHECK
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Freudenthal_triangulation
+ * \brief A class that stores any affine transformation of the Freudenthal-Kuhn
+ * triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \details The data structure is a record that consists of a matrix
+ * that represents the linear transformation of the Freudenthal-Kuhn triangulation
+ * and a vector that represents the offset.
+ *
+ * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation.
+ * Needs to be a model of SimplexInCoxeterTriangulation.
+ */
+template <class Permutahedral_representation_ =
+ Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
+class Freudenthal_triangulation {
+ using Matrix = Eigen::MatrixXd;
+ using Vector = Eigen::VectorXd;
+
+ public:
+ /** \brief Type of the simplices in the triangulation. */
+ using Simplex_handle = Permutahedral_representation_;
+
+ /** \brief Type of the vertices in the triangulation. */
+ using Vertex_handle = typename Permutahedral_representation_::Vertex;
+
+ /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension.
+ * @param[in] dimension The dimension of the triangulation.
+ */
+ Freudenthal_triangulation(std::size_t dimension)
+ : Freudenthal_triangulation(dimension, Matrix::Identity(dimension, dimension), Vector::Zero(dimension)) {
+ is_freudenthal_ = true;
+ }
+
+ /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under
+ * a linear transformation by a given matrix.
+ * @param[in] dimension The dimension of the triangulation.
+ * @param[in] matrix The matrix that defines the linear transformation.
+ * Needs to be invertible.
+ */
+ Freudenthal_triangulation(std::size_t dimension, const Matrix& matrix)
+ : Freudenthal_triangulation(dimension, matrix, Vector::Zero(dimension)) {}
+
+ /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under
+ * an affine transformation by a given matrix and a translation vector.
+ * @param[in] dimension The dimension of the triangulation.
+ * @param[in] matrix The matrix that defines the linear transformation.
+ * Needs to be invertible.
+ * @param[in] offset The offset vector.
+ *
+ * @exception std::invalid_argument In debug mode, if offset size is different from dimension.
+ */
+ Freudenthal_triangulation(unsigned dimension, const Matrix& matrix, const Vector& offset)
+ : dimension_(dimension),
+ matrix_(matrix),
+ offset_(offset),
+ colpivhouseholderqr_(matrix_.colPivHouseholderQr()),
+ is_freudenthal_(false) {
+ GUDHI_CHECK(dimension == offset_.size(), std::invalid_argument("Offset must be of size 'dimension'"));
+ }
+
+ /** \brief Dimension of the triangulation. */
+ unsigned dimension() const { return dimension_; }
+
+ /** \brief Matrix that defines the linear transformation of the triangulation. */
+ const Matrix& matrix() const { return matrix_; }
+
+ /** \brief Vector that defines the offset of the triangulation. */
+ const Vector& offset() const { return offset_; }
+
+ /** \brief Change the linear transformation matrix to a given value.
+ * @param[in] matrix New value of the linear transformation matrix.
+ */
+ void change_matrix(const Eigen::MatrixXd& matrix) {
+ matrix_ = matrix;
+ colpivhouseholderqr_ = matrix.colPivHouseholderQr();
+ is_freudenthal_ = false;
+ }
+
+ /** \brief Change the offset vector to a given value.
+ * @param[in] offset New value of the offset vector.
+ */
+ void change_offset(const Eigen::VectorXd& offset) {
+ offset_ = offset;
+ is_freudenthal_ = false;
+ }
+
+ /** \brief Returns the permutahedral representation of the simplex in the
+ * triangulation that contains a given query point.
+ * \details Using the additional parameter scale, the search can be done in a
+ * triangulation that shares the origin, but is scaled by a given factor.
+ * This parameter can be useful to simulate the point location in a subdivided
+ * triangulation.
+ * The returned simplex is always minimal by inclusion.
+ *
+ * \tparam Point_d A class that represents a point in d-dimensional Euclidean space.
+ * The coordinates should be random-accessible. Needs to provide the method size().
+ *
+ * @param[in] point The query point.
+ * @param[in] scale The scale of the triangulation.
+ *
+ * @exception std::invalid_argument In debug mode, if point dimension is different from triangulation one.
+ */
+ template <class Point_d>
+ Simplex_handle locate_point(const Point_d& point, double scale = 1) const {
+ using Ordered_set_partition = typename Simplex_handle::OrderedSetPartition;
+ using Part = typename Ordered_set_partition::value_type;
+ unsigned d = point.size();
+ GUDHI_CHECK(d == dimension_,
+ std::invalid_argument("The point must be of the same dimension as the triangulation"));
+ double error = 1e-9;
+ Simplex_handle output;
+ std::vector<double> z;
+ if (is_freudenthal_) {
+ for (std::size_t i = 0; i < d; i++) {
+ double x_i = scale * point[i];
+ int y_i = std::floor(x_i);
+ output.vertex().push_back(y_i);
+ z.push_back(x_i - y_i);
+ }
+ } else {
+ Eigen::VectorXd p_vect(d);
+ for (std::size_t i = 0; i < d; i++) p_vect(i) = point[i];
+ Eigen::VectorXd x_vect = colpivhouseholderqr_.solve(p_vect - offset_);
+ for (std::size_t i = 0; i < d; i++) {
+ double x_i = scale * x_vect(i);
+ int y_i = std::floor(x_i);
+ output.vertex().push_back(y_i);
+ z.push_back(x_i - y_i);
+ }
+ }
+ z.push_back(0);
+ Part indices(d + 1);
+ std::iota(indices.begin(), indices.end(), 0);
+ std::sort(indices.begin(), indices.end(), [&z](std::size_t i1, std::size_t i2) { return z[i1] > z[i2]; });
+
+ output.partition().push_back(Part(1, indices[0]));
+ for (std::size_t i = 1; i <= d; ++i)
+ if (z[indices[i - 1]] > z[indices[i]] + error)
+ output.partition().push_back(Part(1, indices[i]));
+ else
+ output.partition().back().push_back(indices[i]);
+ return output;
+ }
+
+ /** \brief Returns the Cartesian coordinates of the given vertex.
+ * \details Using the additional parameter scale, the search can be done in a
+ * triangulation that shares the origin, but is scaled by a given factor.
+ * This parameter can be useful to simulate the computation of Cartesian coordinates
+ * of a vertex in a subdivided triangulation.
+ * @param[in] vertex The query vertex.
+ * @param[in] scale The scale of the triangulation.
+ */
+ Eigen::VectorXd cartesian_coordinates(const Vertex_handle& vertex, double scale = 1) const {
+ Eigen::VectorXd v_vect(dimension_);
+ for (std::size_t j = 0; j < dimension_; j++) v_vect(j) = vertex[j] / scale;
+ return matrix_ * v_vect + offset_;
+ }
+
+ /** \brief Returns the Cartesian coordinates of the barycenter of a given simplex.
+ * \details Using the additional parameter scale, the search can be done in a
+ * triangulation that shares the origin, but is scaled by a given factor.
+ * This parameter can be useful to simulate the computation of Cartesian coordinates
+ * of the barycenter of a simplex in a subdivided triangulation.
+ * @param[in] simplex The query simplex.
+ * @param[in] scale The scale of the triangulation.
+ */
+ Eigen::VectorXd barycenter(const Simplex_handle& simplex, double scale = 1) const {
+ Eigen::VectorXd res_vector(dimension_);
+ res_vector.setZero(dimension_, 1);
+ for (auto v : simplex.vertex_range()) {
+ res_vector += cartesian_coordinates(v, scale);
+ }
+ return (1. / (simplex.dimension() + 1)) * res_vector;
+ }
+
+ protected:
+ unsigned dimension_;
+ Matrix matrix_;
+ Vector offset_;
+ Eigen::ColPivHouseholderQR<Matrix> colpivhouseholderqr_;
+ bool is_freudenthal_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h
new file mode 100644
index 00000000..0533bb83
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h
@@ -0,0 +1,157 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_CARTESIAN_PRODUCT_H_
+#define FUNCTIONS_CARTESIAN_PRODUCT_H_
+
+#include <cstdlib>
+#include <tuple>
+#include <type_traits> // for std::enable_if
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* Get the domain dimension of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof...(T), std::size_t>::type get_amb_d(const std::tuple<T...>& tuple) {
+ return 0;
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof...(T), std::size_t>::type get_amb_d(const std::tuple<T...>& tuple) {
+ return std::get<I>(tuple).amb_d() + get_amb_d<I + 1, T...>(tuple);
+}
+
+/* Get the codomain dimension of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof...(T), std::size_t>::type get_cod_d(const std::tuple<T...>& tuple) {
+ return 0;
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof...(T), std::size_t>::type get_cod_d(const std::tuple<T...>& tuple) {
+ return std::get<I>(tuple).cod_d() + get_cod_d<I + 1, T...>(tuple);
+}
+
+/* Get the seed of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof...(T), void>::type get_seed(const std::tuple<T...>& tuple,
+ Eigen::VectorXd& point, std::size_t i = 0) {}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof...(T), void>::type get_seed(const std::tuple<T...>& tuple,
+ Eigen::VectorXd& point, std::size_t i = 0) {
+ const auto& f = std::get<I>(tuple);
+ std::size_t n = f.amb_d();
+ Eigen::VectorXd seed = f.seed();
+ for (std::size_t j = 0; j < n; ++j) point(i + j) = seed(j);
+ get_seed<I + 1, T...>(tuple, point, i + n);
+}
+
+/* Get the seed of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof...(T), void>::type get_value(const std::tuple<T...>& tuple,
+ const Eigen::VectorXd& x,
+ Eigen::VectorXd& point, std::size_t i = 0,
+ std::size_t j = 0) {}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof...(T), void>::type get_value(const std::tuple<T...>& tuple,
+ const Eigen::VectorXd& x,
+ Eigen::VectorXd& point, std::size_t i = 0,
+ std::size_t j = 0) {
+ const auto& f = std::get<I>(tuple);
+ std::size_t n = f.amb_d();
+ std::size_t k = f.cod_d();
+ Eigen::VectorXd x_i(n);
+ for (std::size_t l = 0; l < n; ++l) x_i(l) = x(i + l);
+ Eigen::VectorXd res = f(x_i);
+ for (std::size_t l = 0; l < k; ++l) point(j + l) = res(l);
+ get_value<I + 1, T...>(tuple, x, point, i + n, j + k);
+}
+
+/**
+ * \class Cartesian_product
+ * \brief Constructs the function the zero-set of which is the Cartesian product
+ * of the zero-sets of some given functions.
+ *
+ * \tparam Functions A pack template parameter for functions. All functions should be models of
+ * the concept FunctionForImplicitManifold.
+ */
+template <class... Functions>
+struct Cartesian_product {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result(cod_d_);
+ get_value(function_tuple_, p, result, 0, 0);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return amb_d_; }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return cod_d_; }
+
+ /** \brief Returns a point on the zero-set. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result(amb_d_);
+ get_seed(function_tuple_, result, 0);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the Cartesian product function.
+ *
+ * @param[in] functions The functions the zero-sets of which are factors in the
+ * Cartesian product of the resulting function.
+ */
+ Cartesian_product(const Functions&... functions) : function_tuple_(std::make_tuple(functions...)) {
+ amb_d_ = get_amb_d(function_tuple_);
+ cod_d_ = get_cod_d(function_tuple_);
+ }
+
+ private:
+ std::tuple<Functions...> function_tuple_;
+ std::size_t amb_d_, cod_d_;
+};
+
+/**
+ * \brief Static constructor of a Cartesian product function.
+ *
+ * @param[in] functions The functions the zero-sets of which are factors in the
+ * Cartesian product of the resulting function.
+ *
+ * \tparam Functions A pack template parameter for functions. All functions should be models of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <typename... Functions>
+Cartesian_product<Functions...> make_product_function(const Functions&... functions) {
+ return Cartesian_product<Functions...>(functions...);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h
new file mode 100644
index 00000000..0603afd8
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h
@@ -0,0 +1,64 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_CONSTANT_FUNCTION_H_
+#define FUNCTIONS_CONSTANT_FUNCTION_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Constant_function
+ * \brief A class that encodes a constant function from R^d to R^k.
+ * This class does not have any implicit manifold in correspondence.
+ */
+struct Constant_function {
+ /** \brief Value of the function at a specified point. The value is constant.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ return value_;
+ }
+
+ /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+ std::size_t amb_d() const { return d_; };
+
+ /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+ std::size_t cod_d() const { return k_; };
+
+ /** \brief No seed point is available. Throws an exception on evocation. */
+ Eigen::VectorXd seed() const { throw "Seed invoked on a constant function.\n"; }
+
+ Constant_function() {}
+
+ /**
+ * \brief Constructor of a constant function from R^d to R^m.
+ *
+ * @param[in] d The domain dimension.
+ * @param[in] k The codomain dimension.
+ * @param[in] value The constant value of the function.
+ */
+ Constant_function(std::size_t d, std::size_t k, const Eigen::VectorXd& value) : d_(d), k_(k), value_(value) {}
+
+ private:
+ std::size_t d_, k_;
+ Eigen::VectorXd value_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h
new file mode 100644
index 00000000..e1fe868f
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h
@@ -0,0 +1,93 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_EMBED_IN_RD_H_
+#define FUNCTIONS_EMBED_IN_RD_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Embed_in_Rd
+ * \brief Embedding of an implicit manifold in a higher dimension.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ */
+template <class Function_>
+struct Embed_in_Rd {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd x = p;
+ Eigen::VectorXd x_k(fun_.amb_d()), x_rest(d_ - fun_.amb_d());
+ for (std::size_t i = 0; i < fun_.amb_d(); ++i) x_k(i) = x(i);
+ for (std::size_t i = fun_.amb_d(); i < d_; ++i) x_rest(i - fun_.amb_d()) = x(i);
+ Eigen::VectorXd result = fun_(x_k);
+ result.conservativeResize(this->cod_d());
+ for (std::size_t i = fun_.cod_d(); i < this->cod_d(); ++i) result(i) = x_rest(i - fun_.cod_d());
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return d_; }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return d_ - (fun_.amb_d() - fun_.cod_d()); }
+
+ /** \brief Returns a point on the zero-set of the embedded function. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = fun_.seed();
+ result.conservativeResize(d_);
+ for (std::size_t l = fun_.amb_d(); l < d_; ++l) result(l) = 0;
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the embedding function.
+ *
+ * @param[in] function The function to be embedded in higher dimension.
+ * @param[in] d Embedding dimension.
+ */
+ Embed_in_Rd(const Function_& function, std::size_t d) : fun_(function), d_(d) {}
+
+ private:
+ Function_ fun_;
+ std::size_t d_;
+};
+
+/**
+ * \brief Static constructor of an embedding function.
+ *
+ * @param[in] function The function to be embedded in higher dimension.
+ * @param[in] d Embedding dimension.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Embed_in_Rd<Function_> make_embedding(const Function_& function, std::size_t d) {
+ return Embed_in_Rd<Function_>(function, d);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h
new file mode 100644
index 00000000..8911f990
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h
@@ -0,0 +1,110 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_SM_IN_RD_H_
+#define FUNCTIONS_FUNCTION_SM_IN_RD_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_Sm_in_Rd
+ * \brief A class for the function that defines an m-dimensional implicit sphere embedded
+ * in the d-dimensional Euclidean space.
+ */
+struct Function_Sm_in_Rd {
+ /** \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd x = p;
+ for (std::size_t i = 0; i < d_; ++i) x(i) -= center_[i];
+ Eigen::VectorXd result = Eigen::VectorXd::Zero(k_);
+ for (std::size_t i = 0; i < m_ + 1; ++i) result(0) += x(i) * x(i);
+ result(0) -= r_ * r_;
+ for (std::size_t j = 1; j < k_; ++j) result(j) = x(m_ + j);
+ return result;
+ }
+
+ /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+ std::size_t amb_d() const { return d_; };
+
+ /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+ std::size_t cod_d() const { return k_; };
+
+ /** \brief Returns a point on the sphere. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = Eigen::VectorXd::Zero(d_);
+ result(0) += r_;
+ for (std::size_t i = 0; i < d_; ++i) result(i) += center_[i];
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+ * in the d-dimensional Euclidean space.
+ *
+ * @param[in] r The radius of the sphere.
+ * @param[in] m The dimension of the sphere.
+ * @param[in] d The ambient dimension of the sphere.
+ * @param[in] center The center of the sphere.
+ */
+ Function_Sm_in_Rd(double r, std::size_t m, std::size_t d, Eigen::VectorXd center)
+ : m_(m), k_(d - m), d_(d), r_(r), center_(center) {}
+
+ /**
+ * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+ * in the d-dimensional Euclidean space centered at the origin.
+ *
+ * @param[in] r The radius of the sphere.
+ * @param[in] m The dimension of the sphere.
+ * @param[in] d The ambient dimension of the sphere.
+ */
+ Function_Sm_in_Rd(double r, std::size_t m, std::size_t d)
+ : m_(m), k_(d - m), d_(d), r_(r), center_(Eigen::VectorXd::Zero(d_)) {}
+
+ /**
+ * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+ * in the (m+1)-dimensional Euclidean space.
+ *
+ * @param[in] r The radius of the sphere.
+ * @param[in] m The dimension of the sphere.
+ * @param[in] center The center of the sphere.
+ */
+ Function_Sm_in_Rd(double r, std::size_t m, Eigen::VectorXd center)
+ : m_(m), k_(1), d_(m_ + 1), r_(r), center_(center) {}
+
+ /**
+ * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+ * in the (m+1)-dimensional Euclidean space centered at the origin.
+ *
+ * @param[in] r The radius of the sphere.
+ * @param[in] m The dimension of the sphere.
+ */
+ Function_Sm_in_Rd(double r, std::size_t m) : m_(m), k_(1), d_(m_ + 1), r_(r), center_(Eigen::VectorXd::Zero(d_)) {}
+
+ Function_Sm_in_Rd(const Function_Sm_in_Rd& rhs) : Function_Sm_in_Rd(rhs.r_, rhs.m_, rhs.d_, rhs.center_) {}
+
+ private:
+ std::size_t m_, k_, d_;
+ double r_;
+ Eigen::VectorXd center_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h
new file mode 100644
index 00000000..b29f0906
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h
@@ -0,0 +1,91 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_
+#define FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_affine_plane_in_Rd
+ * \brief A class for the function that defines an m-dimensional implicit affine plane
+ * embedded in d-dimensional Euclidean space.
+ */
+struct Function_affine_plane_in_Rd {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result = normal_matrix_.transpose() * (p - off_);
+ return result;
+ }
+
+ /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+ std::size_t amb_d() const { return d_; };
+
+ /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+ std::size_t cod_d() const { return k_; };
+
+ /** \brief Returns a point on the affine plane. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = off_;
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines an m-dimensional implicit affine
+ * plane in the d-dimensional Euclidean space.
+ *
+ * @param[in] normal_matrix A normal matrix of the affine plane. The number of rows should
+ * correspond to the ambient dimension, the number of columns should 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();
+ }
+
+ private:
+ 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..620446da
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h
@@ -0,0 +1,80 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_
+#define FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_
+
+#include <cstdlib> // for std::size_t
+#include <cmath> // for std::pow
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_chair_in_R3
+ * \brief A class that encodes the function, the zero-set of which is a so-called
+ * "chair" surface embedded in R^3.
+ */
+struct Function_chair_in_R3 {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2];
+ Eigen::VectorXd result(cod_d());
+ result(0) = std::pow(x * x + y * y + z * z - a_ * k_ * k_, 2) -
+ b_ * ((z - k_) * (z - k_) - 2 * x * x) * ((z + k_) * (z + k_) - 2 * y * y);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return 3; }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return 1; }
+
+ /** \brief Returns a point on the surface. */
+ Eigen::VectorXd seed() const {
+ double t1 = a_ - b_;
+ double discr = t1 * t1 - (1.0 - b_) * (a_ * a_ - b_);
+ double z0 = k_ * std::sqrt((t1 + std::sqrt(discr)) / (1 - b_));
+ Eigen::Vector3d result(off_[0], off_[1], z0 + off_[2]);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines the 'chair' surface
+ * embedded in R^3.
+ *
+ * @param[in] a A numerical parameter.
+ * @param[in] b A numerical parameter.
+ * @param[in] k A numerical parameter.
+ * @param[in] off Offset vector.
+ */
+ Function_chair_in_R3(double a = 0.8, double b = 0.4, double k = 1.0, Eigen::Vector3d off = Eigen::Vector3d::Zero())
+ : a_(a), b_(b), k_(k), off_(off) {}
+
+ protected:
+ double a_, b_, k_;
+ Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
+
+// (x^2 + y^2 + z^2 - a*k^2)^2 - b*((z-k)^2 - 2*x^2)*((z+k)^2 - 2*y^2)
+// sqrt(k/(1-b))*sqrt(a-b + sqrt((a-b)^2 - (1-b)*(a^2 - b)*k^2))
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h
new file mode 100644
index 00000000..f73c4280
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h
@@ -0,0 +1,69 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_IRON_IN_R3_H_
+#define FUNCTIONS_FUNCTION_IRON_IN_R3_H_
+
+#include <cstdlib> // for std::size_t
+#include <cmath> // for std::pow
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_iron_in_R3
+ * \brief A class that encodes the function, the zero-set of which is a surface
+ * embedded in R^3 that ressembles an iron.
+ */
+struct Function_iron_in_R3 {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ double x = p(0), y = p(1), z = p(2);
+ Eigen::VectorXd result(cod_d());
+ result(0) = -std::pow(x, 6) / 300. - std::pow(y, 6) / 300. - std::pow(z, 6) / 300. + x * y * y * z / 2.1 + y * y +
+ std::pow(z - 2, 4) - 1;
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return 3; };
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return 1; };
+
+ /** \brief Returns a point on the surface. */
+ Eigen::VectorXd seed() const {
+ Eigen::Vector3d result(std::pow(4500, 1. / 6), 0, 0);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines a surface embedded in R^3
+ * that ressembles an iron.
+ *
+ * @param[in] off Offset vector.
+ */
+ Function_iron_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero()) : off_(off) {}
+
+ private:
+ Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h
new file mode 100644
index 00000000..beb41e00
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h
@@ -0,0 +1,85 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_
+#define FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_
+
+#include <cstdlib> // for std::size_t
+#include <cmath> // for std::sqrt
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_lemniscate_revolution_in_R3
+ * \brief A class that encodes the function, the zero-set of which is a surface of revolution
+ * around the x axis based on the lemniscate of Bernoulli embedded in R^3.
+ */
+struct Function_lemniscate_revolution_in_R3 {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2];
+ Eigen::VectorXd result(cod_d());
+ double x2 = x * x, y2 = y * y, z2 = z * z, a2 = a_ * a_;
+ double t1 = x2 + y2 + z2;
+ result(0) = t1 * t1 - 2 * a2 * (x2 - y2 - z2);
+ return result;
+ }
+
+ /** \brief Returns the (ambient) domain dimension.*/
+ std::size_t amb_d() const { return 3; };
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return 1; };
+
+ /** \brief Returns a point on the surface. This seed point is only one of
+ * two necessary seed points for the manifold tracing algorithm.
+ * See the method seed2() for the other point.
+ */
+ Eigen::VectorXd seed() const {
+ Eigen::Vector3d result(std::sqrt(2 * a_) + off_[0], off_[1], off_[2]);
+ return result;
+ }
+
+ /** \brief Returns a point on the surface. This seed point is only one of
+ * two necessary seed points for the manifold tracing algorithm.
+ * See the method seed() for the other point.
+ */
+ Eigen::VectorXd seed2() const {
+ Eigen::Vector3d result(-std::sqrt(2 * a_) + off_[0], off_[1], off_[2]);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines a surface of revolution
+ * around the x axis based on the lemniscate of Bernoulli embedded in R^3.
+ *
+ * @param[in] a A numerical parameter.
+ * @param[in] off Offset vector.
+ */
+ Function_lemniscate_revolution_in_R3(double a = 1, Eigen::Vector3d off = Eigen::Vector3d::Zero())
+ : a_(a), off_(off) {}
+
+ private:
+ double a_;
+ Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h
new file mode 100644
index 00000000..11b379f3
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h
@@ -0,0 +1,79 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_
+#define FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_moment_curve_in_Rd
+ * \brief A class for the function that defines an implicit moment curve
+ * in the d-dimensional Euclidean space.
+ */
+struct Function_moment_curve_in_Rd {
+ /** \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result(k_);
+ for (std::size_t i = 1; i < d_; ++i) result(i - 1) = p(i) - p(0) * p(i - 1);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension.. */
+ std::size_t amb_d() const { return d_; };
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return k_; };
+
+ /** \brief Returns a point on the moment curve. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = Eigen::VectorXd::Zero(d_);
+ return result;
+ }
+
+ /**
+ * \brief 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) {}
+
+ private:
+ 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..b54d3c74
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h
@@ -0,0 +1,71 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_TORUS_IN_R3_H_
+#define FUNCTIONS_FUNCTION_TORUS_IN_R3_H_
+
+#include <cstdlib> // for std::size_t
+#include <cmath> // for std::sqrt
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_torus_in_R3
+ * \brief A class that encodes the function, the zero-set of which is a torus
+ * surface embedded in R^3.
+ */
+struct Function_torus_in_R3 {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2];
+ Eigen::VectorXd result(cod_d());
+ result(0) = (z * z + (std::sqrt(x * x + y * y) - r_) * (std::sqrt(x * x + y * y) - r_) - R_ * R_);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return 3; };
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return 1; };
+
+ /** \brief Returns a point on the surface. */
+ Eigen::VectorXd seed() const {
+ Eigen::Vector3d result(R_ + r_ + off_[0], off_[1], off_[2]);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines a torus embedded in R^3.
+ *
+ * @param[in] R The outer radius of the torus.
+ * @param[in] r The inner radius of the torus.
+ * @param[in] off Offset vector.
+ */
+ Function_torus_in_R3(double R = 1, double r = 0.5, Eigen::Vector3d off = Eigen::Vector3d::Zero())
+ : R_(R), r_(r), off_(off) {}
+
+ private:
+ double R_, r_;
+ Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h
new file mode 100644
index 00000000..df1f1eec
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h
@@ -0,0 +1,78 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_
+#define FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Function_whitney_umbrella_in_R3
+ * \brief A class that encodes the function, the zero-set of which is the Whitney umbrella
+ * surface embedded in R^3.
+ */
+struct Function_whitney_umbrella_in_R3 {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ double x = p(0) - off_[0], y = p(1) - off_[1], z = p(2) - off_[2];
+ Eigen::VectorXd result(cod_d());
+ result(0) = x * x - y * y * z;
+ return result;
+ }
+
+ /** \brief Returns the (ambient) domain dimension.*/
+ std::size_t amb_d() const { return 3; };
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return 1; };
+
+ /** \brief Returns a point on the surface. This seed point is only one of
+ * two necessary seed points for the manifold tracing algorithm.
+ * See the method seed2() for the other point.
+ */
+ Eigen::VectorXd seed() const {
+ Eigen::Vector3d result(1 + off_[0], 1 + off_[1], 1 + off_[2]);
+ return result;
+ }
+
+ /** \brief Returns a point on the surface. This seed point is only one of
+ * two necessary seed points for the manifold tracing algorithm.
+ * See the method seed() for the other point.
+ */
+ Eigen::VectorXd seed2() const {
+ Eigen::Vector3d result(-1 + off_[0], -1 + off_[1], 1 + off_[2]);
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the function that defines the Whitney umbrella in R^3.
+ *
+ * @param[in] off Offset vector.
+ */
+ Function_whitney_umbrella_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero()) : off_(off) {}
+
+ private:
+ Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h
new file mode 100644
index 00000000..82e25bb9
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h
@@ -0,0 +1,88 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_LINEAR_TRANSFORMATION_H_
+#define FUNCTIONS_LINEAR_TRANSFORMATION_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Linear_transformation
+ * \brief Transforms the zero-set of the function by a given linear transformation.
+ * The underlying function corresponds to f(M*x), where M is the transformation matrix.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ */
+template <class Function_>
+struct Linear_transformation {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result = fun_(matrix_.householderQr().solve(p));
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return fun_.amb_d(); }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return fun_.cod_d(); }
+
+ /** \brief Returns a point on the zero-set. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = fun_.seed();
+ result = matrix_ * result;
+ return result;
+ }
+
+ /**
+ * \brief Constructor of a linearly transformed function.
+ *
+ * @param[in] function The function to be linearly transformed.
+ * @param[in] matrix The transformation matrix. Its dimension should be d*d,
+ * where d is the domain (ambient) dimension of 'function'.
+ */
+ Linear_transformation(const Function_& function, const Eigen::MatrixXd& matrix) : fun_(function), matrix_(matrix) {}
+
+ private:
+ Function_ fun_;
+ Eigen::MatrixXd matrix_;
+};
+
+/**
+ * \brief Static constructor of a linearly transformed function.
+ *
+ * @param[in] function The function to be linearly transformed.
+ * @param[in] matrix The transformation matrix. Its dimension should be d*d,
+ * where d is the domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Linear_transformation<Function_> make_linear_transformation(const Function_& function, const Eigen::MatrixXd& matrix) {
+ return Linear_transformation<Function_>(function, matrix);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h
new file mode 100644
index 00000000..fdf07f27
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h
@@ -0,0 +1,84 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_NEGATION_H_
+#define FUNCTIONS_NEGATION_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ *\class Negation
+ * \brief Constructs the "minus" function. The zero-set is the same, but
+ * the values at other points are the negative of their original value.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ */
+template <class Function_>
+struct Negation {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result = -fun_(p);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return fun_.amb_d(); }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return fun_.cod_d(); }
+
+ /** \brief Returns a point on the zero-set. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = fun_.seed();
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the negative function.
+ *
+ * @param[in] function The function to be negated.
+ */
+ Negation(const Function_& function) : fun_(function) {}
+
+ private:
+ Function_ fun_;
+};
+
+/**
+ * \brief Static constructor of the negative function.
+ *
+ * @param[in] function The function to be translated.
+ * domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Negation<Function_> negation(const Function_& function) {
+ return Negation<Function_>(function);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h
new file mode 100644
index 00000000..22071d6d
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h
@@ -0,0 +1,111 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_PL_APPROXIMATION_H_
+#define FUNCTIONS_PL_APPROXIMATION_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class PL_approximation
+ * \brief Constructs a piecewise-linear approximation of a function induced by
+ * an ambient triangulation.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ * \tparam Triangulation The triangulation template parameter. Should be a model of
+ * the concept TriangulationForManifoldTracing.
+ */
+template <class Function_, class Triangulation_>
+struct PL_approximation {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ std::size_t cod_d = this->cod_d();
+ std::size_t amb_d = this->amb_d();
+ auto s = tr_.locate_point(p);
+ Eigen::MatrixXd matrix(cod_d, s.dimension() + 1);
+ Eigen::MatrixXd vertex_matrix(amb_d + 1, s.dimension() + 1);
+ for (std::size_t i = 0; i < s.dimension() + 1; ++i) vertex_matrix(0, i) = 1;
+ std::size_t j = 0;
+ for (auto v : s.vertex_range()) {
+ Eigen::VectorXd pt_v = tr_.cartesian_coordinates(v);
+ Eigen::VectorXd fun_v = fun_(pt_v);
+ for (std::size_t i = 1; i < amb_d + 1; ++i) vertex_matrix(i, j) = pt_v(i - 1);
+ for (std::size_t i = 0; i < cod_d; ++i) matrix(i, j) = fun_v(i);
+ j++;
+ }
+ assert(j == s.dimension() + 1);
+ Eigen::VectorXd z(amb_d + 1);
+ z(0) = 1;
+ for (std::size_t i = 1; i < amb_d + 1; ++i) z(i) = p(i - 1);
+ Eigen::VectorXd lambda = vertex_matrix.colPivHouseholderQr().solve(z);
+ Eigen::VectorXd result = matrix * lambda;
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return fun_.amb_d(); }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return fun_.cod_d(); }
+
+ /** \brief Returns a point on the zero-set. */
+ Eigen::VectorXd seed() const {
+ // TODO: not finished. Should use an oracle.
+ return Eigen::VectorXd(amb_d());
+ }
+
+ /**
+ * \brief Constructor of the piecewise-linear approximation of a function
+ * induced by an ambient triangulation.
+ *
+ * @param[in] function The function.
+ * @param[in] triangulation The ambient triangulation.
+ */
+ PL_approximation(const Function_& function, const Triangulation_& triangulation)
+ : fun_(function), tr_(triangulation) {}
+
+ private:
+ Function_ fun_;
+ Triangulation_ tr_;
+};
+
+/**
+ * \brief Static constructor of the piecewise-linear approximation of a function
+ * induced by an ambient triangulation.
+ *
+ * @param[in] function The function.
+ * @param[in] triangulation The ambient triangulation.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_, class Triangulation_>
+PL_approximation<Function_, Triangulation_> make_pl_approximation(const Function_& function,
+ const Triangulation_& triangulation) {
+ return PL_approximation<Function_, Triangulation_>(function, triangulation);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h
new file mode 100644
index 00000000..cbe65abe
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h
@@ -0,0 +1,89 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_TRANSLATE_H_
+#define FUNCTIONS_TRANSLATE_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Translate
+ * \brief Translates the zero-set of the function by a vector.
+ * The underlying function corresponds to f(x-off), where off is the offset vector.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ */
+template <class Function_>
+struct Translate {
+ /**
+ * \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+ Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+ Eigen::VectorXd result = fun_(p - off_);
+ return result;
+ }
+
+ /** \brief Returns the domain (ambient) dimension. */
+ std::size_t amb_d() const { return fun_.amb_d(); }
+
+ /** \brief Returns the codomain dimension. */
+ std::size_t cod_d() const { return fun_.cod_d(); }
+
+ /** \brief Returns a point on the zero-set. */
+ Eigen::VectorXd seed() const {
+ Eigen::VectorXd result = fun_.seed();
+ result += off_;
+ return result;
+ }
+
+ /**
+ * \brief Constructor of the translated function.
+ *
+ * @param[in] function The function to be translated.
+ * @param[in] off The offset vector. The dimension should correspond to the
+ * domain (ambient) dimension of 'function'.
+ */
+ Translate(const Function_& function, const Eigen::VectorXd& off) : fun_(function), off_(off) {}
+
+ private:
+ Function_ fun_;
+ Eigen::VectorXd off_;
+};
+
+/**
+ * \brief Static constructor of a translated function.
+ *
+ * @param[in] function The function to be translated.
+ * @param[in] off The offset vector. The dimension should correspond to the
+ * domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Translate<Function_> translate(const Function_& function, Eigen::VectorXd off) {
+ return Translate<Function_>(function, off);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h
new file mode 100644
index 00000000..6a896e94
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h
@@ -0,0 +1,72 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
+#define FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
+
+#include <cstdlib> // for std::size_t
+#include <cmath> // for std::cos, std::sin
+#include <random> // for std::uniform_real_distribution, std::random_device
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+
+#include <boost/math/constants/constants.hpp> // for PI value
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief Generates a uniform random orthogonal matrix using the "subgroup algorithm" by
+ * Diaconis & Shashahani.
+ * \details Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Uniform_random_rotation_matrices.
+ * The idea: take a random rotation matrix of dimension d-1, embed it
+ * as a d*d matrix M with the last column (0,...,0,1).
+ * Pick a random vector v on a sphere S^d. rotate the matrix M so that its last column is v.
+ * The determinant of the matrix can be either 1 or -1
+ */
+// Note: the householderQR operation at the end seems to take a lot of time at compilation.
+// The CGAL headers are another source of long compilation time.
+Eigen::MatrixXd random_orthogonal_matrix(std::size_t d) {
+ typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
+ typedef typename Kernel::Point_d Point_d;
+ if (d == 1) return Eigen::VectorXd::Constant(1, 1.0);
+ if (d == 2) {
+ // 0. < alpha < 2 Pi
+ std::uniform_real_distribution<double> unif(0., 2 * boost::math::constants::pi<double>());
+ std::random_device rand_dev;
+ std::mt19937 rand_engine(rand_dev());
+ double alpha = unif(rand_engine);
+
+ Eigen::Matrix2d rot;
+ rot << std::cos(alpha), -std::sin(alpha), std::sin(alpha), cos(alpha);
+ return rot;
+ }
+ Eigen::MatrixXd low_dim_rot = random_orthogonal_matrix(d - 1);
+ Eigen::MatrixXd rot(d, d);
+ Point_d v = *CGAL::Random_points_on_sphere_d<Point_d>(d, 1);
+ for (std::size_t i = 0; i < d; ++i) rot(i, 0) = v[i];
+ for (std::size_t i = 0; i < d - 1; ++i)
+ for (std::size_t j = 1; j < d - 1; ++j) rot(i, j) = low_dim_rot(i, j - 1);
+ for (std::size_t j = 1; j < d; ++j) rot(d - 1, j) = 0;
+ rot = rot.householderQr()
+ .householderQ(); // a way to do Gram-Schmidt, see https://forum.kde.org/viewtopic.php?f=74&t=118568#p297246
+ return rot;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
new file mode 100644
index 00000000..ca08f629
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
@@ -0,0 +1,60 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_MESH_MEDIT_H_
+#define IO_MESH_MEDIT_H_
+
+#include <Eigen/Dense>
+
+#include <vector>
+#include <utility> // for std::pair
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Mesh_medit
+ * \brief Structure to store a mesh that can be output in Medit .mesh file format
+ * using the output_meshes_to_medit method.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Mesh_medit {
+ /** \brief Type of a range of vertices. */
+ typedef std::vector<Eigen::VectorXd> Vertex_points;
+ /** \brief Type of a mesh element.
+ * A pair consisting of a vector of vertex indices of type std::size_t
+ * and of an integer that represents the common reference number for
+ * the mesh elements of this type. */
+ typedef std::pair<std::vector<std::size_t>, std::size_t> Mesh_element;
+ /** \brief Type of a range of mesh elements. */
+ typedef std::vector<Mesh_element> Mesh_elements;
+ /** \brief Type of a range of scalar field . */
+ typedef std::vector<double> Scalar_field_range;
+
+ /** \brief Range of vertices of type Eigen::VectorXd to output. */
+ Vertex_points vertex_points;
+ /** \brief Range of edges. */
+ Mesh_elements edges;
+ /** \brief Range of triangles. */
+ Mesh_elements triangles;
+ /** \brief Range of tetrahedra. */
+ Mesh_elements tetrahedra;
+ /** \brief Range of scalar values over triangles. */
+ Scalar_field_range triangles_scalar_range;
+ /** \brief Range of scalar values over tetrahedra. */
+ Scalar_field_range tetrahedra_scalar_range;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
new file mode 100644
index 00000000..9750f366
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
@@ -0,0 +1,171 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+#define IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+
+#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES
+#include <gudhi/IO/Mesh_medit.h>
+
+#include <Eigen/Dense>
+
+#include <cstdlib> // for std::size_t
+#include <map>
+#include <set>
+#include <string>
+#include <utility> // for std::make_pair
+#include <algorithm> // for std::min
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+struct Configuration {
+ Configuration(bool t_edges, bool t_triangles, bool t_tetrahedra, std::size_t r_edges, std::size_t r_triangles,
+ std::size_t r_tetrahedra)
+ : toggle_edges(t_edges),
+ toggle_triangles(t_triangles),
+ toggle_tetrahedra(t_tetrahedra),
+ ref_edges(r_edges),
+ ref_triangles(r_triangles),
+ ref_tetrahedra(r_tetrahedra) {}
+
+ Configuration() {}
+
+ bool toggle_edges = true, toggle_triangles = true, toggle_tetrahedra = true;
+ std::size_t ref_edges = 1, ref_triangles = 1, ref_tetrahedra = 1;
+};
+
+template <class Hasse_cell, class Simplex_cell_map>
+void populate_mesh(Mesh_medit& output, Simplex_cell_map& sc_map, Configuration configuration, std::size_t amb_d,
+ std::map<Hasse_cell*, std::size_t> vi_map) {
+ using Mesh_element_vertices = Mesh_medit::Mesh_elements::value_type::first_type;
+ std::map<Hasse_cell*, std::size_t> ci_map;
+ std::size_t index = vi_map.size() + 1; // current size of output.vertex_points
+ if (sc_map.size() >= 3)
+ for (const auto& sc_pair : sc_map[2]) {
+ Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+ std::set<std::size_t> vertex_indices;
+ Hasse_cell* cell = sc_pair.second;
+ for (const auto& ei_pair : cell->get_boundary())
+ for (const auto& vi_pair : ei_pair.first->get_boundary()) vertex_indices.emplace(vi_map[vi_pair.first]);
+ for (const std::size_t& v : vertex_indices) barycenter += output.vertex_points[v - 1];
+ ci_map.emplace(cell, index++);
+ output.vertex_points.emplace_back((1. / vertex_indices.size()) * barycenter);
+#ifdef DEBUG_TRACES
+ std::string vlist = " (" + std::to_string(index - 1) + ")";
+ for (const std::size_t& v : vertex_indices) vlist += " " + std::to_string(v);
+ cell_vlist_map.emplace(to_string(cell), vlist);
+#endif
+ }
+
+ if (configuration.toggle_edges && sc_map.size() >= 2)
+ for (const auto& sc_pair : sc_map[1]) {
+ Hasse_cell* edge_cell = sc_pair.second;
+ Mesh_element_vertices edge;
+ for (const auto& vi_pair : edge_cell->get_boundary()) edge.push_back(vi_map[vi_pair.first]);
+ output.edges.emplace_back(edge, configuration.ref_edges);
+#ifdef DEBUG_TRACES
+ std::string vlist;
+ for (const std::size_t& v : edge) vlist += " " + std::to_string(v);
+ cell_vlist_map.emplace(to_string(edge_cell), vlist);
+#endif
+ }
+
+ if (configuration.toggle_triangles && sc_map.size() >= 3)
+ for (const auto& sc_pair : sc_map[2]) {
+ for (const auto& ei_pair : sc_pair.second->get_boundary()) {
+ Mesh_element_vertices triangle(1, ci_map[sc_pair.second]);
+ for (const auto& vi_pair : ei_pair.first->get_boundary()) triangle.push_back(vi_map[vi_pair.first]);
+ output.triangles.emplace_back(triangle, configuration.ref_triangles);
+ }
+ }
+
+ if (configuration.toggle_tetrahedra && sc_map.size() >= 4)
+ for (const auto& sc_pair : sc_map[3]) {
+ Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+ std::set<std::size_t> vertex_indices;
+ Hasse_cell* cell = sc_pair.second;
+ for (const auto& ci_pair : cell->get_boundary())
+ for (const auto& ei_pair : ci_pair.first->get_boundary())
+ for (const auto& vi_pair : ei_pair.first->get_boundary()) vertex_indices.emplace(vi_map[vi_pair.first]);
+ for (const std::size_t& v : vertex_indices) barycenter += output.vertex_points[v - 1];
+ output.vertex_points.emplace_back((1. / vertex_indices.size()) * barycenter);
+#ifdef DEBUG_TRACES
+ std::string vlist = " (" + std::to_string(index) + ")";
+ for (const std::size_t& v : vertex_indices) vlist += " " + std::to_string(v);
+ cell_vlist_map.emplace(to_string(cell), vlist);
+#endif
+
+ for (const auto& ci_pair : cell->get_boundary())
+ for (const auto& ei_pair : ci_pair.first->get_boundary()) {
+ Mesh_element_vertices tetrahedron = {index, ci_map[sc_pair.second]};
+ for (const auto& vi_pair : ei_pair.first->get_boundary()) tetrahedron.push_back(vi_map[vi_pair.first]);
+ output.tetrahedra.emplace_back(tetrahedron, configuration.ref_tetrahedra);
+ }
+ index++;
+ }
+}
+
+/** @brief Builds a Gudhi::coxeter_triangulation::Mesh_medit from a Gudhi::coxeter_triangulation::Cell_complex
+ *
+ * @ingroup coxeter_triangulation
+ */
+template <class Cell_complex>
+Mesh_medit build_mesh_from_cell_complex(const Cell_complex& cell_complex,
+ Configuration i_configuration = Configuration(),
+ Configuration b_configuration = Configuration()) {
+ using Hasse_cell = typename Cell_complex::Hasse_cell;
+ Mesh_medit output;
+ std::map<Hasse_cell*, std::size_t> vi_map; // one for vertices, other for 2d-cells
+ std::size_t index = 1; // current size of output.vertex_points
+
+ if (cell_complex.cell_point_map().empty()) return output;
+ std::size_t amb_d = std::min((int)cell_complex.cell_point_map().begin()->second.size(), 3);
+
+ for (const auto& cp_pair : cell_complex.cell_point_map()) {
+#ifdef DEBUG_TRACES
+ std::string vlist;
+ vlist += " " + std::to_string(index);
+ cell_vlist_map.emplace(to_string(cp_pair.first), vlist);
+#endif
+ vi_map.emplace(cp_pair.first, index++);
+ output.vertex_points.push_back(cp_pair.second);
+ output.vertex_points.back().conservativeResize(amb_d);
+ }
+
+ populate_mesh(output, cell_complex.interior_simplex_cell_maps(), i_configuration, amb_d, vi_map);
+#ifdef DEBUG_TRACES
+ for (const auto& sc_map : cell_complex.interior_simplex_cell_maps())
+ for (const auto& sc_pair : sc_map) {
+ std::string simplex = "I" + to_string(sc_pair.first);
+ std::string cell = to_string(sc_pair.second);
+ std::string vlist = cell_vlist_map.at(cell).substr(1);
+ simplex_vlist_map.emplace(simplex, vlist);
+ }
+#endif
+ populate_mesh(output, cell_complex.boundary_simplex_cell_maps(), b_configuration, amb_d, vi_map);
+#ifdef DEBUG_TRACES
+ for (const auto& sc_map : cell_complex.boundary_simplex_cell_maps())
+ for (const auto& sc_pair : sc_map) {
+ std::string simplex = "B" + to_string(sc_pair.first);
+ std::string cell = to_string(sc_pair.second);
+ std::string vlist = cell_vlist_map.at(cell).substr(1);
+ simplex_vlist_map.emplace(simplex, vlist);
+ }
+#endif
+ return output;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
new file mode 100644
index 00000000..a2995738
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
@@ -0,0 +1,550 @@
+#ifndef IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_
+#define IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_
+
+#ifdef DEBUG_TRACES // All this part of code can be skipped if DEBUG_TRACES are not ON - cmake -DDEBUG_TRACES=ON .
+
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <list>
+#include <string>
+#include <regex>
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+template <class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector) {
+ os << "(";
+ if (vector.empty()) {
+ os << ")";
+ return os;
+ }
+ auto v_it = vector.begin();
+ os << *v_it++;
+ for (; v_it != vector.end(); ++v_it) os << ", " << *v_it;
+ os << ")";
+ return os;
+}
+
+/* A class to make the vector horizontal instead of vertical */
+struct Straighten {
+ Straighten(const Eigen::VectorXd& vector) : vector_(vector) {}
+ const Eigen::VectorXd& vector_;
+};
+
+std::ostream& operator<<(std::ostream& os, const Straighten& str) {
+ std::size_t size = str.vector_.size();
+ os << "(" << str.vector_(0);
+ if (size == 0) {
+ os << ")";
+ return os;
+ }
+ for (std::size_t i = 1; i < size; ++i) os << ", " << str.vector_(i);
+ os << ")";
+ return os;
+}
+
+std::string id_from_simplex(const std::string& simplex) {
+ std::regex r("\\s+"), r2("\\(|\\)|\\{|\\}"), r3(","), r4("\\["), r5("\\]");
+ std::string output = std::regex_replace(simplex, r, "");
+ output = std::regex_replace(output, r2, ":");
+ output = std::regex_replace(output, r3, ".");
+ output = std::regex_replace(output, r4, "_");
+ output = std::regex_replace(output, r5, "");
+ return output;
+}
+
+template <typename T>
+std::string to_string(const T& t) {
+ std::ostringstream oss;
+ oss << t;
+ return oss.str();
+}
+
+struct MT_inserted_info {
+ std::string qr_face_, init_face_, qr_intersection_;
+ bool qr_success_, is_boundary_;
+ template <class Query_result, class Simplex_handle>
+ MT_inserted_info(const Query_result& qr, const Simplex_handle& face, bool is_boundary)
+ : qr_face_(to_string(face)),
+ init_face_(to_string(face)),
+ qr_intersection_(to_string(qr.intersection)),
+ qr_success_(qr.success),
+ is_boundary_(is_boundary) {}
+};
+std::list<MT_inserted_info> mt_seed_inserted_list, mt_inserted_list;
+
+struct CC_summary_info {
+ std::string face_, cell_;
+ template <class SC_pair>
+ CC_summary_info(const SC_pair& sc_pair) : face_(to_string(sc_pair.first)), cell_(to_string(sc_pair.second)) {}
+};
+using CC_summary_list = std::list<CC_summary_info>;
+std::vector<CC_summary_list> cc_interior_summary_lists, cc_boundary_summary_lists;
+
+struct CC_detail_info {
+ enum class Result_type { self, face, coface, inserted, join_single, join_is_face };
+ std::string simplex_, trigger_, init_simplex_;
+ Result_type status_;
+ bool join_trigger_ = false;
+ std::list<std::string> faces_, post_faces_, cofaces_;
+ template <class Simplex_handle>
+ CC_detail_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)) {}
+};
+using CC_detail_list = std::list<CC_detail_info>;
+std::vector<CC_detail_list> cc_interior_detail_lists, cc_boundary_detail_lists;
+std::vector<CC_detail_list> cc_interior_insert_detail_lists, cc_boundary_insert_detail_lists;
+
+struct CC_prejoin_info {
+ enum class Result_type { join_single, join_is_face, join_different, join_same };
+ std::string simplex_, join_;
+ std::vector<std::string> faces_;
+ std::size_t dimension_;
+ Result_type status_;
+ template <class Simplex_handle>
+ CC_prejoin_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)), dimension_(simplex.dimension()) {}
+};
+using CC_prejoin_list = std::list<CC_prejoin_info>;
+std::vector<CC_prejoin_list> cc_interior_prejoin_lists, cc_boundary_prejoin_lists;
+
+struct CC_join_info {
+ enum class Result_type { self, face, coface, inserted, join_single, join_is_face };
+ std::string simplex_, join_, trigger_;
+ Result_type status_;
+ std::list<std::string> boundary_faces_;
+ std::list<std::string> faces_, post_faces_, cofaces_;
+ template <class Simplex_handle>
+ CC_join_info(const Simplex_handle& simplex) : simplex_(to_string(simplex)) {}
+};
+bool join_switch = false;
+std::vector<CC_detail_list> cc_interior_join_detail_lists, cc_boundary_join_detail_lists;
+
+std::map<std::string, std::string> cell_vlist_map;
+std::map<std::string, std::string> simplex_vlist_map;
+
+std::ostringstream mt_ostream, vis_ostream;
+std::vector<std::ostringstream> cc_summary_ostream, cc_traces_ostream;
+
+std::string simplex_format(const std::string& simplex, bool is_boundary) {
+ std::string b_simplex = (is_boundary ? "B" : "I") + simplex;
+ std::string tooltiptext;
+ auto it = simplex_vlist_map.find(b_simplex);
+ if (it == simplex_vlist_map.end())
+ tooltiptext = "deleted";
+ else
+ tooltiptext = simplex_vlist_map.at(b_simplex);
+ return (std::string) "<a class=\"" + (is_boundary ? "boundary" : "interior") + "\" href=\"#" +
+ id_from_simplex(b_simplex) + "\">" + b_simplex + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+std::string simplex_format(const std::string& b_simplex) {
+ bool is_boundary = b_simplex[0] == 'B';
+ std::string tooltiptext;
+ auto it = simplex_vlist_map.find(b_simplex);
+ if (it == simplex_vlist_map.end())
+ tooltiptext = "deleted";
+ else
+ tooltiptext = simplex_vlist_map.at(b_simplex);
+ return (std::string) "<a class=\"" + (is_boundary ? "boundary" : "interior") + "\" href=\"#" +
+ id_from_simplex(b_simplex) + "\">" + b_simplex + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+void write_head(std::ofstream& ofs) {
+ ofs << " <head>\n"
+ << " <title>Cell complex debug trace</title>\n"
+ << " <style>\n"
+ << " a.boundary {\n"
+ << " position: relative;\n"
+ << " display: inline-block;\n"
+ << " color: darkred;\n"
+ << " background-color: lightgreen\n"
+ << " }\n"
+ << " a.interior {\n"
+ << " position: relative;\n"
+ << " display: inline-block;\n"
+ << " color: navy;\n"
+ << " background-color: yellow\n"
+ << " }\n"
+ << " .tooltiptext {\n"
+ << " visibility: hidden;\n"
+ << " width: 120px;\n"
+ << " background-color: #555;\n"
+ << " color: #fff;\n"
+ << " text-align: center;\n"
+ << " padding: 5px 0;\n"
+ << " border-radius: 6px;\n"
+ << " position: absolute;\n"
+ << " z-index: 1;\n"
+ << " bottom: 125%;\n"
+ << " left: 50%;\n"
+ << " margin-left: -60px;\n"
+ << " opacity: 0;\n"
+ << " transition: opacity 0.3s;\n"
+ << " }\n"
+ << " .boundary .tooltiptext::after {\n"
+ << " content: \"\";\n"
+ << " position: absolute;\n"
+ << " top: 100%;\n"
+ << " left: 50%;\n"
+ << " margin-left: -5px;\n"
+ << " border-width: 5px;\n"
+ << " border-style: solid;\n"
+ << " border-color: #555 transparent transparent transparent;\n"
+ << " }\n"
+ << " .interior .tooltiptext::after {\n"
+ << " content: \"\";\n"
+ << " position: absolute;\n"
+ << " top: 100%;\n"
+ << " left: 50%;\n"
+ << " margin-left: -5px;\n"
+ << " border-width: 5px;\n"
+ << " border-style: solid;\n"
+ << " border-color: #555 transparent transparent transparent;\n"
+ << " }\n"
+ << " .boundary:hover .tooltiptext {\n"
+ << " visibility: visible;\n"
+ << " opacity: 1;\n"
+ << " }\n"
+ << " .interior:hover .tooltiptext {\n"
+ << " visibility: visible;\n"
+ << " opacity: 1;\n"
+ << " }\n"
+ << " ul.nav {\n"
+ << " list-style-type: none;\n"
+ << " margin: 0;\n"
+ << " padding: 0;\n"
+ << " overflow: auto;\n"
+ << " background-color: #333;\n"
+ << " position: fixed;\n"
+ << " height: 100%;\n"
+ << " width: 15%;\n"
+ << " }\n"
+ << " ul.nav li a {\n"
+ << " display: block;\n"
+ << " color: white;\n"
+ << " text-align: left;\n"
+ << " padding: 14px 16px;\n"
+ << " text-decoration: none;\n"
+ << " }\n"
+ << " .active {\n"
+ << " background-color: #4CAF50;\n"
+ << " }\n"
+ << " div {\n"
+ << " margin-left: 15%;\n"
+ << " padding: 1px 16px\n"
+ << " }\n"
+ << " div.navi {\n"
+ << " margin-left: 0%;\n"
+ << " padding: 0px 0px\n"
+ << " }\n"
+ << " h1 {\n"
+ << " margin-left: 15%;\n"
+ << " padding: 1px 16px\n"
+ << " }\n"
+ << " </style>\n"
+ << " </head>\n";
+}
+
+void write_nav(std::ofstream& ofs) {
+ ofs << " <div class=\"navi\" style=\"margin-top:30px;background-color:#1abc9c;\">\n"
+ << " <ul class=\"nav\">\n"
+ << " <li><a href=\"#mant\">Manifold tracing</a></li>\n"
+ << " <li><a href=\"#cell\">Cell complex</a>\n"
+ << " <ul>\n";
+ for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+ ofs << " <li><a href=\"#dim" << i << "\">Dimension " << i << "</a>\n";
+ ofs << " <ul>\n";
+ ofs << " <li><a href=\"#dim" << i << "i\">Interior</a></li>\n";
+ if (i < cc_boundary_summary_lists.size()) {
+ ofs << " <li><a href=\"#dim" << i << "b\">Boundary</a></li>\n";
+ }
+ ofs << " </ul>\n";
+ ofs << " </li>\n";
+ }
+ ofs << " </ul>\n"
+ << " </li>\n"
+ << " <li><a href=\"#visu\">Visualization details</a></li>\n"
+ << " </ul>\n"
+ << " </div>\n";
+}
+
+void write_mt(std::ofstream& ofs) {
+ ofs << " <div id=\"mant\">\n";
+ ofs << " <h2> Manifold debug trace </h2>\n";
+ ofs << " <h3> Simplices inserted during the seed phase </h3>\n";
+ ofs << " <ul>\n";
+ for (const MT_inserted_info& mt_info : mt_seed_inserted_list) {
+ if (mt_info.qr_success_) {
+ ofs << " <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+ if (mt_info.qr_face_ != mt_info.init_face_)
+ ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+ ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+ } else
+ ofs << " <li>Failed to insert " << mt_info.init_face_ << "</li>\n";
+ }
+ ofs << " </ul>\n";
+ ofs << " <h3> Simplices inserted during the while loop phase </h3>\n";
+ ofs << " <ul>\n";
+ for (const MT_inserted_info& mt_info : mt_inserted_list) {
+ if (mt_info.qr_success_) {
+ ofs << " <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+ if (mt_info.qr_face_ != mt_info.init_face_)
+ ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+ ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+ } else
+ ofs << " <li>Failed to insert " << mt_info.init_face_ << ")</li>\n";
+ }
+ ofs << " </ul>\n";
+ ofs << " </div>\n";
+}
+
+void write_cc(std::ofstream& ofs) {
+ ofs << " <div id=\"cell\">\n"
+ << " <h2> Cell complex debug trace </h2>\n"
+ << " <p>Go to:</p>\n"
+ << " <ul>\n";
+ for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+ ofs << " <li><a href=\"#dim" << i << "\">Dimension " << i << "</a></li>\n";
+ }
+ ofs << " </ul>\n";
+ for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+ ofs << " <h3 id=\"dim" << i << "\"> Dimension " << i << "</h3>\n";
+ ofs << " <h4 id=\"dim" << i << "i\"> Summary for interior simplices</h4>\n";
+ if (i < cc_boundary_summary_lists.size()) ofs << " <p><a href=\"#dim" << i << "b\">Go to boundary</a></p>\n";
+ ofs << " <ul>\n";
+ for (const CC_summary_info& cc_info : cc_interior_summary_lists[i])
+ ofs << " <li id = \"" << id_from_simplex("I" + cc_info.face_) << "\">"
+ << simplex_format(cc_info.face_, false) << " cell =" << cc_info.cell_ << "</li>\n";
+ ofs << " </ul>\n";
+ ofs << " <h4> Prejoin state of the interior cells of dimension " << i << "</h4>\n";
+ auto prejoin_it = cc_interior_prejoin_lists[i].begin();
+ while (prejoin_it != cc_interior_prejoin_lists[i].end()) {
+ std::size_t j = prejoin_it->dimension_;
+ ofs << " <h5>" << j << "-dimensional ambient simplices</h5>\n";
+ ofs << " <ul>\n";
+ for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+ ofs << " <li>" << simplex_format(prejoin_it->simplex_, false)
+ << " join = " << simplex_format(prejoin_it->join_, false) << " boundary:\n"
+ << " <ul>\n";
+ for (const auto& face : prejoin_it->faces_) ofs << " <li>" << simplex_format(face) << "</li>";
+ ofs << " </ul>\n";
+ switch (prejoin_it->status_) {
+ case (CC_prejoin_info::Result_type::join_single):
+ ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, false)
+ << " as it has a single face.</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_is_face):
+ ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, false)
+ << " as its join " << simplex_format(prejoin_it->join_, false) << " is one of the faces.</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_different):
+ ofs << " <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, false)
+ << " and replaced by its join " << simplex_format(prejoin_it->join_, false) << ".</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_same):
+ ofs << " <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, false)
+ << ".</p>";
+ }
+ ofs << " </li>";
+ }
+ ofs << " </ul>\n";
+ }
+ ofs << " <h4> Details for interior simplices</h4>\n";
+ ofs << " <ul>\n";
+ for (const CC_detail_info& cc_info : cc_interior_detail_lists[i]) {
+ if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+ ofs << " <li style=\"color:magenta\" id = \"" << id_from_simplex("I" + cc_info.simplex_)
+ << "\"> Simplex " << simplex_format(cc_info.simplex_, false) << " has only one face ("
+ << simplex_format(cc_info.trigger_, false) << ") and is deleted.";
+ continue;
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+ ofs << " <li style=\"color:darkmagenta\" id = \"" << id_from_simplex("I" + cc_info.simplex_)
+ << "\"> The join of the simplex " << simplex_format(cc_info.simplex_, false) << " is one of its faces ("
+ << simplex_format(cc_info.trigger_, false) << "), hence it is is deleted.";
+ continue;
+ }
+ ofs << " <li> Insert_cell called for " << simplex_format(cc_info.simplex_, false) << "\n";
+ ofs << " <ul>\n";
+ // for (const std::string& cof: cc_info.faces_)
+ // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, false)
+ // << " is a face of " << simplex_format(cof, false) << "\n";
+ ofs << " </ul>\n";
+ ofs << " <ul>\n";
+ if (cc_info.status_ == CC_detail_info::Result_type::self) {
+ ofs << " <p><span style=\"color:blue\">The simplex " << simplex_format(cc_info.simplex_, false)
+ << " already exists in the cell complex!</span></p>\n";
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::face) {
+ ofs << " <p><span style=\"color:red\">The simplex " << simplex_format(cc_info.simplex_, false)
+ << " is a face of the simplex " << simplex_format(cc_info.trigger_, false) << "!</span><br>\n";
+ ofs << " <ul>\n";
+ for (const std::string post_face : cc_info.post_faces_)
+ ofs << " <li id = \"" << id_from_simplex("I" + post_face) << "\">"
+ << "Post deleting " << simplex_format(post_face, false) << "</li>\n";
+ ofs << " </ul>\n";
+ ofs << " </p>\n";
+ ofs << " <p id = \"" << id_from_simplex("I" + cc_info.trigger_) << "\">"
+ << "Deleting " << simplex_format(cc_info.trigger_, false) << "</p>\n";
+ }
+ // for (const std::string& fac: cc_info.cofaces_)
+ // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, false)
+ // << " is a coface of " << simplex_format(fac, false) << "\n";
+ if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+ ofs << " <p><span style=\"color:darkorange\">The simplex " << simplex_format(cc_info.simplex_, false)
+ << " is a coface of the simplex " << simplex_format(cc_info.trigger_, false) << "!</span><p>\n";
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+ ofs << " <p><span style=\"color:green\">Successfully inserted "
+ << simplex_format(cc_info.simplex_, false) << "!</span><p>\n";
+ }
+ ofs << " </ul>\n";
+ ofs << " </li>\n";
+ }
+ ofs << " </ul>\n";
+
+ if (i < cc_boundary_summary_lists.size()) {
+ ofs << " <h4 id=\"dim" << i << "b\"> Summary for boundary simplices</h4>\n";
+ ofs << " <p><a href=\"#dim" << i << "i\">Go to interior</a></p>\n";
+ ofs << " <ul>\n";
+ for (const CC_summary_info& cc_info : cc_boundary_summary_lists[i])
+ ofs << " <li id = \"" << id_from_simplex("B" + cc_info.face_) << "\">"
+ << simplex_format(cc_info.face_, true) << " cell =" << cc_info.cell_ << "</li>\n";
+ ofs << " </ul>\n";
+ ofs << " <h4> Prejoin state of the boundary cells of dimension " << i << "</h4>\n";
+ auto prejoin_it = cc_boundary_prejoin_lists[i].begin();
+ while (prejoin_it != cc_boundary_prejoin_lists[i].end()) {
+ std::size_t j = prejoin_it->dimension_;
+ ofs << " <h5>" << j << "-dimensional ambient simplices</h5>\n";
+ ofs << " <ul>\n";
+ for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+ ofs << " <li>" << simplex_format(prejoin_it->simplex_, true)
+ << " join = " << simplex_format(prejoin_it->join_, true) << " boundary:\n"
+ << " <ul>\n";
+ for (const auto& face : prejoin_it->faces_) ofs << " <li>" << simplex_format(face) << "</li>";
+ ofs << " </ul>\n";
+ switch (prejoin_it->status_) {
+ case (CC_prejoin_info::Result_type::join_single):
+ ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, true)
+ << " as it has a single face.</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_is_face):
+ ofs << " <p style=\"color: red\">Deleted " << simplex_format(prejoin_it->simplex_, true)
+ << " as its join " << simplex_format(prejoin_it->join_, true) << " is one of the faces.</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_different):
+ ofs << " <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, true)
+ << " and replaced by its join " << simplex_format(prejoin_it->join_, true) << ".</p>";
+ break;
+ case (CC_prejoin_info::Result_type::join_same):
+ ofs << " <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, true)
+ << ".</p>";
+ }
+ ofs << " </li>";
+ }
+ ofs << " </ul>\n";
+ }
+ }
+ if (i < cc_boundary_detail_lists.size()) {
+ ofs << " <h4> Details for boundary simplices</h4>\n"
+ << " <ul>\n";
+ for (const CC_detail_info& cc_info : cc_boundary_detail_lists[i]) {
+ if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+ ofs << " <li style=\"color:magenta\" id = \"" << id_from_simplex("B" + cc_info.simplex_)
+ << "\"> Simplex " << simplex_format(cc_info.simplex_, true) << " has only one face ("
+ << simplex_format(cc_info.trigger_, true) << ") and is deleted.";
+ continue;
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+ ofs << " <li style=\"color:darkmagenta\" id = \"" << id_from_simplex("B" + cc_info.simplex_)
+ << "\"> The join of the simplex " << simplex_format(cc_info.simplex_, true) << " is one of its faces ("
+ << simplex_format(cc_info.trigger_, true) << "), hence it is is deleted.";
+ continue;
+ }
+ ofs << " <li> Insert_simplex called on " << simplex_format(cc_info.simplex_, true);
+ ofs << " <ul>\n";
+ // for (const std::string& cof: cc_info.faces_)
+ // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, true)
+ // << " is a face of " << simplex_format(cof, true) << "\n";
+ ofs << " </ul>\n";
+ ofs << " <ul>\n";
+ if (cc_info.status_ == CC_detail_info::Result_type::self) {
+ ofs << " <p><span style=\"color:blue\">The simplex " << simplex_format(cc_info.simplex_, true)
+ << " already exists in the cell complex!</span></p>\n";
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::face) {
+ ofs << " <p><span style=\"color:red\">The simplex " << simplex_format(cc_info.simplex_, true)
+ << " is a face of the simplex " << simplex_format(cc_info.trigger_, true) << "!</span><br>\n";
+ ofs << " <ul>\n";
+ for (const std::string post_face : cc_info.post_faces_)
+ ofs << " <li id=\"" << id_from_simplex("B" + post_face) << "\">Post deleting "
+ << simplex_format(post_face, true) << "</li>\n";
+ ofs << " </ul>\n";
+ ofs << " </p>\n";
+ ofs << " <p id=\"" << id_from_simplex(cc_info.trigger_) << "\">Deleting "
+ << simplex_format(cc_info.trigger_, true) << "</p>\n";
+ }
+ // for (const std::string& fac: cc_info.cofaces_)
+ // ofs << " <li>Checking if " << simplex_format(cc_info.simplex_, true)
+ // << " is a coface of " << simplex_format(fac, true) << "\n";
+ ofs << " </ul>\n";
+ ofs << " </li>\n";
+ if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+ ofs << " <p><span style=\"color:darkorange\">The simplex "
+ << simplex_format(cc_info.simplex_, true) << " is a coface of the simplex "
+ << simplex_format(cc_info.trigger_, true) << "!</span><p>\n";
+ }
+ if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+ ofs << " <p><span style=\"color:green\">Successfully inserted "
+ << simplex_format(cc_info.simplex_, true) << "!</span><p>\n";
+ }
+ }
+ ofs << " </ul>\n";
+ }
+ }
+ ofs << " </div>\n";
+}
+
+void write_visu(std::ofstream& ofs) {
+ ofs << " <div id=\"visu\">\n"
+ << " <h2> Visualization details debug trace </h2>\n";
+ // std::vector<std::map<std::string, std::string> > vs_maps(cc_interior_summary_lists.size());
+ std::map<std::string, std::string> vs_map;
+ for (const auto& sv_pair : simplex_vlist_map) vs_map.emplace(sv_pair.second, sv_pair.first);
+ ofs << " <ul>\n";
+ for (const auto& vs_pair : vs_map) {
+ std::string w_simplex = vs_pair.second.substr(1);
+ bool is_boundary = vs_pair.second[0] == 'B';
+ ofs << " <li><b>" << vs_pair.first << "</b>: " << simplex_format(w_simplex, is_boundary) << "</li>\n";
+ }
+ ofs << " </ul>\n";
+ ofs << " </div>\n";
+}
+
+void write_to_html(std::string file_name) {
+ std::ofstream ofs(file_name + ".html", std::ofstream::out);
+ ofs << "<!DOCTYPE html>\n"
+ << "<html>\n";
+ write_head(ofs);
+ ofs << " <body>\n";
+ write_nav(ofs);
+ ofs << " <h1> Debug traces for " << file_name << " </h1>\n";
+ write_mt(ofs);
+ write_cc(ofs);
+ write_visu(ofs);
+ ofs << " </body>\n";
+ ofs << "</html>\n";
+
+ ofs.close();
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif // DEBUG_TRACES
+#endif // IO_OUTPUT_DEBUG_TRACES_TO_HTML_H_
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
new file mode 100644
index 00000000..f69d8b29
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
@@ -0,0 +1,154 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_OUTPUT_MESHES_TO_MEDIT_H_
+#define IO_OUTPUT_MESHES_TO_MEDIT_H_
+
+#include <gudhi/IO/Mesh_medit.h>
+
+#include <Eigen/Dense>
+
+#include <cstdlib> // for std::size_t
+#include <fstream> // for std::ofstream
+#include <vector>
+#include <type_traits> // for std::enable_if
+#include <tuple> // for std::get
+#include <utility> // for std::make_pair
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+using Vertex_points = Mesh_medit::Vertex_points;
+using Mesh_elements = Mesh_medit::Mesh_elements;
+using Scalar_field_range = Mesh_medit::Scalar_field_range;
+
+template <std::size_t I = 0, typename... Meshes>
+typename std::enable_if<I == sizeof...(Meshes), void>::type fill_meshes(Vertex_points& vertex_points,
+ Mesh_elements& edges, Mesh_elements& triangles,
+ Mesh_elements& tetrahedra,
+ Scalar_field_range& triangles_scalar_range,
+ Scalar_field_range& tetrahedra_scalar_range,
+ std::size_t index, const Meshes&... meshes) {}
+
+template <std::size_t I = 0, typename... Meshes>
+typename std::enable_if<I != sizeof...(Meshes), void>::type fill_meshes(Vertex_points& vertex_points,
+ Mesh_elements& edges, Mesh_elements& triangles,
+ Mesh_elements& tetrahedra,
+ Scalar_field_range& triangles_scalar_range,
+ Scalar_field_range& tetrahedra_scalar_range,
+ std::size_t index, const Meshes&... meshes) {
+ auto mesh = std::get<I>(std::forward_as_tuple(meshes...));
+ for (const auto& v : mesh.vertex_points) vertex_points.push_back(v);
+ for (const auto& e : mesh.edges) {
+ std::vector<std::size_t> edge;
+ for (const auto& v_i : e.first) edge.push_back(v_i + index);
+ edges.emplace_back(edge, e.second);
+ }
+ for (const auto& t : mesh.triangles) {
+ std::vector<std::size_t> triangle;
+ for (const auto& v_i : t.first) triangle.push_back(v_i + index);
+ triangles.emplace_back(triangle, t.second);
+ }
+ for (const auto& t : mesh.tetrahedra) {
+ std::vector<std::size_t> tetrahedron;
+ for (const auto& v_i : t.first) tetrahedron.push_back(v_i + index);
+ tetrahedra.emplace_back(tetrahedron, t.second);
+ }
+ for (const auto& b : mesh.triangles_scalar_range) triangles_scalar_range.push_back(b);
+ for (const auto& b : mesh.tetrahedra_scalar_range) tetrahedra_scalar_range.push_back(b);
+ fill_meshes<I + 1, Meshes...>(vertex_points, edges, triangles, tetrahedra, triangles_scalar_range,
+ tetrahedra_scalar_range, index + mesh.vertex_points.size(), meshes...);
+}
+
+/** \brief Outputs a text file with specified meshes that can be visualized in
+ * <a target="_blank" href="https://www.ljll.math.upmc.fr/frey/software.html">Medit</a>.
+ *
+ * @param[in] amb_d Ambient dimension. Can be 2 or 3.
+ * @param[in] file_name The name of the output file.
+ * @param[in] meshes A pack of meshes to be specified separated by commas.
+ *
+ * @ingroup coxeter_triangulation
+ */
+template <typename... Meshes>
+void output_meshes_to_medit(std::size_t amb_d, std::string file_name, const Meshes&... meshes) {
+ Vertex_points vertex_points;
+ Mesh_elements edges, triangles, tetrahedra;
+ Scalar_field_range triangles_scalar_range, tetrahedra_scalar_range;
+ fill_meshes(vertex_points, edges, triangles, tetrahedra, triangles_scalar_range, tetrahedra_scalar_range, 0,
+ meshes...);
+
+ std::ofstream ofs(file_name + ".mesh", std::ofstream::out);
+ std::ofstream ofs_bb(file_name + ".bb", std::ofstream::out);
+
+ if (amb_d == 2) {
+ ofs << "MeshVersionFormatted 1\nDimension 2\n";
+ ofs_bb << "2 1 ";
+ ofs << "Vertices\n" << vertex_points.size() << "\n";
+ for (auto p : vertex_points) {
+ ofs << p[0] << " " << p[1] << " 2\n";
+ }
+ ofs << "Edges " << edges.size() << "\n";
+ for (auto e : edges) {
+ for (auto v : e.first) ofs << v << " ";
+ ofs << e.second << std::endl;
+ }
+ ofs << "Triangles " << triangles.size() << "\n";
+ for (auto s : triangles) {
+ for (auto v : s.first) {
+ ofs << v << " ";
+ }
+ ofs << s.second << std::endl;
+ }
+
+ ofs_bb << triangles_scalar_range.size() << " 1\n";
+ for (auto& b : triangles_scalar_range) ofs_bb << b << "\n";
+
+ } else {
+ ofs << "MeshVersionFormatted 1\nDimension 3\n";
+ ofs_bb << "3 1 ";
+ ofs << "Vertices\n" << vertex_points.size() << "\n";
+ for (auto p : vertex_points) {
+ ofs << p[0] << " " << p[1] << " " << p[2] << " 2\n";
+ }
+ ofs << "Edges " << edges.size() << "\n";
+ for (auto e : edges) {
+ for (auto v : e.first) ofs << v << " ";
+ ofs << e.second << std::endl;
+ }
+ ofs << "Triangles " << triangles.size() << "\n";
+ for (auto s : triangles) {
+ for (auto v : s.first) {
+ ofs << v << " ";
+ }
+ ofs << s.second << std::endl;
+ }
+ ofs << "Tetrahedra " << tetrahedra.size() << "\n";
+ for (auto s : tetrahedra) {
+ for (auto v : s.first) {
+ ofs << v << " ";
+ }
+ ofs << s.second << std::endl;
+ }
+
+ ofs_bb << triangles_scalar_range.size() + tetrahedra_scalar_range.size() << " 1\n";
+ for (auto& b : triangles_scalar_range) ofs_bb << b << "\n";
+ for (auto& b : tetrahedra_scalar_range) ofs_bb << b << "\n";
+ }
+
+ ofs.close();
+ ofs_bb.close();
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
new file mode 100644
index 00000000..277f8b6c
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
@@ -0,0 +1,261 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
+#define IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
+
+#include <Eigen/Dense>
+
+#include <gudhi/Permutahedral_representation/face_from_indices.h>
+#include <gudhi/Functions/Constant_function.h>
+#include <gudhi/Functions/PL_approximation.h>
+#include <gudhi/Coxeter_triangulation/Query_result.h>
+#include <gudhi/Debug_utils.h> // for GUDHI_CHECK
+
+#include <vector>
+#include <limits> // for std::numeric_limits<>
+#include <cmath> // for std::fabs
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Implicit_manifold_intersection_oracle
+ * \brief An oracle that supports the intersection query on an implicit manifold.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ * \tparam Domain_function_ The domain function template parameter. Should be a model of
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_, class Domain_function_ = Constant_function>
+class Implicit_manifold_intersection_oracle {
+ /* Computes the affine coordinates of the intersection point of the implicit manifold
+ * and the affine hull of the simplex. */
+ template <class Simplex_handle, class Triangulation>
+ Eigen::VectorXd compute_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const {
+ std::size_t cod_d = this->cod_d();
+ Eigen::MatrixXd matrix(cod_d + 1, cod_d + 1);
+ for (std::size_t i = 0; i < cod_d + 1; ++i) matrix(0, i) = 1;
+ std::size_t j = 0;
+ for (auto v : simplex.vertex_range()) {
+ Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
+ for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1);
+ j++;
+ }
+ Eigen::VectorXd z(cod_d + 1);
+ z(0) = 1;
+ for (std::size_t i = 1; i < cod_d + 1; ++i) z(i) = 0;
+ Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+ if (!z.isApprox(matrix*lambda)) {
+ // NaN non valid results
+ for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
+ std::numeric_limits<double>::quiet_NaN();
+ }
+ return lambda;
+ }
+
+ /* Computes the affine coordinates of the intersection point of the boundary
+ * of the implicit manifold and the affine hull of the simplex. */
+ template <class Simplex_handle, class Triangulation>
+ Eigen::VectorXd compute_boundary_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const {
+ std::size_t cod_d = this->cod_d();
+ Eigen::MatrixXd matrix(cod_d + 2, cod_d + 2);
+ for (std::size_t i = 0; i < cod_d + 2; ++i) matrix(0, i) = 1;
+ std::size_t j = 0;
+ for (auto v : simplex.vertex_range()) {
+ Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
+ for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1);
+ Eigen::VectorXd bv_coords = domain_fun_(triangulation.cartesian_coordinates(v));
+ matrix(cod_d + 1, j) = bv_coords(0);
+ j++;
+ }
+ Eigen::VectorXd z(cod_d + 2);
+ z(0) = 1;
+ for (std::size_t i = 1; i < cod_d + 2; ++i) z(i) = 0;
+ Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+ if (!z.isApprox(matrix*lambda)) {
+ // NaN non valid results
+ for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
+ std::numeric_limits<double>::quiet_NaN();
+ }
+ return lambda;
+ }
+
+ /* Computes the intersection result for a given simplex in a triangulation. */
+ template <class Simplex_handle, class Triangulation>
+ Query_result<Simplex_handle> intersection_result(const Eigen::VectorXd& lambda, const Simplex_handle& simplex,
+ const Triangulation& triangulation) const {
+ using QR = Query_result<Simplex_handle>;
+ std::size_t amb_d = triangulation.dimension();
+ std::size_t cod_d = simplex.dimension();
+ for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) {
+ if (std::isnan(lambda(i))) return QR({Eigen::VectorXd(), false});
+ GUDHI_CHECK((std::fabs(lambda(i) - 1.) > std::numeric_limits<double>::epsilon() &&
+ std::fabs(lambda(i) - 0.) > std::numeric_limits<double>::epsilon()),
+ std::invalid_argument("A vertex of the triangulation lies exactly on the manifold"));
+ if (lambda(i) < 0. || lambda(i) > 1.) return QR({Eigen::VectorXd(), false});
+ }
+ Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d);
+ auto v_range = simplex.vertex_range();
+ auto v_it = v_range.begin();
+ for (std::size_t i = 0; i < cod_d + 1 && v_it != v_range.end(); ++v_it, ++i) {
+ Eigen::VectorXd v_coords = triangulation.cartesian_coordinates(*v_it);
+ for (std::size_t j = 0; j < amb_d; ++j) vertex_matrix(i, j) = v_coords(j);
+ }
+ Eigen::VectorXd intersection = lambda.transpose() * vertex_matrix;
+ return QR({intersection, true});
+ }
+
+ public:
+ /** \brief Ambient dimension of the implicit manifold. */
+ std::size_t amb_d() const { return fun_.amb_d(); }
+
+ /** \brief Codimension of the implicit manifold. */
+ std::size_t cod_d() const { return fun_.cod_d(); }
+
+ /** \brief The seed point of the implicit manifold. */
+ Eigen::VectorXd seed() const { return fun_.seed(); }
+
+ /** \brief Intersection query with the relative interior of the manifold.
+ *
+ * \details The returned structure Query_result contains the boolean value
+ * that is true only if the intersection point of the query simplex and
+ * the relative interior of the manifold exists, the intersection point
+ * and the face of the query simplex that contains
+ * the intersection point.
+ *
+ * \tparam Simplex_handle The class of the query simplex.
+ * Needs to be a model of the concept SimplexInCoxeterTriangulation.
+ * \tparam Triangulation The class of the triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ *
+ * @param[in] simplex The query simplex. The dimension of the simplex
+ * should be the same as the codimension of the manifold
+ * (the codomain dimension of the function).
+ * @param[in] triangulation The ambient triangulation. The dimension of
+ * the triangulation should be the same as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Simplex_handle, class Triangulation>
+ Query_result<Simplex_handle> intersects(const Simplex_handle& simplex, const Triangulation& triangulation) const {
+ Eigen::VectorXd lambda = compute_lambda(simplex, triangulation);
+ return intersection_result(lambda, simplex, triangulation);
+ }
+
+ /** \brief Intersection query with the boundary of the manifold.
+ *
+ * \details The returned structure Query_result contains the boolean value
+ * that is true only if the intersection point of the query simplex and
+ * the boundary of the manifold exists, the intersection point
+ * and the face of the query simplex that contains
+ * the intersection point.
+ *
+ * \tparam Simplex_handle The class of the query simplex.
+ * Needs to be a model of the concept SimplexInCoxeterTriangulation.
+ * \tparam Triangulation The class of the triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ *
+ * @param[in] simplex The query simplex. The dimension of the simplex
+ * should be the same as the codimension of the boundary of the manifold
+ * (the codomain dimension of the function + 1).
+ * @param[in] triangulation The ambient triangulation. The dimension of
+ * the triangulation should be the same as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Simplex_handle, class Triangulation>
+ Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex,
+ const Triangulation& triangulation) const {
+ //std::cout << "intersects_boundary" << std::endl;
+ Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation);
+ return intersection_result(lambda, simplex, triangulation);
+ }
+
+ /** \brief Returns true if the input point lies inside the piecewise-linear
+ * domain induced by the given ambient triangulation that defines the relative
+ * interior of the piecewise-linear approximation of the manifold.
+ *
+ * @param p The input point. Needs to have the same dimension as the ambient
+ * dimension of the manifold (the domain dimension of the function).
+ * @param triangulation The ambient triangulation. Needs to have the same
+ * dimension as the ambient dimension of the manifold
+ * (the domain dimension of the function).
+ */
+ template <class Triangulation>
+ bool lies_in_domain(const Eigen::VectorXd& p, const Triangulation& triangulation) const {
+ Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p);
+ return pl_p(0) < 0;
+ }
+
+ /** \brief Returns the function that defines the interior of the manifold */
+ const Function_& function() const { return fun_; }
+
+ /** \brief Constructs an intersection oracle for an implicit manifold potentially
+ * with boundary from given function and domain.
+ *
+ * @param function The input function that represents the implicit manifold
+ * before the restriction with the domain.
+ * @param domain_function The input domain function that can be used to define an implicit
+ * manifold with boundary.
+ */
+ Implicit_manifold_intersection_oracle(const Function_& function, const Domain_function_& domain_function)
+ : fun_(function), domain_fun_(domain_function) {}
+
+ /** \brief Constructs an intersection oracle for an implicit manifold
+ * without boundary from a given function.
+ *
+ * \details To use this constructor, the template Domain_function_ needs to be left
+ * at its default value (Gudhi::coxeter_triangulation::Constant_function).
+ *
+ * @param function The input function that represents the implicit manifold
+ * without boundary.
+ */
+ Implicit_manifold_intersection_oracle(const Function_& function)
+ : fun_(function), domain_fun_(function.amb_d(), 1, Eigen::VectorXd::Constant(1, -1)) {}
+
+ private:
+ Function_ fun_;
+ Domain_function_ domain_fun_;
+};
+
+/** \brief Static constructor of an intersection oracle from a function with a domain.
+ *
+ * @param function The input function that represents the implicit manifold
+ * before the restriction with the domain.
+ * @param domain_function The input domain function that can be used to define an implicit
+ * manifold with boundary.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_, class Domain_function_>
+Implicit_manifold_intersection_oracle<Function_, Domain_function_> make_oracle(
+ const Function_& function, const Domain_function_& domain_function) {
+ return Implicit_manifold_intersection_oracle<Function_, Domain_function_>(function, domain_function);
+}
+
+/** \brief Static constructor of an intersection oracle from a function without a domain.
+ *
+ * @param function The input function that represents the implicit manifold
+ * without boundary.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Implicit_manifold_intersection_oracle<Function_> make_oracle(const Function_& function) {
+ return Implicit_manifold_intersection_oracle<Function_>(function);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h
new file mode 100644
index 00000000..d61bbed7
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h
@@ -0,0 +1,270 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef MANIFOLD_TRACING_H_
+#define MANIFOLD_TRACING_H_
+
+#include <gudhi/IO/output_debug_traces_to_html.h> // for DEBUG_TRACES
+#include <gudhi/Coxeter_triangulation/Query_result.h>
+
+#include <boost/functional/hash.hpp>
+
+#include <Eigen/Dense>
+
+#include <queue>
+#include <unordered_map>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \ingroup coxeter_triangulation
+ */
+
+/** \class Manifold_tracing
+ * \brief A class that assembles methods for manifold tracing algorithm.
+ *
+ * \tparam Triangulation_ The type of the ambient triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ */
+template <class Triangulation_>
+class Manifold_tracing {
+ public:
+ using Simplex_handle = typename Triangulation_::Simplex_handle;
+
+ struct Simplex_hash {
+ typedef Simplex_handle argument_type;
+ typedef std::size_t result_type;
+ result_type operator()(const argument_type& s) const noexcept {
+ return boost::hash<typename Simplex_handle::Vertex>()(s.vertex());
+ }
+ };
+
+ public:
+ /** \brief Type of the output simplex map with keys of type Triangulation_::Simplex_handle
+ * and values of type Eigen::VectorXd.
+ * This type should be used for the output in the method manifold_tracing_algorithm.
+ */
+ typedef std::unordered_map<Simplex_handle, Eigen::VectorXd, Simplex_hash> Out_simplex_map;
+
+ /**
+ * \brief Computes the set of k-simplices that intersect
+ * a boundaryless implicit manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ * Needs to be a model of the concept IntersectionOracle.
+ *
+ * \param[in] seed_points The range of points on the manifold from which
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] out_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the input manifold and the mapped values
+ * are the intersection points.
+ */
+ template <class Point_range, class Intersection_oracle>
+ void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation,
+ const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) {
+ std::size_t cod_d = oracle.cod_d();
+ std::queue<Simplex_handle> queue;
+
+ for (const auto& p : seed_points) {
+ Simplex_handle full_simplex = triangulation.locate_point(p);
+ for (Simplex_handle face : full_simplex.face_range(cod_d)) {
+ Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
+ if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) {
+#ifdef DEBUG_TRACES
+ mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif
+ queue.emplace(face);
+ break;
+ }
+ }
+ }
+
+ while (!queue.empty()) {
+ Simplex_handle s = queue.front();
+ queue.pop();
+ for (auto cof : s.coface_range(cod_d + 1)) {
+ for (auto face : cof.face_range(cod_d)) {
+ Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
+ if (qr.success && out_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
+ }
+ }
+ }
+ }
+
+ /**
+ * \brief Computes the set of k-simplices that intersect
+ * the dimensional manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ * Needs to be a model of the concept IntersectionOracle.
+ *
+ * \param[in] seed_points The range of points on the manifold from which
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] interior_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the relative interior of the input manifold
+ * and the mapped values are the intersection points.
+ * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the boundary of the input manifold
+ * and the mapped values are the intersection points.
+ */
+ template <class Point_range, class Intersection_oracle>
+ void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation_& triangulation,
+ const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map,
+ Out_simplex_map& boundary_simplex_map) {
+ std::size_t cod_d = oracle.cod_d();
+ std::queue<Simplex_handle> queue;
+
+ for (const auto& p : seed_points) {
+ Simplex_handle full_simplex = triangulation.locate_point(p);
+ for (Simplex_handle face : full_simplex.face_range(cod_d)) {
+ auto qr = oracle.intersects(face, triangulation);
+#ifdef DEBUG_TRACES
+ mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif
+ if (qr.success) {
+ if (oracle.lies_in_domain(qr.intersection, triangulation)) {
+ if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
+ } else {
+ for (Simplex_handle cof : face.coface_range(cod_d + 1)) {
+ auto qrb = oracle.intersects_boundary(cof, triangulation);
+#ifdef DEBUG_TRACES
+ mt_seed_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
+#endif
+ if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection);
+ }
+ }
+ // break;
+ }
+ }
+ }
+
+ while (!queue.empty()) {
+ Simplex_handle s = queue.front();
+ queue.pop();
+ for (auto cof : s.coface_range(cod_d + 1)) {
+ for (auto face : cof.face_range(cod_d)) {
+ auto qr = oracle.intersects(face, triangulation);
+#ifdef DEBUG_TRACES
+ mt_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif
+ if (qr.success) {
+ if (oracle.lies_in_domain(qr.intersection, triangulation)) {
+ if (interior_simplex_map.emplace(face, qr.intersection).second) queue.emplace(face);
+ } else {
+ auto qrb = oracle.intersects_boundary(cof, triangulation);
+#ifdef DEBUG_TRACES
+ mt_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
+#endif
+ if (qrb.success) boundary_simplex_map.emplace(cof, qrb.intersection);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /** \brief Empty constructor */
+ Manifold_tracing() {}
+};
+
+/**
+ * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm
+ * that computes the set of k-simplices that intersect
+ * a boundaryless implicit manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ * \tparam Triangulation_ The type of the ambient triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ * Needs to be a model of the concept IntersectionOracle.
+ * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map.
+ *
+ * \param[in] seed_points The range of points on the manifold from which
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] out_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the input manifold and the mapped values
+ * are the intersection points.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map>
+void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation,
+ const Intersection_oracle& oracle, Out_simplex_map& out_simplex_map) {
+ Manifold_tracing<Triangulation> mt;
+ mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map);
+}
+
+/**
+ * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm
+ * the dimensional manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ * \tparam Triangulation_ The type of the ambient triangulation.
+ * Needs to be a model of the concept TriangulationForManifoldTracing.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ * Needs to be a model of the concept IntersectionOracle.
+ * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map.
+ *
+ * \param[in] seed_points The range of points on the manifold from which
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] interior_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the relative interior of the input manifold
+ * and the mapped values are the intersection points.
+ * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the boundary of the input manifold
+ * and the mapped values are the intersection points.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Point_range, class Triangulation, class Intersection_oracle, class Out_simplex_map>
+void manifold_tracing_algorithm(const Point_range& seed_points, const Triangulation& triangulation,
+ const Intersection_oracle& oracle, Out_simplex_map& interior_simplex_map,
+ Out_simplex_map& boundary_simplex_map) {
+ Manifold_tracing<Triangulation> mt;
+ mt.manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h
new file mode 100644
index 00000000..76438c91
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h
@@ -0,0 +1,216 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_H_
+#define PERMUTAHEDRAL_REPRESENTATION_H_
+
+#include <gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h>
+
+#include <utility> // for std::make_pair
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Permutahedral_representation
+ * \brief A class that stores the permutahedral representation of a simplex
+ * in a Coxeter triangulation or a Freudenthal-Kuhn triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \details The data structure is a record consisting of a range that
+ * represents the vertex and a range that represents the ordered set
+ * partition, both of which identify the simplex in the triangulation.
+ *
+ * \tparam Vertex_ needs to be a random-access range.
+ * \tparam Ordered_set_partition_ needs to be a a random-access range that consists of
+ * random-access ranges.
+ */
+template <class Vertex_, class Ordered_set_partition_>
+class Permutahedral_representation {
+ typedef Permutahedral_representation<Vertex_, Ordered_set_partition_> Self;
+
+ public:
+ /** \brief Type of the vertex. */
+ typedef Vertex_ Vertex;
+
+ /** \brief Type of the ordered partition. */
+ typedef Ordered_set_partition_ OrderedSetPartition;
+
+ /** \brief Permutahedral_representation constructor from a vertex and an ordered set partition.
+ *
+ * @param[in] vertex Vertex.
+ * @param[in] partition Ordered set partition.
+ *
+ * \details If the size of vertex is d, the ranges in partition must consist
+ * of the integers 0,...,d without repetition or collision between the ranges.
+ */
+ Permutahedral_representation(const Vertex& vertex, const OrderedSetPartition& partition)
+ : vertex_(vertex), partition_(partition) {}
+
+ /** \brief Constructor for an empty permutahedral representation that does not correspond
+ * to any simplex.
+ */
+ Permutahedral_representation() {}
+
+ /** \brief Dimension of the simplex. */
+ std::size_t dimension() const { return partition_.size() - 1; }
+
+ /** \brief Lexicographically-minimal vertex. */
+ Vertex& vertex() { return vertex_; }
+
+ /** \brief Lexicographically-minimal vertex. */
+ const Vertex& vertex() const { return vertex_; }
+
+ /** \brief Ordered set partition. */
+ OrderedSetPartition& partition() { return partition_; }
+
+ /** \brief Identifying vertex. */
+ const OrderedSetPartition& partition() const { return partition_; }
+
+ /** \brief Equality operator.
+ * Returns true if an only if both vertex and the ordered set partition coincide.
+ */
+ bool operator==(const Permutahedral_representation& other) const {
+ if (dimension() != other.dimension()) return false;
+ if (vertex_ != other.vertex_) return false;
+ for (std::size_t k = 0; k < partition_.size(); ++k)
+ if (partition_[k] != other.partition_[k]) return false;
+ return true;
+ }
+
+ /** \brief Inequality operator.
+ * Returns true if an only if either vertex or the ordered set partition are different.
+ */
+ bool operator!=(const Permutahedral_representation& other) const { return !(*this == other); }
+
+ typedef Gudhi::coxeter_triangulation::Vertex_iterator<Self> Vertex_iterator;
+ typedef boost::iterator_range<Vertex_iterator> Vertex_range;
+ /** \brief Returns a range of vertices of the simplex.
+ * The type of vertices is Vertex.
+ */
+ Vertex_range vertex_range() const { return Vertex_range(Vertex_iterator(*this), Vertex_iterator()); }
+
+ typedef Gudhi::coxeter_triangulation::Face_iterator<Self> Face_iterator;
+ typedef boost::iterator_range<Face_iterator> Face_range;
+ /** \brief Returns a range of permutahedral representations of faces of the simplex.
+ * @param[in] value_dim The dimension of the faces. Must be between 0 and the dimension of the simplex.
+ */
+ Face_range face_range(std::size_t value_dim) const {
+ return Face_range(Face_iterator(*this, value_dim), Face_iterator());
+ }
+
+ /** \brief Returns a range of permutahedral representations of facets of the simplex.
+ * The dimension of the simplex must be strictly positive.
+ */
+ Face_range facet_range() const { return Face_range(Face_iterator(*this, dimension() - 1), Face_iterator()); }
+
+ typedef Gudhi::coxeter_triangulation::Coface_iterator<Self> Coface_iterator;
+ typedef boost::iterator_range<Coface_iterator> Coface_range;
+ /** \brief Returns a range of permutahedral representations of cofaces of the simplex.
+ * @param[in] value_dim The dimension of the cofaces. Must be between the dimension of the simplex and the ambient
+ * dimension (the size of the vertex).
+ */
+ Coface_range coface_range(std::size_t value_dim) const {
+ return Coface_range(Coface_iterator(*this, value_dim), Coface_iterator());
+ }
+
+ /** \brief Returns a range of permutahedral representations of cofacets of the simplex.
+ * The dimension of the simplex must be strictly different from the ambient dimension (the size of the vertex).
+ */
+ Coface_range cofacet_range() const {
+ return Coface_range(Coface_iterator(*this, dimension() + 1), Coface_iterator());
+ }
+
+ /** \brief Returns true, if the simplex is a face of other simplex.
+ *
+ * @param[in] other A simplex that is potential a coface of the current simplex.
+ */
+ bool is_face_of(const Permutahedral_representation& other) const {
+ using Part = typename OrderedSetPartition::value_type;
+
+ if (other.dimension() < dimension()) return false;
+ if (other.vertex_.size() != vertex_.size())
+ std::cerr << "Error: Permutahedral_representation::is_face_of: incompatible ambient dimensions.\n";
+
+ Vertex v_self = vertex_, v_other = other.vertex_;
+ auto self_partition_it = partition_.begin();
+ auto other_partition_it = other.partition_.begin();
+ while (self_partition_it != partition_.end()) {
+ while (other_partition_it != other.partition_.end() && v_self != v_other) {
+ const Part& other_part = *other_partition_it++;
+ if (other_partition_it == other.partition_.end()) return false;
+ for (const auto& k : other_part) v_other[k]++;
+ }
+ if (other_partition_it == other.partition_.end()) return false;
+ const Part& self_part = *self_partition_it++;
+ if (self_partition_it == partition_.end()) return true;
+ for (const auto& k : self_part) v_self[k]++;
+ }
+ return true;
+ }
+
+ private:
+ Vertex vertex_;
+ OrderedSetPartition partition_;
+};
+
+/** \brief Print a permutahedral representation to a stream.
+ * \ingroup coxeter_triangulation
+ *
+ * @param[in] os The output stream.
+ * @param[in] simplex A simplex represented by its permutahedral representation.
+ */
+template <class Vertex, class OrderedSetPartition>
+std::ostream& operator<<(std::ostream& os, const Permutahedral_representation<Vertex, OrderedSetPartition>& simplex) {
+ // vertex part
+ os << "(";
+ if (simplex.vertex().empty()) {
+ os << ")";
+ return os;
+ }
+ auto v_it = simplex.vertex().begin();
+ os << *v_it++;
+ for (; v_it != simplex.vertex().end(); ++v_it) os << ", " << *v_it;
+ os << ")";
+
+ // ordered partition part
+ using Part = typename OrderedSetPartition::value_type;
+ auto print_part = [&os](const Part& p) {
+ os << "{";
+ if (p.empty()) {
+ os << "}";
+ }
+ auto p_it = p.begin();
+ os << *p_it++;
+ for (; p_it != p.end(); ++p_it) os << ", " << *p_it;
+ os << "}";
+ };
+ os << " [";
+ if (simplex.partition().empty()) {
+ os << "]";
+ return os;
+ }
+ auto o_it = simplex.partition().begin();
+ print_part(*o_it++);
+ for (; o_it != simplex.partition().end(); ++o_it) {
+ os << ", ";
+ print_part(*o_it);
+ }
+ os << "]";
+ return os;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif // PERMUTAHEDRAL_REPRESENTATION_H_
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h
new file mode 100644
index 00000000..5f382e31
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h
@@ -0,0 +1,83 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_
+
+#include <vector>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate combinations of
+ * k elements in a set of n elements.
+ * Based on the algorithm by Mifsud.
+ */
+class Combination_iterator
+ : public boost::iterator_facade<Combination_iterator, std::vector<uint> const, boost::forward_traversal_tag> {
+ typedef std::vector<uint> value_t;
+
+ protected:
+ friend class boost::iterator_core_access;
+
+ bool equal(Combination_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void increment() {
+ if (value_[0] == n_ - k_) {
+ is_end_ = true;
+ return;
+ }
+ uint j = k_ - 1;
+ if (value_[j] < n_ - 1) {
+ value_[j]++;
+ return;
+ }
+ for (; j > 0; --j)
+ if (value_[j - 1] < n_ - k_ + j - 1) {
+ value_[j - 1]++;
+ for (uint s = j; s < k_; s++) value_[s] = value_[j - 1] + s - (j - 1);
+ return;
+ }
+ }
+
+ public:
+ Combination_iterator(const uint& n, const uint& k) : value_(k), is_end_(n == 0), n_(n), k_(k) {
+ for (uint i = 0; i < k; ++i) value_[i] = i;
+ }
+
+ // Used for the creating an end iterator
+ Combination_iterator() : is_end_(true), n_(0), k_(0) {}
+
+ void reinitialize() {
+ if (n_ > 0) {
+ is_end_ = false;
+ for (uint i = 0; i < n_; ++i) value_[i] = i;
+ }
+ }
+
+ private:
+ value_t value_; // the dereference value
+ bool is_end_; // is true when the current permutation is the final one
+
+ uint n_;
+ uint k_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h
new file mode 100644
index 00000000..3ee73754
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h
@@ -0,0 +1,114 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_
+
+#include <vector>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate combinations of
+ * k elements in a set of n elements.
+ * Based on the algorithm by Mifsud.
+ */
+class Integer_combination_iterator
+ : public boost::iterator_facade<Integer_combination_iterator, std::vector<uint> const,
+ boost::forward_traversal_tag> {
+ using value_t = std::vector<uint>;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ bool equal(Integer_combination_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void increment() {
+ uint j1 = 0;
+ uint s = 0;
+ while (value_[j1] == 0 && j1 < k_) j1++;
+ uint j2 = j1 + 1;
+ while (value_[j2] == bounds_[j2]) {
+ if (bounds_[j2] != 0) {
+ s += value_[j1];
+ value_[j1] = 0;
+ j1 = j2;
+ }
+ j2++;
+ }
+ if (j2 >= k_) {
+ is_end_ = true;
+ return;
+ }
+ s += value_[j1] - 1;
+ value_[j1] = 0;
+ value_[j2]++;
+ uint i = 0;
+ while (s >= bounds_[i]) {
+ value_[i] = bounds_[i];
+ s -= bounds_[i];
+ i++;
+ }
+ value_[i++] = s;
+ }
+
+ public:
+ template <class Bound_range>
+ Integer_combination_iterator(const uint& n, const uint& k, const Bound_range& bounds)
+ : value_(k + 2), is_end_(n == 0 || k == 0), 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) {}
+
+ private:
+ 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<uint> bounds_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h
new file mode 100644
index 00000000..866079fa
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h
@@ -0,0 +1,93 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_
+
+#include <vector>
+#include <limits>
+
+#include <gudhi/Permutahedral_representation/Permutation_iterator.h>
+#include <gudhi/Permutahedral_representation/Set_partition_iterator.h>
+
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that represents an ordered set partition of a set {0,...,n-1} in k parts as
+ * a pair of an unordered set partition given in lexicographic order and
+ * a permutation of the parts.
+ */
+struct Ordered_set_partition {
+ Set_partition_iterator s_it_;
+ Permutation_iterator p_it_;
+
+ // Ordered_set_partition(const Set_partition_iterator& s_it, const Permutation_iterator& p_it)
+ // : s_it_(s_it), p_it_(p_it) {}
+
+ const std::vector<uint> operator[](const uint& i) const { return (*s_it_)[(*p_it_)[i]]; }
+
+ std::size_t size() const { return s_it_->size(); }
+};
+
+/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts.
+ *
+ */
+class Ordered_set_partition_iterator
+ : public boost::iterator_facade<Ordered_set_partition_iterator, Ordered_set_partition const,
+ boost::forward_traversal_tag> {
+ using value_t = Ordered_set_partition;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ bool equal(Ordered_set_partition_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void increment() {
+ if (++value_.p_it_ == p_end_) {
+ if (++value_.s_it_ == s_end_) {
+ is_end_ = true;
+ return;
+ } else
+ value_.p_it_.reinitialize();
+ }
+ }
+
+ public:
+ Ordered_set_partition_iterator(const uint& n, const uint& k)
+ : value_({Set_partition_iterator(n, k), Permutation_iterator(k)}), is_end_(n == 0) {}
+
+ // Used for the creating an end iterator
+ Ordered_set_partition_iterator() : is_end_(true) {}
+
+ void reinitialize() {
+ is_end_ = false;
+ value_.p_it_.reinitialize();
+ value_.s_it_.reinitialize();
+ }
+
+ private:
+ Set_partition_iterator s_end_; // Set partition iterator and the corresponding end iterator
+ Permutation_iterator p_end_; // Permutation iterator and the corresponding end iterator
+ value_t value_; // the dereference value
+ bool is_end_; // is true when the current permutation is the final one
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h
new file mode 100644
index 00000000..db145741
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h
@@ -0,0 +1,254 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_
+#define PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_
+
+#include <gudhi/Permutahedral_representation/Size_range.h>
+#include <gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h>
+#include <gudhi/Permutahedral_representation/Integer_combination_iterator.h>
+#include <gudhi/Permutahedral_representation/Combination_iterator.h>
+#include <gudhi/Permutahedral_representation/face_from_indices.h>
+#include <boost/iterator/iterator_facade.hpp>
+
+#include <vector>
+#include <iostream>
+#include <algorithm> // for std::find
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* \addtogroup coxeter_triangulation
+ * Iterator types for Permutahedral_representation
+ * @{
+ */
+
+/* \brief Iterator over the vertices of a simplex
+ * represented by its permutahedral representation.
+ *
+ * Forward iterator, 'value_type' is Permutahedral_representation::Vertex.*/
+template <class Permutahedral_representation>
+class Vertex_iterator
+ : public boost::iterator_facade<Vertex_iterator<Permutahedral_representation>,
+ typename Permutahedral_representation::Vertex const, boost::forward_traversal_tag> {
+ private:
+ friend class boost::iterator_core_access;
+
+ using Vertex = typename Permutahedral_representation::Vertex;
+ using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+ using value_t = Vertex;
+
+ bool equal(Vertex_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void update_value() {
+ std::size_t d = value_.size();
+ for (auto i : *o_it_)
+ if (i != d)
+ value_[i]++;
+ else
+ for (std::size_t j = 0; j < d; j++) value_[j]--;
+ }
+
+ void increment() {
+ if (is_end_) return;
+ update_value();
+ if (++o_it_ == o_end_) is_end_ = true;
+ }
+
+ public:
+ Vertex_iterator(const Permutahedral_representation& simplex)
+ : o_it_(simplex.partition().begin()),
+ o_end_(simplex.partition().end()),
+ value_(simplex.vertex()),
+ is_end_(o_it_ == o_end_) {}
+
+ Vertex_iterator() : is_end_(true) {}
+
+ private:
+ typename Ordered_partition::const_iterator o_it_, o_end_;
+ value_t value_;
+ bool is_end_;
+
+}; // Vertex_iterator
+
+/*---------------------------------------------------------------------------*/
+/* \brief Iterator over the k-faces of a simplex
+ * given by its permutahedral representation.
+ *
+ * Forward iterator, value_type is Permutahedral_representation. */
+template <class Permutahedral_representation>
+class Face_iterator : public boost::iterator_facade<Face_iterator<Permutahedral_representation>,
+ Permutahedral_representation const, boost::forward_traversal_tag> {
+ using value_t = Permutahedral_representation;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ using Vertex = typename Permutahedral_representation::Vertex;
+ using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+ bool equal(Face_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void increment() {
+ if (++c_it_ == c_end_) {
+ is_end_ = true;
+ return;
+ }
+ update_value();
+ }
+
+ void update_value() {
+ // Combination *c_it_ is supposed to be sorted in increasing order
+ value_ = face_from_indices<Permutahedral_representation>(simplex_, *c_it_);
+ }
+
+ public:
+ Face_iterator(const Permutahedral_representation& simplex, const uint& k)
+ : simplex_(simplex),
+ k_(k),
+ l_(simplex.dimension()),
+ c_it_(l_ + 1, k_ + 1),
+ is_end_(k_ > l_),
+ value_({Vertex(simplex.vertex().size()), Ordered_partition(k + 1)}) {
+ update_value();
+ }
+
+ // Used for the creating an end iterator
+ Face_iterator() : is_end_(true) {}
+
+ private:
+ Permutahedral_representation simplex_; // Input simplex
+ uint k_;
+ uint l_; // Dimension of the input simplex
+ Combination_iterator c_it_, c_end_; // indicates the vertices in the current face
+
+ bool is_end_; // is true when the current permutation is the final one
+ value_t value_; // the dereference value
+
+}; // Face_iterator
+
+/*---------------------------------------------------------------------------*/
+/* \brief Iterator over the k-cofaces of a simplex
+ * given by its permutahedral representation.
+ *
+ * Forward iterator, value_type is Permutahedral_representation. */
+template <class Permutahedral_representation>
+class Coface_iterator
+ : public boost::iterator_facade<Coface_iterator<Permutahedral_representation>, Permutahedral_representation const,
+ boost::forward_traversal_tag> {
+ using value_t = Permutahedral_representation;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ using Vertex = typename Permutahedral_representation::Vertex;
+ using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+ bool equal(Coface_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void increment() {
+ uint i = 0;
+ for (; i < k_ + 1; i++) {
+ if (++(o_its_[i]) != o_end_) break;
+ }
+ if (i == k_ + 1) {
+ if (++i_it_ == i_end_) {
+ is_end_ = true;
+ return;
+ }
+ o_its_.clear();
+ for (uint j = 0; j < k_ + 1; j++)
+ o_its_.emplace_back(simplex_.partition()[j].size(), (*i_it_)[j] + 1);
+ } else
+ for (uint j = 0; j < i; j++) o_its_[j].reinitialize();
+ update_value();
+ }
+
+ void update_value() {
+ value_.vertex() = simplex_.vertex();
+ for (auto& p : value_.partition()) p.clear();
+ uint u_ = 0; // the part in o_its_[k_] that contains t_
+ for (; u_ <= (*i_it_)[k_]; u_++) {
+ auto range = (*o_its_[k_])[u_];
+ if (std::find(range.begin(), range.end(), t_) != range.end()) break;
+ }
+ uint i = 0;
+ for (uint j = u_ + 1; j <= (*i_it_)[k_]; j++, i++)
+ for (uint b : (*o_its_[k_])[j]) {
+ uint c = simplex_.partition()[k_][b];
+ value_.partition()[i].push_back(c);
+ value_.vertex()[c]--;
+ }
+ for (uint h = 0; h < k_; h++)
+ for (uint j = 0; j <= (*i_it_)[h]; j++, i++) {
+ for (uint b : (*o_its_[h])[j]) value_.partition()[i].push_back(simplex_.partition()[h][b]);
+ }
+ for (uint j = 0; j <= u_; j++, i++)
+ for (uint b : (*o_its_[k_])[j]) value_.partition()[i].push_back(simplex_.partition()[k_][b]);
+ // sort the values in each part (probably not needed)
+ for (auto& part : value_.partition()) std::sort(part.begin(), part.end());
+ }
+
+ public:
+ Coface_iterator(const Permutahedral_representation& simplex, const uint& l)
+ : simplex_(simplex),
+ d_(simplex.vertex().size()),
+ l_(l),
+ k_(simplex.dimension()),
+ i_it_(l_ - k_, k_ + 1, Size_range<Ordered_partition>(simplex.partition())),
+ is_end_(k_ > l_),
+ value_({Vertex(d_), Ordered_partition(l_ + 1)}) {
+ uint j = 0;
+ for (; j < simplex_.partition()[k_].size(); j++)
+ if (simplex_.partition()[k_][j] == d_) {
+ t_ = j;
+ break;
+ }
+ if (j == simplex_.partition()[k_].size()) {
+ std::cerr << "Coface iterator: the argument simplex is not a permutahedral representation\n";
+ is_end_ = true;
+ return;
+ }
+ for (uint i = 0; i < k_ + 1; i++)
+ o_its_.emplace_back(simplex_.partition()[i].size(), (*i_it_)[i] + 1);
+ update_value();
+ }
+
+ // Used for the creating an end iterator
+ Coface_iterator() : is_end_(true) {}
+
+ private:
+ Permutahedral_representation simplex_; // Input simplex
+ uint d_; // Ambient dimension
+ uint l_; // Dimension of the coface
+ uint k_; // Dimension of the input simplex
+ uint t_; // The position of d in simplex_.partition()[k_]
+ Integer_combination_iterator i_it_, i_end_; // indicates in how many parts each simplex_[i] is subdivided
+ std::vector<Ordered_set_partition_iterator> o_its_; // indicates subdivision for each simplex_[i]
+ Ordered_set_partition_iterator o_end_; // one end for all o_its_
+
+ bool is_end_; // is true when the current permutation is the final one
+ value_t value_; // the dereference value
+
+}; // Coface_iterator
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h
new file mode 100644
index 00000000..0f91d41c
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h
@@ -0,0 +1,120 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_
+
+#include <cstdlib> // for std::size_t
+#include <vector>
+
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate permutations.
+ * Based on the optimization of the Heap's algorithm by Sedgewick.
+ */
+class Permutation_iterator
+ : public boost::iterator_facade<Permutation_iterator, std::vector<uint> const, boost::forward_traversal_tag> {
+ using value_t = std::vector<uint>;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ bool equal(Permutation_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void swap_two_indices(std::size_t i, std::size_t j) {
+ uint t = value_[i];
+ value_[i] = value_[j];
+ value_[j] = t;
+ }
+
+ void elementary_increment() {
+ uint j = 0;
+ while (d_[j] == j + 1) {
+ d_[j] = 0;
+ ++j;
+ }
+ if (j == n_ - 1) {
+ is_end_ = true;
+ return;
+ }
+ uint k = j + 1;
+ uint x = (k % 2 ? d_[j] : 0);
+ swap_two_indices(k, x);
+ ++d_[j];
+ }
+
+ void elementary_increment_optim_3() {
+ if (ct_ != 0) {
+ --ct_;
+ swap_two_indices(1 + (ct_ % 2), 0);
+ } else {
+ ct_ = 5;
+ uint j = 2;
+ while (d_[j] == j + 1) {
+ d_[j] = 0;
+ ++j;
+ }
+ if (j == n_ - 1) {
+ is_end_ = true;
+ return;
+ }
+ uint k = j + 1;
+ uint x = (k % 2 ? d_[j] : 0);
+ swap_two_indices(k, x);
+ ++d_[j];
+ }
+ }
+
+ void increment() {
+ if (optim_3_)
+ elementary_increment_optim_3();
+ else
+ elementary_increment();
+ }
+
+ public:
+ Permutation_iterator(const uint& n) : value_(n), is_end_(n == 0), optim_3_(n >= 3), n_(n), d_(n), ct_(5) {
+ for (uint i = 0; i < n; ++i) {
+ value_[i] = i;
+ d_[i] = 0;
+ }
+ if (n > 0) d_[n - 1] = -1;
+ }
+
+ // Used for the creating an end iterator
+ Permutation_iterator() : is_end_(true), n_(0) {}
+
+ void reinitialize() {
+ if (n_ > 0) is_end_ = false;
+ }
+
+ private:
+ value_t value_; // the dereference value
+ bool is_end_; // is true when the current permutation is the final one
+ bool optim_3_; // true if n>=3. for n >= 3, the algorithm is optimized
+
+ uint n_;
+ std::vector<uint> d_; // mix radix digits with radix [2 3 4 ... n-1 (sentinel=-1)]
+ uint ct_; // counter with values in {0,...,5} used in the n>=3 optimization.
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h
new file mode 100644
index 00000000..94ac10c2
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h
@@ -0,0 +1,111 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_
+
+#include <vector>
+#include <limits>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts.
+ *
+ */
+class Set_partition_iterator
+ : public boost::iterator_facade<Set_partition_iterator, std::vector<std::vector<uint>> const,
+ boost::forward_traversal_tag> {
+ using value_t = std::vector<std::vector<uint>>;
+
+ private:
+ friend class boost::iterator_core_access;
+
+ bool equal(Set_partition_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ value_t const& dereference() const { return value_; }
+
+ void update_value() {
+ for (uint i = 0; i < k_; i++) value_[i].clear();
+ for (uint i = 0; i < n_; i++) value_[rgs_[i]].push_back(i);
+ }
+
+ void increment() {
+ if (k_ <= 1) {
+ is_end_ = true;
+ return;
+ }
+ uint i = n_ - 1;
+ while (rgs_[i] + 1 > max_[i] || rgs_[i] + 1 >= k_) i--;
+ if (i == 0) {
+ is_end_ = true;
+ return;
+ }
+ rgs_[i]++;
+ uint mm = max_[i];
+ mm += (rgs_[i] >= mm);
+ max_[i + 1] = mm;
+ while (++i < n_) {
+ rgs_[i] = 0;
+ max_[i + 1] = mm;
+ }
+ uint p = k_;
+ if (mm < p) do {
+ max_[i] = p;
+ --i;
+ --p;
+ rgs_[i] = p;
+ } while (max_[i] < p);
+ update_value();
+ }
+
+ public:
+ Set_partition_iterator(const uint& n, const uint& k)
+ : value_(k), rgs_(n, 0), max_(n + 1), is_end_(n == 0), n_(n), k_(k) {
+ max_[0] = std::numeric_limits<uint>::max();
+ for (uint i = 0; i <= n - k; ++i) value_[0].push_back(i);
+ for (uint i = n - k + 1, j = 1; i < n; ++i, ++j) {
+ rgs_[i] = j;
+ value_[j].push_back(i);
+ }
+ for (uint i = 1; i <= n; i++) max_[i] = rgs_[i - 1] + 1;
+ update_value();
+ }
+
+ // Used for creating an end iterator
+ Set_partition_iterator() : is_end_(true), n_(0), k_(0) {}
+
+ void reinitialize() {
+ if (n_ > 0) is_end_ = false;
+ for (uint i = 0; i <= n_ - k_; ++i) rgs_[i] = 0;
+ for (uint i = n_ - k_ + 1, j = 1; i < n_; ++i, ++j) rgs_[i] = j;
+ for (uint i = 1; i <= n_; i++) max_[i] = rgs_[i - 1] + 1;
+ update_value();
+ }
+
+ private:
+ value_t value_; // the dereference value
+ std::vector<uint> rgs_; // restricted growth string
+ std::vector<uint> max_; // max_[i] = max(rgs_[0],...,rgs[i-1]) + 1
+ bool is_end_; // is true when the current permutation is the final one
+
+ uint n_;
+ uint k_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h
new file mode 100644
index 00000000..905d68d5
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h
@@ -0,0 +1,54 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Simplex_comparator
+ * \brief A comparator class for Permutahedral_representation.
+ * The comparison is in lexicographic order first on
+ * vertices and then on ordered partitions with sorted parts.
+ * The lexicographic order forces that any face is larger than
+ * a coface.
+ *
+ * \tparam Permutahdral_representation_ Needs to be
+ * Permutahedral_representation<Vertex_, Ordered_set_partition_>
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Permutahedral_representation_>
+struct Simplex_comparator {
+ /** \brief Comparison between two permutahedral representations.
+ * Both permutahedral representations need to be valid and
+ * the vertices of both permutahedral representations need to be of the same size.
+ */
+ bool operator()(const Permutahedral_representation_& lhs, const Permutahedral_representation_& rhs) const {
+ if (lhs.vertex() < rhs.vertex()) return true;
+ if (lhs.vertex() > rhs.vertex()) return false;
+
+ if (lhs.partition().size() > rhs.partition().size()) return true;
+ if (lhs.partition().size() < rhs.partition().size()) return false;
+
+ if (lhs.partition() < rhs.partition()) return true;
+ if (lhs.partition() > rhs.partition()) return false;
+
+ return false;
+ }
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h
new file mode 100644
index 00000000..c43effc8
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h
@@ -0,0 +1,73 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_
+
+#include <cstdlib> // for std::size_t
+
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief Auxillary iterator class for sizes of parts in an ordered set partition.
+ */
+template <class T_it>
+class Size_iterator
+ : public boost::iterator_facade<Size_iterator<T_it>, std::size_t const, boost::forward_traversal_tag> {
+ friend class boost::iterator_core_access;
+
+ private:
+ bool equal(Size_iterator const& other) const { return (is_end_ && other.is_end_); }
+
+ std::size_t const& dereference() const { return value_; }
+
+ void increment() {
+ if (++t_it_ == t_end_) {
+ is_end_ = true;
+ return;
+ }
+ value_ = t_it_->size() - 1;
+ }
+
+ public:
+ Size_iterator(const T_it& t_begin, const T_it& t_end) : t_it_(t_begin), t_end_(t_end), is_end_(t_begin == t_end) {
+ if (!is_end_) value_ = t_it_->size() - 1;
+ }
+
+ private:
+ T_it t_it_, t_end_;
+ bool is_end_;
+ std::size_t value_;
+};
+
+template <class T>
+class Size_range {
+ const T& t_;
+
+ public:
+ typedef Size_iterator<typename T::const_iterator> iterator;
+
+ Size_range(const T& t) : t_(t) {}
+
+ std::size_t operator[](std::size_t i) const { return t_[i].size() - 1; }
+
+ iterator begin() const { return iterator(t_.begin(), t_.end()); }
+
+ iterator end() const { return iterator(t_.end(), t_.end()); }
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h
new file mode 100644
index 00000000..47120689
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h
@@ -0,0 +1,66 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_
+#define PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_
+
+#include <cstdlib> // for std::size_t
+#include <algorithm>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief Computes the permutahedral representation of a face of a given simplex
+ * and a range of the vertex indices that compose the face.
+ *
+ * \tparam Permutahedral_representation has to be Permutahedral_representation
+ * \tparam Index_range is a range of unsigned integers taking values in 0,...,k,
+ * where k is the dimension of the simplex simplex.
+ *
+ * @param[in] simplex Input simplex.
+ * @param[in] indices Input range of indices.
+ */
+template <class Permutahedral_representation, class Index_range>
+Permutahedral_representation face_from_indices(const Permutahedral_representation& simplex,
+ const Index_range& indices) {
+ using range_index = typename Index_range::value_type;
+ using Ordered_set_partition = typename Permutahedral_representation::OrderedSetPartition;
+ using Part = typename Ordered_set_partition::value_type;
+ using part_index = typename Part::value_type;
+ Permutahedral_representation value;
+ std::size_t d = simplex.vertex().size();
+ value.vertex() = simplex.vertex();
+ std::size_t k = indices.size() - 1;
+ value.partition().resize(k + 1);
+ std::size_t l = simplex.partition().size() - 1;
+ for (std::size_t h = 1; h < k + 1; h++)
+ for (range_index i = indices[h - 1]; i < indices[h]; i++)
+ for (part_index j : simplex.partition()[i]) value.partition()[h - 1].push_back(j);
+ for (range_index i = indices[k]; i < l + 1; i++)
+ for (part_index j : simplex.partition()[i]) value.partition()[k].push_back(j);
+ for (range_index i = 0; i < indices[0]; i++)
+ for (part_index j : simplex.partition()[i]) {
+ if (j != d)
+ value.vertex()[j]++;
+ else
+ for (std::size_t l = 0; l < d; l++) value.vertex()[l]--;
+ value.partition()[k].push_back(j);
+ }
+ // sort the values in each part (probably not needed)
+ for (auto& part : value.partition()) std::sort(part.begin(), part.end());
+ return value;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/test/CMakeLists.txt b/src/Coxeter_triangulation/test/CMakeLists.txt
new file mode 100644
index 00000000..74ded91e
--- /dev/null
+++ b/src/Coxeter_triangulation/test/CMakeLists.txt
@@ -0,0 +1,30 @@
+project(Coxeter_triangulation_test)
+
+include(GUDHI_boost_test)
+
+if (NOT EIGEN3_VERSION VERSION_LESS 3.1.0)
+ add_executable ( Coxeter_triangulation_permutahedral_representation_test perm_rep_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_permutahedral_representation_test)
+
+ add_executable ( Coxeter_triangulation_freudenthal_triangulation_test freud_triang_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_freudenthal_triangulation_test)
+
+ add_executable ( Coxeter_triangulation_functions_test function_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_functions_test)
+
+ # because of random_orthogonal_matrix inclusion
+ if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+ add_executable ( Coxeter_triangulation_random_orthogonal_matrix_function_test random_orthogonal_matrix_function_test.cpp )
+ target_link_libraries(Coxeter_triangulation_random_orthogonal_matrix_function_test ${CGAL_LIBRARY})
+ gudhi_add_boost_test(Coxeter_triangulation_random_orthogonal_matrix_function_test)
+ endif()
+
+ add_executable ( Coxeter_triangulation_oracle_test oracle_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_oracle_test)
+
+ add_executable ( Coxeter_triangulation_manifold_tracing_test manifold_tracing_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_manifold_tracing_test)
+
+ add_executable ( Coxeter_triangulation_cell_complex_test cell_complex_test.cpp )
+ gudhi_add_boost_test(Coxeter_triangulation_cell_complex_test)
+endif() \ No newline at end of file
diff --git a/src/Coxeter_triangulation/test/cell_complex_test.cpp b/src/Coxeter_triangulation/test/cell_complex_test.cpp
new file mode 100644
index 00000000..4f7f3ec5
--- /dev/null
+++ b/src/Coxeter_triangulation/test/cell_complex_test.cpp
@@ -0,0 +1,59 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "cell_complex"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <gudhi/Debug_utils.h>
+#include <gudhi/IO/output_debug_traces_to_html.h>
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Function_torus_in_R3.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+#include <gudhi/Manifold_tracing.h>
+#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(cell_complex) {
+ double radius = 1.1111;
+ Function_torus_in_R3 fun_torus(radius, 3 * radius);
+ Eigen::VectorXd seed = fun_torus.seed();
+ Function_Sm_in_Rd fun_bound(2.5 * radius, 2, seed);
+
+ auto oracle = make_oracle(fun_torus, fun_bound);
+ double lambda = 0.2;
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+ cox_tr.change_matrix(lambda * cox_tr.matrix());
+
+ using MT = Manifold_tracing<Coxeter_triangulation<> >;
+ using Out_simplex_map = typename MT::Out_simplex_map;
+ std::vector<Eigen::VectorXd> seed_points(1, seed);
+ Out_simplex_map interior_simplex_map, boundary_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map, boundary_simplex_map);
+
+ std::size_t intr_d = oracle.amb_d() - oracle.cod_d();
+ Cell_complex<Out_simplex_map> cell_complex(intr_d);
+ cell_complex.construct_complex(interior_simplex_map, boundary_simplex_map);
+
+ std::size_t interior_sc_map_size0 = cell_complex.interior_simplex_cell_map(0).size();
+ std::size_t interior_sc_map_size1 = cell_complex.interior_simplex_cell_map(1).size();
+ std::size_t interior_sc_map_size2 = cell_complex.interior_simplex_cell_map(2).size();
+ std::size_t boundary_sc_map_size0 = cell_complex.boundary_simplex_cell_map(0).size();
+ std::size_t boundary_sc_map_size1 = cell_complex.boundary_simplex_cell_map(1).size();
+ BOOST_CHECK(interior_simplex_map.size() == interior_sc_map_size0);
+ BOOST_CHECK(boundary_sc_map_size0 - boundary_sc_map_size1 == 0);
+ BOOST_CHECK(interior_sc_map_size0 - interior_sc_map_size1 + interior_sc_map_size2 == 0);
+}
diff --git a/src/Coxeter_triangulation/test/freud_triang_test.cpp b/src/Coxeter_triangulation/test/freud_triang_test.cpp
new file mode 100644
index 00000000..2cf8f00e
--- /dev/null
+++ b/src/Coxeter_triangulation/test/freud_triang_test.cpp
@@ -0,0 +1,114 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "freudenthal_triangulation"
+#include <boost/test/unit_test.hpp>
+
+#include <gudhi/Unitary_tests_utils.h>
+#include <gudhi/Freudenthal_triangulation.h>
+#include <gudhi/Coxeter_triangulation.h>
+
+BOOST_AUTO_TEST_CASE(freudenthal_triangulation) {
+ // Point location check
+ typedef std::vector<double> Point;
+ typedef Gudhi::coxeter_triangulation::Freudenthal_triangulation<> FK_triangulation;
+ typedef typename FK_triangulation::Simplex_handle Simplex_handle;
+ typedef typename FK_triangulation::Vertex_handle Vertex_handle;
+ typedef typename Simplex_handle::OrderedSetPartition Ordered_set_partition;
+ typedef typename Ordered_set_partition::value_type Part;
+
+ FK_triangulation tr(3);
+
+ // Point location check
+ {
+ Point point({3, -1, 0});
+ Simplex_handle s = tr.locate_point(point);
+ BOOST_CHECK(s.vertex() == Vertex_handle({3, -1, 0}));
+ BOOST_CHECK(s.partition() == Ordered_set_partition({Part({0, 1, 2, 3})}));
+ }
+
+ {
+ Point point({3.5, -1.5, 0.5});
+ Simplex_handle s = tr.locate_point(point);
+ BOOST_CHECK(s.vertex() == Vertex_handle({3, -2, 0}));
+ BOOST_CHECK(s.partition() == Ordered_set_partition({Part({0, 1, 2}), Part({3})}));
+ }
+
+ {
+ Point point({3.5, -1.8, 0.5});
+ Simplex_handle s = tr.locate_point(point);
+ BOOST_CHECK(s.vertex() == Vertex_handle({3, -2, 0}));
+ BOOST_CHECK(s.partition() == Ordered_set_partition({Part({0, 2}), Part({1}), Part({3})}));
+ }
+
+ {
+ Point point({3.5, -1.8, 0.3});
+ Simplex_handle s = tr.locate_point(point);
+ BOOST_CHECK(s.vertex() == Vertex_handle({3, -2, 0}));
+ BOOST_CHECK(s.partition() == Ordered_set_partition({Part({0}), Part({2}), Part({1}), Part({3})}));
+ }
+
+ // Dimension check
+ BOOST_CHECK(tr.dimension() == 3);
+ // Matrix check
+ Eigen::MatrixXd default_matrix = Eigen::MatrixXd::Identity(3, 3);
+ BOOST_CHECK(tr.matrix() == default_matrix);
+ // Vector check
+ Eigen::MatrixXd default_offset = Eigen::VectorXd::Zero(3);
+ BOOST_CHECK(tr.offset() == default_offset);
+
+ // Barycenter check
+ Point point({3.5, -1.8, 0.3});
+ Simplex_handle s = tr.locate_point(point);
+ Eigen::Vector3d barycenter_cart = Eigen::Vector3d::Zero();
+ for (auto v : s.vertex_range())
+ for (std::size_t i = 0; i < v.size(); i++) barycenter_cart(i) += v[i];
+ barycenter_cart /= 4.; // simplex is three-dimensional
+ Eigen::Vector3d barycenter = tr.barycenter(s);
+ for (std::size_t i = 0; (long int)i < barycenter.size(); i++)
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(barycenter(i), barycenter_cart(i), 1e-7);
+
+ // Barycenter check for twice the scale
+ s = tr.locate_point(point, 2);
+ barycenter_cart = Eigen::Vector3d::Zero();
+ for (auto v : s.vertex_range())
+ for (std::size_t i = 0; i < v.size(); i++) barycenter_cart(i) += v[i];
+ barycenter_cart /= 3.; // simplex is now a two-dimensional face
+ barycenter_cart /= 2.; // scale
+ barycenter = tr.barycenter(s, 2);
+ for (std::size_t i = 0; (long int)i < barycenter.size(); i++)
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(barycenter(i), barycenter_cart(i), 1e-7);
+
+ // Matrix and offset change check
+ Eigen::MatrixXd new_matrix(3, 3);
+ new_matrix << 1, 0, 0, -1, 1, 0, -1, 0, 1;
+ Eigen::Vector3d new_offset(1.5, 1, 0.5);
+ tr.change_matrix(new_matrix);
+ tr.change_offset(new_offset);
+
+ BOOST_CHECK(tr.matrix() == new_matrix);
+ BOOST_CHECK(tr.offset() == new_offset);
+}
+
+#ifdef GUDHI_DEBUG
+BOOST_AUTO_TEST_CASE(freudenthal_triangulation_exceptions_in_debug_mode) {
+ // Point location check
+ typedef Gudhi::coxeter_triangulation::Freudenthal_triangulation<> FK_triangulation;
+
+ BOOST_CHECK_THROW (FK_triangulation tr(3, Eigen::MatrixXd::Identity(3, 3), Eigen::VectorXd::Zero(4)),
+ std::invalid_argument);
+
+ FK_triangulation tr(3);
+ // Point of dimension 4
+ std::vector<double> point({3.5, -1.8, 0.3, 4.1});
+ BOOST_CHECK_THROW (tr.locate_point(point), std::invalid_argument);
+}
+#endif
diff --git a/src/Coxeter_triangulation/test/function_test.cpp b/src/Coxeter_triangulation/test/function_test.cpp
new file mode 100644
index 00000000..43dbcb75
--- /dev/null
+++ b/src/Coxeter_triangulation/test/function_test.cpp
@@ -0,0 +1,158 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+// workaround for the annoying boost message in boost 1.69
+#define BOOST_PENDING_INTEGER_LOG2_HPP
+#include <boost/integer/integer_log2.hpp>
+// end workaround
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "function"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Function_affine_plane_in_Rd.h>
+#include <gudhi/Functions/Constant_function.h>
+#include <gudhi/Functions/Function_chair_in_R3.h>
+#include <gudhi/Functions/Function_torus_in_R3.h>
+#include <gudhi/Functions/Function_whitney_umbrella_in_R3.h>
+#include <gudhi/Functions/Function_lemniscate_revolution_in_R3.h>
+#include <gudhi/Functions/Function_iron_in_R3.h>
+#include <gudhi/Functions/Function_moment_curve_in_Rd.h>
+#include <gudhi/Functions/Embed_in_Rd.h>
+#include <gudhi/Functions/Translate.h>
+#include <gudhi/Functions/Linear_transformation.h>
+#include <gudhi/Functions/Negation.h>
+#include <gudhi/Functions/Cartesian_product.h>
+#include <gudhi/Functions/PL_approximation.h>
+
+#include <gudhi/Coxeter_triangulation.h>
+
+#include <string>
+
+#include <random>
+#include <cstdlib>
+
+using namespace Gudhi::coxeter_triangulation;
+
+template <class Function>
+void test_function(const Function& fun) {
+ Eigen::VectorXd seed = fun.seed();
+ Eigen::VectorXd res_seed = fun(fun.seed());
+ BOOST_CHECK(seed.size() == (long int)fun.amb_d());
+ BOOST_CHECK(res_seed.size() == (long int)fun.cod_d());
+ for (std::size_t i = 0; i < fun.cod_d(); i++) GUDHI_TEST_FLOAT_EQUALITY_CHECK(res_seed(i), 0., 1e-10);
+}
+
+BOOST_AUTO_TEST_CASE(function) {
+ {
+ // the sphere testing part
+ std::size_t m = 3, d = 5;
+ Eigen::VectorXd center(d);
+ center << 2, 1.5, -0.5, 4.5, -1;
+ double radius = 5;
+ typedef Function_Sm_in_Rd Function_sphere;
+ Function_sphere fun_sphere(radius, m, d, center);
+ test_function(fun_sphere);
+ }
+ {
+ // the affine plane testing part
+ std::size_t m = 0, d = 5;
+ Eigen::MatrixXd normal_matrix = Eigen::MatrixXd::Zero(d, d - m);
+ for (std::size_t i = 0; i < d - m; ++i) normal_matrix(i, i) = 1;
+ typedef Function_affine_plane_in_Rd Function_plane;
+ Function_plane fun_plane(normal_matrix);
+ test_function(fun_plane);
+ }
+ {
+ // the constant function testing part
+ std::size_t k = 2, d = 5;
+ auto x = Eigen::VectorXd::Constant(k, 1);
+ Constant_function fun_const(d, k, x);
+ Eigen::VectorXd res_zero = fun_const(Eigen::VectorXd::Zero(d));
+ for (std::size_t i = 0; i < k; ++i) GUDHI_TEST_FLOAT_EQUALITY_CHECK(res_zero(i), x(i), 1e-10);
+ }
+ {
+ // the chair function
+ Function_chair_in_R3 fun_chair;
+ test_function(fun_chair);
+ }
+ {
+ // the torus function
+ Function_torus_in_R3 fun_torus;
+ test_function(fun_torus);
+ }
+ {
+ // the whitney umbrella function
+ Function_whitney_umbrella_in_R3 fun_umbrella;
+ test_function(fun_umbrella);
+ }
+ {
+ // the lemniscate revolution function
+ Function_lemniscate_revolution_in_R3 fun_lemniscate;
+ test_function(fun_lemniscate);
+ }
+ {
+ // the iron function
+ Function_iron_in_R3 fun_iron;
+ test_function(fun_iron);
+ }
+ {
+ Function_moment_curve_in_Rd fun_moment_curve(3, 5);
+ test_function(fun_moment_curve);
+ }
+ {
+ // function embedding
+ Function_iron_in_R3 fun_iron;
+ auto fun_embed = make_embedding(fun_iron, 5);
+ test_function(fun_iron);
+
+ // function translation
+ Eigen::VectorXd off = Eigen::VectorXd::Random(5);
+ auto fun_trans = translate(fun_embed, off);
+ test_function(fun_trans);
+
+ // function linear transformation
+ Eigen::MatrixXd matrix = Eigen::MatrixXd::Random(5, 5);
+ BOOST_CHECK(matrix.determinant() != 0.);
+ auto fun_lin = make_linear_transformation(fun_trans, matrix);
+ test_function(fun_lin);
+
+ // function negative
+ auto fun_neg = negation(fun_lin);
+ test_function(fun_neg);
+
+ // function product
+ typedef Function_Sm_in_Rd Function_sphere;
+ Function_sphere fun_sphere(1, 1);
+ auto fun_prod = make_product_function(fun_sphere, fun_sphere, fun_sphere);
+ test_function(fun_prod);
+
+ // function PL approximation
+ Coxeter_triangulation<> cox_tr(6);
+ typedef Coxeter_triangulation<>::Vertex_handle Vertex_handle;
+ auto fun_pl = make_pl_approximation(fun_prod, cox_tr);
+ Vertex_handle v0 = Vertex_handle(cox_tr.dimension(), 0);
+ Eigen::VectorXd x0 = cox_tr.cartesian_coordinates(v0);
+ Eigen::VectorXd value0 = fun_prod(x0);
+ Eigen::VectorXd pl_value0 = fun_pl(x0);
+ for (std::size_t i = 0; i < fun_pl.cod_d(); i++) GUDHI_TEST_FLOAT_EQUALITY_CHECK(value0(i), pl_value0(i), 1e-10);
+ Vertex_handle v1 = v0;
+ v1[0] += 1;
+ Eigen::VectorXd x1 = cox_tr.cartesian_coordinates(v1);
+ Eigen::VectorXd value1 = fun_prod(x1);
+ Eigen::VectorXd pl_value1 = fun_pl(x1);
+ for (std::size_t i = 0; i < fun_pl.cod_d(); i++) GUDHI_TEST_FLOAT_EQUALITY_CHECK(value1(i), pl_value1(i), 1e-10);
+ Eigen::VectorXd pl_value_mid = fun_pl(0.5 * x0 + 0.5 * x1);
+ for (std::size_t i = 0; i < fun_pl.cod_d(); i++)
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(0.5 * value0(i) + 0.5 * value1(i), pl_value_mid(i), 1e-10);
+ }
+}
diff --git a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
new file mode 100644
index 00000000..63497f5a
--- /dev/null
+++ b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
@@ -0,0 +1,62 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "manifold_tracing"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+#include <gudhi/Manifold_tracing.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(manifold_tracing) {
+ // manifold without boundary
+ Function_Sm_in_Rd fun_sph(5.1111, 2);
+ auto oracle = make_oracle(fun_sph);
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ // cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+
+ using MT = Manifold_tracing<Coxeter_triangulation<> >;
+ Eigen::VectorXd seed = fun_sph.seed();
+ std::vector<Eigen::VectorXd> seed_points(1, seed);
+ typename MT::Out_simplex_map out_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle, out_simplex_map);
+
+ for (auto si_pair : out_simplex_map) {
+ BOOST_CHECK(si_pair.first.dimension() == oracle.function().cod_d());
+ BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d());
+ }
+ std::clog << "out_simplex_map.size() = " << out_simplex_map.size() << "\n";
+ BOOST_CHECK(out_simplex_map.size() == 1118);
+
+ // manifold with boundary
+ Function_Sm_in_Rd fun_boundary(3.0, 2, fun_sph.seed());
+ auto oracle_with_boundary = make_oracle(fun_sph, fun_boundary);
+ typename MT::Out_simplex_map interior_simplex_map, boundary_simplex_map;
+ manifold_tracing_algorithm(seed_points, cox_tr, oracle_with_boundary, interior_simplex_map, boundary_simplex_map);
+ for (auto si_pair : interior_simplex_map) {
+ BOOST_CHECK(si_pair.first.dimension() == oracle.function().cod_d());
+ BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d());
+ }
+ std::clog << "interior_simplex_map.size() = " << interior_simplex_map.size() << "\n";
+ BOOST_CHECK(interior_simplex_map.size() == 96);
+ for (auto si_pair : boundary_simplex_map) {
+ BOOST_CHECK(si_pair.first.dimension() == oracle.function().cod_d() + 1);
+ BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d());
+ }
+ std::clog << "boundary_simplex_map.size() = " << boundary_simplex_map.size() << "\n";
+ BOOST_CHECK(boundary_simplex_map.size() == 54);
+}
diff --git a/src/Coxeter_triangulation/test/oracle_test.cpp b/src/Coxeter_triangulation/test/oracle_test.cpp
new file mode 100644
index 00000000..ed2042f5
--- /dev/null
+++ b/src/Coxeter_triangulation/test/oracle_test.cpp
@@ -0,0 +1,56 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "oracle"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <string>
+
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Cartesian_product.h>
+
+#include <gudhi/Coxeter_triangulation.h>
+
+#include <random>
+#include <cstdlib>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(oracle) {
+ Function_Sm_in_Rd fun_sph(5.1111, 2);
+ auto oracle = make_oracle(fun_sph);
+ Coxeter_triangulation<> cox_tr(oracle.amb_d());
+ // cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+
+ Eigen::VectorXd seed = fun_sph.seed();
+ auto s = cox_tr.locate_point(seed);
+
+ std::size_t num_intersected_edges = 0;
+ for (auto f : s.face_range(oracle.cod_d())) {
+ auto qr = oracle.intersects(f, cox_tr);
+ if (qr.success) num_intersected_edges++;
+ auto vertex_it = f.vertex_range().begin();
+ Eigen::Vector3d p1 = cox_tr.cartesian_coordinates(*vertex_it++);
+ Eigen::Vector3d p2 = cox_tr.cartesian_coordinates(*vertex_it++);
+ BOOST_CHECK(vertex_it == f.vertex_range().end());
+ Eigen::MatrixXd m(3, 3);
+ if (qr.success) {
+ m.col(0) = qr.intersection;
+ m.col(1) = p1;
+ m.col(2) = p2;
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(m.determinant(), 0.0, 1e-10);
+ }
+ }
+ BOOST_CHECK(num_intersected_edges == 3 || num_intersected_edges == 4);
+}
diff --git a/src/Coxeter_triangulation/test/perm_rep_test.cpp b/src/Coxeter_triangulation/test/perm_rep_test.cpp
new file mode 100644
index 00000000..a668fc66
--- /dev/null
+++ b/src/Coxeter_triangulation/test/perm_rep_test.cpp
@@ -0,0 +1,61 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "permutahedral_representation"
+#include <boost/test/unit_test.hpp>
+
+#include <gudhi/Permutahedral_representation.h>
+
+BOOST_AUTO_TEST_CASE(permutahedral_representation) {
+ typedef std::vector<int> Vertex;
+ typedef std::vector<std::size_t> Part;
+ typedef std::vector<Part> Partition;
+ typedef Gudhi::coxeter_triangulation::Permutahedral_representation<Vertex, Partition> Simplex_handle;
+ Vertex v0(10, 0);
+ Partition omega = {Part({5}), Part({2}), Part({3, 7}), Part({4, 9}), Part({0, 6, 8}), Part({1, 10})};
+ Simplex_handle s(v0, omega);
+
+ // Dimension check
+ BOOST_CHECK(s.dimension() == 5);
+
+ // Vertex number check
+ std::vector<Vertex> vertices;
+ for (auto& v : s.vertex_range()) vertices.push_back(v);
+ BOOST_CHECK(vertices.size() == 6);
+
+ // Facet number check
+ std::vector<Simplex_handle> facets;
+ for (auto& f : s.facet_range()) facets.push_back(f);
+ BOOST_CHECK(facets.size() == 6);
+
+ // Face of dim 3 number check
+ std::vector<Simplex_handle> faces3;
+ for (auto& f : s.face_range(3)) faces3.push_back(f);
+ BOOST_CHECK(faces3.size() == 15);
+
+ // Cofacet number check
+ std::vector<Simplex_handle> cofacets;
+ for (auto& f : s.cofacet_range()) cofacets.push_back(f);
+ BOOST_CHECK(cofacets.size() == 12);
+
+ // Is face check
+ Vertex v1(10, 0);
+ Partition omega1 = {Part({5}), Part({0, 1, 2, 3, 4, 6, 7, 8, 9, 10})};
+ Simplex_handle s1(v1, omega1);
+ Vertex v2(10, 0);
+ v2[1] = -1;
+ Partition omega2 = {Part({1}), Part({5}), Part({2}), Part({3, 7}), Part({4, 9}), Part({0, 6, 8}), Part({10})};
+ Simplex_handle s2(v2, omega2);
+ BOOST_CHECK(s.is_face_of(s));
+ BOOST_CHECK(s1.is_face_of(s));
+ BOOST_CHECK(!s2.is_face_of(s));
+ BOOST_CHECK(s.is_face_of(s2));
+}
diff --git a/src/Coxeter_triangulation/test/random_orthogonal_matrix_function_test.cpp b/src/Coxeter_triangulation/test/random_orthogonal_matrix_function_test.cpp
new file mode 100644
index 00000000..84178741
--- /dev/null
+++ b/src/Coxeter_triangulation/test/random_orthogonal_matrix_function_test.cpp
@@ -0,0 +1,36 @@
+/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ * Author(s): Siargey Kachanovich
+ *
+ * Copyright (C) 2019 Inria
+ *
+ * Modification(s):
+ * - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "random_orthogonal_matrix_function"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <gudhi/Functions/random_orthogonal_matrix.h>
+
+#include <string>
+
+#include <random>
+#include <cstdlib>
+
+using namespace Gudhi::coxeter_triangulation;
+
+// this test is separated as it requires CGAL
+BOOST_AUTO_TEST_CASE(random_orthogonal_matrix_function) {
+ // random orthogonal matrix
+ Eigen::MatrixXd matrix = random_orthogonal_matrix(5);
+ Eigen::MatrixXd id_matrix = matrix.transpose() * matrix;
+ for (std::size_t i = 0; i < 5; ++i)
+ for (std::size_t j = 0; j < 5; ++j)
+ if (i == j)
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(id_matrix(i, j), 1.0, 1e-10);
+ else
+ GUDHI_TEST_FLOAT_EQUALITY_CHECK(id_matrix(i, j), 0.0, 1e-10);
+}
diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h
index a2de4335..879fb96a 100644
--- a/src/common/doc/examples.h
+++ b/src/common/doc/examples.h
@@ -1,130 +1,134 @@
// List of GUDHI examples and utils - Doxygen needs at least a file tag to analyse comments
// Generated from scripts/cpp_examples_for_doxygen.py
/*! @file Examples
- * \section Skeleton_blocker_example_section Skeleton_blocker
- * @example Skeleton_blocker_iteration.cpp
- * @example Skeleton_blocker_from_simplices.cpp
- * @example Skeleton_blocker_link.cpp
- * \section Alpha_complex_example_section Alpha_complex
- * @example alpha_complex_persistence.cpp
- * @example alpha_complex_3d_persistence.cpp
- * @example Alpha_complex_from_points.cpp
- * @example Alpha_complex_3d_from_points.cpp
- * @example Weighted_alpha_complex_3d_from_points.cpp
- * @example Fast_alpha_complex_from_off.cpp
- * @example Alpha_complex_from_off.cpp
- * @example Weighted_alpha_complex_from_points.cpp
+ * \section Witness_complex_example_section Witness_complex
+ * @example strong_witness_persistence.cpp
+ * @example weak_witness_persistence.cpp
+ * @example example_witness_complex_off.cpp
+ * @example example_strong_witness_complex_off.cpp
+ * @example example_nearest_landmark_table.cpp
+ * @example example_witness_complex_sphere.cpp
+ * \section Contraction_example_section Contraction
+ * @example Rips_contraction.cpp
+ * @example Garland_heckbert.cpp
* \section Simplex_tree_example_section Simplex_tree
- * @example cech_complex_cgal_mini_sphere_3d.cpp
* @example mini_simplex_tree.cpp
- * @example example_alpha_shapes_3_simplex_tree_from_off_file.cpp
+ * @example cech_complex_cgal_mini_sphere_3d.cpp
* @example graph_expansion_with_blocker.cpp
* @example simple_simplex_tree.cpp
* @example simplex_tree_from_cliques_of_graph.cpp
- * \section Collapse_example_section Collapse
- * @example point_cloud_edge_collapse_rips_persistence.cpp
- * @example distance_matrix_edge_collapse_rips_persistence.cpp
- * @example edge_collapse_basic_example.cpp
- * @example edge_collapse_conserve_persistence.cpp
+ * @example example_alpha_shapes_3_simplex_tree_from_off_file.cpp
* \section Persistent_cohomology_example_section Persistent_cohomology
- * @example rips_persistence_via_boundary_matrix.cpp
- * @example rips_persistence_step_by_step.cpp
- * @example plain_homology.cpp
* @example custom_persistence_sort.cpp
- * @example persistence_from_simple_simplex_tree.cpp
+ * @example rips_persistence_step_by_step.cpp
* @example persistence_from_file.cpp
+ * @example rips_persistence_via_boundary_matrix.cpp
+ * @example plain_homology.cpp
* @example rips_multifield_persistence.cpp
+ * @example persistence_from_simple_simplex_tree.cpp
+ * \section Subsampling_example_section Subsampling
+ * @example example_sparsify_point_set.cpp
+ * @example example_choose_n_farthest_points.cpp
+ * @example example_custom_distance.cpp
+ * @example example_pick_n_random_points.cpp
+ * \section Toplex_map_example_section Toplex_map
+ * @example simple_toplex_map.cpp
+ * \section Collapse_example_section Collapse
+ * @example distance_matrix_edge_collapse_rips_persistence.cpp
+ * @example point_cloud_edge_collapse_rips_persistence.cpp
+ * @example edge_collapse_conserve_persistence.cpp
+ * @example edge_collapse_basic_example.cpp
+ * \section Cech_complex_example_section Cech_complex
+ * @example cech_persistence.cpp
+ * @example cech_complex_step_by_step.cpp
+ * @example cech_complex_example_from_points.cpp
+ * \section Bitmap_cubical_complex_example_section Bitmap_cubical_complex
+ * @example periodic_cubical_complex_persistence.cpp
+ * @example cubical_complex_persistence.cpp
+ * @example Random_bitmap_cubical_complex.cpp
+ * \section Coxeter_triangulation_example_section Coxeter_triangulation
+ * @example cell_complex_from_basic_circle_manifold.cpp
+ * @example manifold_tracing_flat_torus_with_boundary.cpp
+ * @example manifold_tracing_custom_function.cpp
* \section Nerve_GIC_example_section Nerve_GIC
* @example VoronoiGIC.cpp
* @example Nerve.cpp
- * @example FuncGIC.cpp
* @example CoordGIC.cpp
- * \section Rips_complex_example_section Rips_complex
- * @example rips_persistence.cpp
- * @example sparse_rips_persistence.cpp
- * @example rips_correlation_matrix_persistence.cpp
- * @example rips_distance_matrix_persistence.cpp
- * @example example_one_skeleton_rips_from_correlation_matrix.cpp
- * @example example_rips_complex_from_csv_distance_matrix_file.cpp
- * @example example_sparse_rips.cpp
- * @example example_one_skeleton_rips_from_distance_matrix.cpp
- * @example example_rips_complex_from_off_file.cpp
- * @example example_one_skeleton_rips_from_points.cpp
+ * @example FuncGIC.cpp
+ * \section Tangential_complex_example_section Tangential_complex
+ * @example example_basic.cpp
+ * @example example_with_perturb.cpp
* \section Persistence_representations_example_section Persistence_representations
+ * @example persistence_vectors/create_persistence_vectors.cpp
+ * @example persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
+ * @example persistence_vectors/plot_persistence_vectors.cpp
+ * @example persistence_vectors/average_persistence_vectors.cpp
+ * @example persistence_vectors/compute_distance_of_persistence_vectors.cpp
* @example persistence_landscapes_on_grid/average_landscapes_on_grid.cpp
* @example persistence_landscapes_on_grid/create_landscapes_on_grid.cpp
* @example persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp
- * @example persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp
* @example persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp
+ * @example persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp
* @example persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp
- * @example persistence_intervals/plot_persistence_intervals.cpp
- * @example persistence_intervals/plot_histogram_of_intervals_lengths.cpp
- * @example persistence_intervals/compute_bottleneck_distance.cpp
* @example persistence_intervals/plot_persistence_Betti_numbers.cpp
+ * @example persistence_intervals/compute_bottleneck_distance.cpp
* @example persistence_intervals/compute_number_of_dominant_intervals.cpp
- * @example persistence_heat_maps/average_persistence_heat_maps.cpp
- * @example persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp
- * @example persistence_heat_maps/create_pssk.cpp
- * @example persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp
+ * @example persistence_intervals/plot_histogram_of_intervals_lengths.cpp
+ * @example persistence_intervals/plot_persistence_intervals.cpp
* @example persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp
+ * @example persistence_heat_maps/create_pssk.cpp
* @example persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp
- * @example persistence_heat_maps/plot_persistence_heat_map.cpp
+ * @example persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp
+ * @example persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp
* @example persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp
+ * @example persistence_heat_maps/plot_persistence_heat_map.cpp
* @example persistence_heat_maps/create_persistence_heat_maps.cpp
- * @example persistence_vectors/average_persistence_vectors.cpp
- * @example persistence_vectors/plot_persistence_vectors.cpp
- * @example persistence_vectors/create_persistence_vectors.cpp
- * @example persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp
- * @example persistence_vectors/compute_distance_of_persistence_vectors.cpp
- * @example persistence_landscapes/create_landscapes.cpp
- * @example persistence_landscapes/average_landscapes.cpp
+ * @example persistence_heat_maps/average_persistence_heat_maps.cpp
* @example persistence_landscapes/compute_distance_of_landscapes.cpp
- * @example persistence_landscapes/plot_landscapes.cpp
* @example persistence_landscapes/compute_scalar_product_of_landscapes.cpp
- * @example persistence_heat_maps.cpp
- * @example sliced_wasserstein.cpp
+ * @example persistence_landscapes/average_landscapes.cpp
+ * @example persistence_landscapes/plot_landscapes.cpp
+ * @example persistence_landscapes/create_landscapes.cpp
+ * @example persistence_landscape_on_grid.cpp
* @example persistence_intervals.cpp
- * @example persistence_vectors.cpp
* @example persistence_landscape.cpp
- * @example persistence_landscape_on_grid.cpp
- * \section common_example_section common
- * @example off_file_from_shape_generator.cpp
- * @example example_CGAL_points_off_reader.cpp
- * @example example_vector_double_points_off_reader.cpp
- * @example example_CGAL_3D_points_off_reader.cpp
- * \section Subsampling_example_section Subsampling
- * @example example_pick_n_random_points.cpp
- * @example example_choose_n_farthest_points.cpp
- * @example example_sparsify_point_set.cpp
- * @example example_custom_distance.cpp
- * \section Contraction_example_section Contraction
- * @example Rips_contraction.cpp
- * @example Garland_heckbert.cpp
+ * @example persistence_vectors.cpp
+ * @example sliced_wasserstein.cpp
+ * @example persistence_heat_maps.cpp
+ * \section Spatial_searching_example_section Spatial_searching
+ * @example example_spatial_searching.cpp
* \section Bottleneck_distance_example_section Bottleneck_distance
* @example bottleneck_distance.cpp
* @example bottleneck_basic_example.cpp
* @example alpha_rips_persistence_bottleneck_distance.cpp
- * \section Tangential_complex_example_section Tangential_complex
- * @example example_basic.cpp
- * @example example_with_perturb.cpp
- * \section Cech_complex_example_section Cech_complex
- * @example cech_persistence.cpp
- * @example cech_complex_example_from_points.cpp
- * @example cech_complex_step_by_step.cpp
- * \section Spatial_searching_example_section Spatial_searching
- * @example example_spatial_searching.cpp
- * \section Bitmap_cubical_complex_example_section Bitmap_cubical_complex
- * @example periodic_cubical_complex_persistence.cpp
- * @example cubical_complex_persistence.cpp
- * @example Random_bitmap_cubical_complex.cpp
- * \section Toplex_map_example_section Toplex_map
- * @example simple_toplex_map.cpp
- * \section Witness_complex_example_section Witness_complex
- * @example weak_witness_persistence.cpp
- * @example strong_witness_persistence.cpp
- * @example example_witness_complex_sphere.cpp
- * @example example_strong_witness_complex_off.cpp
- * @example example_nearest_landmark_table.cpp
- * @example example_witness_complex_off.cpp
+ * \section common_example_section common
+ * @example off_file_from_shape_generator.cpp
+ * @example example_vector_double_points_off_reader.cpp
+ * @example example_CGAL_points_off_reader.cpp
+ * @example example_CGAL_3D_points_off_reader.cpp
+ * \section Alpha_complex_example_section Alpha_complex
+ * @example alpha_complex_3d_persistence.cpp
+ * @example alpha_complex_persistence.cpp
+ * @example Fast_alpha_complex_from_off.cpp
+ * @example Alpha_complex_3d_from_points.cpp
+ * @example Alpha_complex_from_off.cpp
+ * @example Weighted_alpha_complex_3d_from_points.cpp
+ * @example Weighted_alpha_complex_from_points.cpp
+ * @example Alpha_complex_from_points.cpp
+ * \section Skeleton_blocker_example_section Skeleton_blocker
+ * @example Skeleton_blocker_from_simplices.cpp
+ * @example Skeleton_blocker_link.cpp
+ * @example Skeleton_blocker_iteration.cpp
+ * \section Rips_complex_example_section Rips_complex
+ * @example rips_persistence.cpp
+ * @example rips_correlation_matrix_persistence.cpp
+ * @example sparse_rips_persistence.cpp
+ * @example rips_distance_matrix_persistence.cpp
+ * @example example_sparse_rips.cpp
+ * @example example_rips_complex_from_csv_distance_matrix_file.cpp
+ * @example example_one_skeleton_rips_from_correlation_matrix.cpp
+ * @example example_one_skeleton_rips_from_distance_matrix.cpp
+ * @example example_one_skeleton_rips_from_points.cpp
+ * @example example_rips_complex_from_off_file.cpp
*/
diff --git a/src/common/doc/installation.h b/src/common/doc/installation.h
index 609f2f2f..ef668dfb 100644
--- a/src/common/doc/installation.h
+++ b/src/common/doc/installation.h
@@ -126,6 +126,8 @@ make \endverbatim
* Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp</a>
* \li <a href="alpha_complex_3d_persistence_8cpp-example.html">
* Alpha_complex/alpha_complex_3d_persistence.cpp</a>
+ * \li <a href="_coxeter_triangulation_2manifold_tracing_flat_torus_with_boundary_8cpp-example.html">
+ * Coxeter_triangulation/manifold_tracing_flat_torus_with_boundary.cpp</a>
*
* \subsection eigen Eigen
* Some GUDHI modules (cf. \ref main_page "modules list"), and few examples require
@@ -170,6 +172,12 @@ make \endverbatim
* Witness_complex/example_witness_complex_off.cpp</a>
* \li <a href="example_witness_complex_sphere_8cpp-example.html">
* Witness_complex/example_witness_complex_sphere.cpp</a>
+ * \li <a href="_coxeter_triangulation_2cell_complex_from_basic_circle_manifold_8cpp-example.html">
+ * Coxeter_triangulation/cell_complex_from_basic_circle_manifold.cpp</a>
+ * \li <a href="_coxeter_triangulation_2manifold_tracing_custom_function_8cpp-example.html">
+ * Coxeter_triangulation/manifold_tracing_custom_function.cpp</a>
+ * \li <a href="_coxeter_triangulation_2manifold_tracing_flat_torus_with_boundary_8cpp-example.html">
+ * Coxeter_triangulation/manifold_tracing_flat_torus_with_boundary.cpp</a>
*
* \subsection tbb Threading Building Blocks
* <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel&reg; TBB</a> lets you easily write parallel
diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md
index e19af537..17354179 100644
--- a/src/common/doc/main_page.md
+++ b/src/common/doc/main_page.md
@@ -135,7 +135,7 @@
</tr>
</table>
-## Filtrations and reconstructions {#FiltrationsReconstructions}
+## Filtrations
### Alpha complex
<table>
@@ -298,6 +298,32 @@
</tr>
</table>
+## Manifold reconstructions
+### Coxeter triangulation
+
+<table>
+ <tr>
+ <td width="35%" rowspan=2>
+ \image html "manifold_tracing_on_custom_function_example.png"
+ </td>
+ <td width="50%">
+ Coxeter triangulation module is designed to provide tools for constructing a piecewise-linear approximation of an
+ \f$m\f$-dimensional smooth manifold embedded in \f$ \mathbb{R}^d \f$ using an ambient triangulation.
+ </td>
+ <td width="15%">
+ <b>Author:</b> Siargey Kachanovich<br>
+ <b>Introduced in:</b> GUDHI 3.4.0<br>
+ <b>Copyright:</b> MIT [(LGPL v3)](../../licensing/)<br>
+ <b>Requires:</b> \ref eigen &ge; 3.1.0
+ </td>
+ </tr>
+ <tr>
+ <td colspan=2 height="25">
+ <b>User manual:</b> \ref coxeter_triangulation
+ </td>
+ </tr>
+</table>
+
### Tangential complex
<table>