diff options
72 files changed, 6631 insertions, 137 deletions
diff --git a/.github/next_release.md b/.github/next_release.md index 26143b0e..170c7494 100644 --- a/.github/next_release.md +++ b/.github/next_release.md @@ -6,8 +6,8 @@ We are now using GitHub to develop the GUDHI library, do not hesitate to [fork t Below is a list of changes made since GUDHI 3.4.0: -- [Module](link) - - ... +- [Coxeter triangulation](https://gudhi.inria.fr/doc/latest/group__coxeter__triangulation.html) + - constructs a piecewise-linear approximation of an m-dimensional smooth manifold embedded in R^d using an ambient triangulation. - [Module](link) - ... diff --git a/CMakeLists.txt b/CMakeLists.txt index 4257a025..d0cf6a25 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,7 @@ add_gudhi_module(Bitmap_cubical_complex) add_gudhi_module(Bottleneck_distance) add_gudhi_module(Collapse) add_gudhi_module(Contraction) +add_gudhi_module(Coxeter_triangulation) add_gudhi_module(Cech_complex) add_gudhi_module(Hasse_complex) add_gudhi_module(Persistence_representations) diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index 16fa29d0..b5afff52 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -15,7 +15,6 @@ title = {{Statistical analysis and parameter selection for Mapper}}, volume = {19}, year = {2018}, url = {http://jmlr.org/papers/v19/17-291.html}, -doi = {10.5555/3291125.3291137} } @inproceedings{Dey13, @@ -1324,3 +1323,18 @@ year = "2011" doi = {10.4230/LIPIcs.SoCG.2020.19}, annote = {Keywords: Computational Topology, Topological Data Analysis, Edge Collapse, Simple Collapse, Persistent homology} } + +@phdthesis{KachanovichThesis, + TITLE = {{Meshing submanifolds using Coxeter triangulations}}, + AUTHOR = {Kachanovich, Siargey}, + URL = {https://hal.inria.fr/tel-02419148}, + NUMBER = {2019AZUR4072}, + SCHOOL = {{COMUE Universit{\'e} C{\^o}te d'Azur (2015 - 2019)}}, + YEAR = {2019}, + MONTH = Oct, + KEYWORDS = {Mesh generation ; Coxeter triangulations ; Simplex quality ; Triangulations of the Euclidean space ; Freudenthal-Kuhn triangulations ; G{\'e}n{\'e}ration de maillages ; Triangulations de Coxeter ; Qualit{\'e} des simplexes ; Triangulations de l'espace euclidien ; Triangulations de Freudenthal-Kuhn}, + TYPE = {Theses}, + PDF = {https://hal.inria.fr/tel-02419148v2/file/2019AZUR4072.pdf}, + HAL_ID = {tel-02419148}, + HAL_VERSION = {v2}, +}
\ No newline at end of file diff --git a/scripts/cpp_examples_for_doxygen.py b/scripts/cpp_examples_for_doxygen.py new file mode 100644 index 00000000..5c091c4f --- /dev/null +++ b/scripts/cpp_examples_for_doxygen.py @@ -0,0 +1,16 @@ +import os +import glob + +for gd_mod in glob.glob("src/*/"): + mod_files = [] + for paths in [gd_mod + 'utilities', gd_mod + 'example']: + if os.path.isdir(paths): + for root, dirs, files in os.walk(paths): + for file in files: + if file.endswith(".cpp"): + mod_files.append(str(os.path.join(root, file)).split(paths)[1][1:]) + if len(mod_files) > 0: + mod = str(gd_mod).split('/')[1] + print(' * \section ' + mod + '_example_section ' + mod) + for file in mod_files: + print(' * @example ' + file) diff --git a/src/Bottleneck_distance/utilities/bottleneckdistance.md b/src/Bottleneck_distance/utilities/bottleneckdistance.md index a81426cf..2f5dedc9 100644 --- a/src/Bottleneck_distance/utilities/bottleneckdistance.md +++ b/src/Bottleneck_distance/utilities/bottleneckdistance.md @@ -10,14 +10,14 @@ Leave the lines above as it is required by the web site generator 'Jekyll' {:/comment}
-## bottleneck_read_file_example ##
+## bottleneck_distance ##
This program computes the Bottleneck distance between two persistence diagram files.
**Usage**
```
- bottleneck_read_file_example <file_1.pers> <file_2.pers> [<tolerance>]
+ bottleneck_distance <file_1.pers> <file_2.pers> [<tolerance>]
```
where
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 Binary files differnew file mode 100644 index 00000000..8bb8ba9a --- /dev/null +++ b/src/Coxeter_triangulation/doc/custom_function.png diff --git a/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png Binary files differnew file mode 100644 index 00000000..338b39fe --- /dev/null +++ b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png 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 Binary files differnew file mode 100644 index 00000000..04912729 --- /dev/null +++ b/src/Coxeter_triangulation/doc/manifold_tracing_on_custom_function_example.png diff --git a/src/Coxeter_triangulation/doc/two_triangulations.png b/src/Coxeter_triangulation/doc/two_triangulations.png Binary files differnew file mode 100644 index 00000000..055d93e7 --- /dev/null +++ b/src/Coxeter_triangulation/doc/two_triangulations.png 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 b557727b..879fb96a 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -1,98 +1,134 @@ -// List of GUDHI examples - Doxygen needs at least a file tag to analyse comments -// In user_version, `find . -name "*.cpp"` in example and utilities folders +// 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 Alpha_complex_examples Alpha complex - * @example Alpha_complex_from_off.cpp - * @example Alpha_complex_from_points.cpp - * \section bottleneck_examples bottleneck - * @example bottleneck_basic_example.cpp - * @example alpha_rips_persistence_bottleneck_distance.cpp - * @example example_nearest_landmark_table.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_witness_complex_sphere.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 mini_simplex_tree.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 * @example example_alpha_shapes_3_simplex_tree_from_off_file.cpp - * @example cech_complex_cgal_mini_sphere_3d.cpp - * @example plain_homology.cpp - * @example persistence_from_file.cpp + * \section Persistent_cohomology_example_section Persistent_cohomology + * @example custom_persistence_sort.cpp * @example rips_persistence_step_by_step.cpp + * @example persistence_from_file.cpp * @example rips_persistence_via_boundary_matrix.cpp - * @example custom_persistence_sort.cpp - * @example persistence_from_simple_simplex_tree.cpp + * @example plain_homology.cpp * @example rips_multifield_persistence.cpp - * @example Skeleton_blocker_from_simplices.cpp - * @example Skeleton_blocker_iteration.cpp - * @example Skeleton_blocker_link.cpp - * @example Garland_heckbert.cpp - * @example Rips_contraction.cpp - * @example Random_bitmap_cubical_complex.cpp - * @example example_CGAL_3D_points_off_reader.cpp - * @example example_vector_double_points_off_reader.cpp - * @example example_CGAL_points_off_reader.cpp - * @example example_one_skeleton_rips_from_distance_matrix.cpp - * @example example_one_skeleton_rips_from_points.cpp - * @example example_rips_complex_from_csv_distance_matrix_file.cpp - * @example example_rips_complex_from_off_file.cpp - * @example persistence_intervals.cpp - * @example persistence_vectors.cpp - * @example persistence_heat_maps.cpp - * @example persistence_landscape_on_grid.cpp - * @example persistence_landscape.cpp - * @example example_basic.cpp - * @example example_with_perturb.cpp - * @example example_custom_distance.cpp - * @example example_choose_n_farthest_points.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 - * @example CoordGIC.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 CoordGIC.cpp * @example FuncGIC.cpp - * @example VoronoiGIC.cpp - * @example example_spatial_searching.cpp - * @example alpha_complex_3d_persistence.cpp - * @example alpha_complex_persistence.cpp - * @example Weighted_alpha_complex_3d_from_points.cpp - * @example bottleneck_distance.cpp - * @example weak_witness_persistence.cpp - * @example strong_witness_persistence.cpp - * @example cubical_complex_persistence.cpp - * @example periodic_cubical_complex_persistence.cpp - * @example off_file_from_shape_generator.cpp - * @example rips_distance_matrix_persistence.cpp - * @example rips_persistence.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/plot_landscapes_on_grid.cpp - * @example persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp * @example persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp - * @example persistence_landscapes_on_grid/average_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/compute_number_of_dominant_intervals.cpp * @example persistence_intervals/plot_persistence_Betti_numbers.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/compute_number_of_dominant_intervals.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/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_distance_of_persistence_heat_maps.cpp * @example persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp - * @example persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp - * @example persistence_heat_maps/average_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/plot_persistence_vectors.cpp - * @example persistence_vectors/compute_distance_of_persistence_vectors.cpp - * @example persistence_vectors/average_persistence_vectors.cpp - * @example persistence_vectors/create_persistence_vectors.cpp - * @example persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp - * @example persistence_landscapes/average_landscapes.cpp - * @example persistence_landscapes/compute_scalar_product_of_landscapes.cpp - * @example persistence_landscapes/create_landscapes.cpp + * @example persistence_heat_maps/average_persistence_heat_maps.cpp * @example persistence_landscapes/compute_distance_of_landscapes.cpp + * @example persistence_landscapes/compute_scalar_product_of_landscapes.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_landscape.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 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 5d40a48e..ef668dfb 100644 --- a/src/common/doc/installation.h +++ b/src/common/doc/installation.h @@ -88,9 +88,9 @@ 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="alpha_complex_from_off_8cpp-example.html"> + * \li <a href="_alpha_complex_from_off_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="alpha_complex_from_points_8cpp-example.html"> + * \li <a href="_alpha_complex_from_points_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_points.cpp</a> * \li <a href="alpha_complex_persistence_8cpp-example.html"> * Alpha_complex/alpha_complex_persistence.cpp</a> @@ -100,15 +100,15 @@ make \endverbatim * Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp.cpp</a> * \li <a href="bottleneck_basic_example_8cpp-example.html"> * Bottleneck_distance/bottleneck_basic_example.cpp</a> - * \li <a href="bottleneck_read_file_8cpp-example.html"> + * \li <a href="bottleneck_distance_8cpp-example.html"> * Bottleneck_distance/bottleneck_distance.cpp</a> - * \li <a href="coord_g_i_c_8cpp-example.html"> + * \li <a href="_coord_g_i_c_8cpp-example.html"> * Nerve_GIC/CoordGIC.cpp</a> - * \li <a href="func_g_i_c_8cpp-example.html"> + * \li <a href="_func_g_i_c_8cpp-example.html"> * Nerve_GIC/FuncGIC.cpp</a> - * \li <a href="nerve_8cpp-example.html"> + * \li <a href="_nerve_8cpp-example.html"> * Nerve_GIC/Nerve.cpp</a> - * \li <a href="voronoi_g_i_c_8cpp-example.html"> + * \li <a href="_voronoi_g_i_c_8cpp-example.html"> * Nerve_GIC/VoronoiGIC.cpp</a> * \li <a href="example_spatial_searching_8cpp-example.html"> * Spatial_searching/example_spatial_searching.cpp</a> @@ -122,10 +122,12 @@ make \endverbatim * Tangential_complex/example_basic.cpp</a> * \li <a href="example_with_perturb_8cpp-example.html"> * Tangential_complex/example_with_perturb.cpp</a> - * \li <a href="weighted_alpha_complex_3d_from_points_8cpp-example.html"> + * \li <a href="_weighted_alpha_complex_3d_from_points_8cpp-example.html"> * 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 @@ -134,15 +136,15 @@ make \endverbatim * * The following examples/utilities require the <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen</a> and will not be * built if Eigen is not installed: - * \li <a href="alpha_complex_from_off_8cpp-example.html"> + * \li <a href="_alpha_complex_from_off_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="alpha_complex_from_points_8cpp-example.html"> + * \li <a href="_alpha_complex_from_points_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_points.cpp</a> * \li <a href="alpha_complex_persistence_8cpp-example.html"> * Alpha_complex/alpha_complex_persistence.cpp</a> * \li <a href="alpha_complex_3d_persistence_8cpp-example.html"> * Alpha_complex/alpha_complex_3d_persistence.cpp</a> - * \li <a href="weighted_alpha_complex_3d_from_points_8cpp-example.html"> + * \li <a href="_weighted_alpha_complex_3d_from_points_8cpp-example.html"> * Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp</a> * \li <a href="alpha_rips_persistence_bottleneck_distance_8cpp-example.html"> * Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp.cpp</a> @@ -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® TBB</a> lets you easily write parallel @@ -179,27 +187,27 @@ make \endverbatim * Having Intel® TBB installed is recommended to parallelize and accelerate some GUDHI computations. * * The following examples/utilities are using Intel® TBB if installed: - * \li <a href="alpha_complex_from_off_8cpp-example.html"> + * \li <a href="_alpha_complex_from_off_8cpp-example.html"> * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="alpha_complex_from_points_8cpp-example.html"> + * \li <a href="_alpha_complex_from_points_8cpp-example.html"> * Alpha_complex/Alpha_complex_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="alpha_complex_persistence_8cpp-example.html"> * Alpha_complex/alpha_complex_persistence.cpp</a> - * \li <a href="bitmap_cubical_complex_8cpp-example.html"> + * \li <a href="cubical_complex_persistence_8cpp-example.html"> * Bitmap_cubical_complex/cubical_complex_persistence.cpp</a> - * \li <a href="bitmap_cubical_complex_periodic_boundary_conditions_8cpp-example.html"> + * \li <a href="periodic_cubical_complex_persistence_8cpp-example.html"> * Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp</a> - * \li <a href="random_bitmap_cubical_complex_8cpp-example.html"> + * \li <a href="_random_bitmap_cubical_complex_8cpp-example.html"> * Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp</a> - * \li <a href="coord_g_i_c_8cpp-example.html"> + * \li <a href="_coord_g_i_c_8cpp-example.html"> * Nerve_GIC/CoordGIC.cpp</a> - * \li <a href="func_g_i_c_8cpp-example.html"> + * \li <a href="_func_g_i_c_8cpp-example.html"> * Nerve_GIC/FuncGIC.cpp</a> - * \li <a href="nerve_8cpp-example.html"> + * \li <a href="_nerve_8cpp-example.html"> * Nerve_GIC/Nerve.cpp</a> - * \li <a href="voronoi_g_i_c_8cpp-example.html"> + * \li <a href="_voronoi_g_i_c_8cpp-example.html"> * Nerve_GIC/VoronoiGIC.cpp</a> * \li <a href="simple_simplex_tree_8cpp-example.html"> * Simplex_tree/simple_simplex_tree.cpp</a> 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 ≥ 3.1.0 + </td> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref coxeter_triangulation + </td> + </tr> +</table> + ### Tangential complex <table> diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index f534fc2a..96107cfe 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -41,14 +41,15 @@ function( add_gudhi_debug_info DEBUG_INFO ) endfunction( add_gudhi_debug_info ) if(PYTHONINTERP_FOUND) - if(PYBIND11_FOUND) + if(PYBIND11_FOUND AND CYTHON_FOUND) add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}") + # PyBind11 modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'clustering', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'datasets', ") - endif() - if(CYTHON_FOUND) + + # Cython modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_reader', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ") @@ -278,6 +279,7 @@ if(PYTHONINTERP_FOUND) # Some files for pip package file(COPY "introduction.rst" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") + file(COPY "pyproject.toml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") add_custom_command( OUTPUT gudhi.so @@ -298,35 +300,30 @@ if(PYTHONINTERP_FOUND) if(SCIPY_FOUND) if(SKLEARN_FOUND) if(OT_FOUND) - if(PYBIND11_FOUND) - if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") - # User warning - Sphinx is a static pages generator, and configured to work fine with user_version - # Images and biblio warnings because not found on developper version - if (GUDHI_PYTHON_PATH STREQUAL "src/python") - set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss") - endif() - # sphinx target requires gudhi.so, because conf.py reads gudhi version from it - add_custom_target(sphinx - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx - DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" - COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) - add_test(NAME sphinx_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) - # Set missing or not modules - set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") - else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 4.11.0") - set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - else(PYBIND11_FOUND) - message("++ Python documentation module will not be compiled because pybind11 was not found") + if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") + # User warning - Sphinx is a static pages generator, and configured to work fine with user_version + # Images and biblio warnings because not found on developper version + if (GUDHI_PYTHON_PATH STREQUAL "src/python") + set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss") + endif() + # sphinx target requires gudhi.so, because conf.py reads gudhi version from it + add_custom_target(sphinx + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx + DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" + COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) + add_test(NAME sphinx_py_test + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) + # Set missing or not modules + set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") + else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 4.11.0") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(PYBIND11_FOUND) + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) else(OT_FOUND) message("++ Python documentation module will not be compiled because POT was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") @@ -402,9 +399,7 @@ if(PYTHONINTERP_FOUND) COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py") - if (PYBIND11_FOUND) - add_gudhi_py_test(test_bottleneck_distance) - endif() + add_gudhi_py_test(test_bottleneck_distance) # Cover complex file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) @@ -518,14 +513,14 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_reader_utils) # Wasserstein - if(OT_FOUND AND PYBIND11_FOUND) + if(OT_FOUND) # EagerPy dependency because of enable_autodiff=True if(EAGERPY_FOUND) add_gudhi_py_test(test_wasserstein_distance) endif() + add_gudhi_py_test(test_wasserstein_barycenter) - endif() - if(OT_FOUND) + if(TORCH_FOUND AND TENSORFLOW_FOUND AND EAGERPY_FOUND) add_gudhi_py_test(test_wasserstein_with_tensors) endif() @@ -546,7 +541,7 @@ if(PYTHONINTERP_FOUND) endif() # Tomato - if(SCIPY_FOUND AND SKLEARN_FOUND AND PYBIND11_FOUND) + if(SCIPY_FOUND AND SKLEARN_FOUND) add_gudhi_py_test(test_tomato) endif() @@ -563,10 +558,10 @@ if(PYTHONINTERP_FOUND) # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") - else(CYTHON_FOUND) - message("++ Python module will not be compiled because cython was not found") + else(PYBIND11_FOUND AND CYTHON_FOUND) + message("++ Python module will not be compiled because cython and/or pybind11 was/were not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(CYTHON_FOUND) + endif(PYBIND11_FOUND AND CYTHON_FOUND) else(PYTHONINTERP_FOUND) message("++ Python module will not be compiled because no Python interpreter was found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index cb7700ce..35c344e3 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -359,7 +359,7 @@ Python Optimal Transport ------------------------ The :doc:`Wasserstein distance </wasserstein_distance_user>` -module requires `POT <https://pot.readthedocs.io/>`_, a library that provides +module requires `POT <https://pythonot.github.io/>`_, a library that provides several solvers for optimization problems related to Optimal Transport. PyTorch diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml new file mode 100644 index 00000000..a9fb4985 --- /dev/null +++ b/src/python/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools", "wheel", "numpy>=1.15.0", "cython", "pybind11"] +build-backend = "setuptools.build_meta" diff --git a/src/python/setup.py.in b/src/python/setup.py.in index c400b601..23746998 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -82,10 +82,10 @@ setup( }, description='The Gudhi library is an open source library for ' \ 'Computational Topology and Topological Data Analysis (TDA).', + data_files=[('.', ['./introduction.rst'])], long_description_content_type='text/x-rst', long_description=long_description, ext_modules = ext_modules, install_requires = ['numpy >= 1.15.0',], - setup_requires = ['cython','numpy >= 1.15.0','pybind11',], package_data={"": ["*.dll"], }, ) |