diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/Persistent_cohomology/doc/Intro_persistent_cohomology.h | 162 | ||||
-rw-r--r-- | src/Persistent_cohomology/example/CMakeLists.txt | 5 | ||||
-rw-r--r-- | src/Persistent_cohomology/example/custom_persistence_sort.cpp | 116 | ||||
-rw-r--r-- | src/Persistent_cohomology/example/plain_homology.cpp | 14 | ||||
-rw-r--r-- | src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h | 244 | ||||
-rw-r--r-- | src/Persistent_cohomology/test/CMakeLists.txt | 8 | ||||
-rw-r--r-- | src/Persistent_cohomology/test/betti_numbers_unit_test.cpp | 234 | ||||
-rw-r--r-- | src/common/doc/main_page.h | 7 |
8 files changed, 644 insertions, 146 deletions
diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h new file mode 100644 index 00000000..c8081cac --- /dev/null +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -0,0 +1,162 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef DOC_PERSISTENT_COHOMOLOGY_INTRO_PERSISTENT_COHOMOLOGY_H_ +#define DOC_PERSISTENT_COHOMOLOGY_INTRO_PERSISTENT_COHOMOLOGY_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +// needs namespace for Doxygen to link on classes +namespace persistent_cohomology { + +/** \defgroup persistent_cohomology Persistent Cohomology + + \author Clément Maria + + Computation of persistent cohomology using the algorithm of + \cite DBLP:journals/dcg/SilvaMV11 and \cite DBLP:journals/corr/abs-1208-5018 + and the Compressed Annotation Matrix + implementation of \cite DBLP:conf/esa/BoissonnatDM13 + + The theory of homology consists in attaching to a topological space a sequence of + (homology) groups, + capturing global topological features + like connected components, holes, cavities, etc. Persistent homology studies the evolution + -- birth, life and death -- of + these features when the topological space is changing. Consequently, the theory is essentially + composed of three elements: + topological spaces, their homology groups and an evolution scheme. + + <DT>Topological Spaces:</DT> + Topological spaces are represented by simplicial complexes. + Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>. + A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices + \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$ + on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$, + \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K} + \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ + is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is + a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq + f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$. + + We define the concept FilteredComplex which enumerates the requirements for a class + to represent a filtered complex from which persistent homology may be computed. + We use the vocabulary of simplicial complexes, but the concept + is valid for any type of cell complex. The main requirements + are the definition of: + \li type <CODE>Indexing_tag</CODE>, which is a model of the concept + <CODE>IndexingTag</CODE>, + describing the nature of the indexing scheme, + \li type Simplex_handle to manipulate simplices, + \li method <CODE>int dimension(Simplex_handle)</CODE> returning + the dimension of a simplex, + \li type and method <CODE>Boundary_simplex_range + boundary_simplex_range(Simplex_handle)</CODE> that returns + a range giving access to the codimension 1 subsimplices of the + input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the + definition of the operator \f$\partial\f$. The iterators have + value type <CODE>Simplex_handle</CODE>, + \li type and method + <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE> + that returns a range giving + access to all the simplices of the complex read in the order + assigned by the indexing scheme, + \li type and method + <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of + the filtration on the simplex represented by the handle. + + <DT>Homology:</DT> + For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>, + denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the + group of formal sums of + n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a + linear operator + \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$ + such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] = + \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$, + where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain + groups form a sequence: + + \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R}) + \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ } + \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ } \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f] + + of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms + \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$. + The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$ + for every \f$n > 0\f$ + and we define the homology groups: + + \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \ \partial_{n+1}\f] + + We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology + theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology. + + <DT>Indexing Scheme:</DT> + "Changing" a simplicial complex consists in applying a simplicial map. + An <EM>indexing scheme</EM> is a directed graph together with a traversal + order, such that two + consecutive nodes in the graph are connected by an arrow (either forward or backward). + The nodes represent simplicial complexes and the directed edges simplicial maps. + + From the computational point of view, there are two types of indexing schemes of + interest + in persistent homology: <EM>linear</EM> ones + \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet + \longrightarrow \bullet\f$ + in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , + and <EM>zigzag</EM> ones + \f$\bullet \longrightarrow \bullet \longleftarrow \cdots + \longrightarrow \bullet + \longleftarrow \bullet \f$ in zigzag persistent + homology \cite DBLP:journals/focm/CarlssonS10. + These indexing schemes have a natural left-to-right traversal order, and we + describe them with ranges and iterators. + In the current release of the Gudhi library, only the linear case is implemented. + + In the following, we consider the case where the indexing scheme is induced + by a filtration. + Ordering the simplices + by increasing filtration values (breaking ties so as a simplex appears after + its subsimplices of same filtration value) provides an indexing scheme. + +\section Examples + We provide several example files: run these examples with -h for details on their use, and read the README file. + +\li <CODE>rips_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram. + +\li <CODE>rips_multifield_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram +with a family of field coefficients. + +\li <CODE>performance_rips_persistence.cpp</CODE> provides timings for the construction of the Rips complex on a set of +points sampling a Klein bottle in \f$\mathbb{R}^5\f$ with a simplex tree, its conversion to a +Hasse diagram and the computation of persistent homology and multi-field persistent homology for the +different representations. + + \copyright GNU General Public License v3. + */ + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // DOC_PERSISTENT_COHOMOLOGY_INTRO_PERSISTENT_COHOMOLOGY_H_ diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 1f32da17..186a6c33 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -74,12 +74,17 @@ if(GMPXX_FOUND AND GMP_FOUND) add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp) target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + add_executable(custom_persistence_sort custom_persistence_sort.cpp) + target_link_libraries(custom_persistence_sort ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + if (TBB_FOUND) target_link_libraries(alpha_complex_persistence ${TBB_RELEASE_LIBRARY}) target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_RELEASE_LIBRARY}) + target_link_libraries(custom_persistence_sort ${TBB_RELEASE_LIBRARY}) endif() add_test(alpha_complex_persistence_2_0_45 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -m 0.45 -p 2) add_test(periodic_alpha_complex_3d_persistence_2_0 ${CMAKE_CURRENT_BINARY_DIR}/periodic_alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off 2 0) + add_test(custom_persistence_sort ${CMAKE_CURRENT_BINARY_DIR}/custom_persistence_sort) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") diff --git a/src/Persistent_cohomology/example/custom_persistence_sort.cpp b/src/Persistent_cohomology/example/custom_persistence_sort.cpp new file mode 100644 index 00000000..d591a0d0 --- /dev/null +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -0,0 +1,116 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 INRIA Saclay (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <fstream> // for std::ofstream +#include <algorithm> // for std::sort + + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Point = Kernel::Point_d; +using Alpha_complex = Gudhi::alphacomplex::Alpha_complex<Kernel>; + +std::vector<Point> random_points() { + // Instanciate a random point generator + CGAL::Random rng(0); + + // Generate "points_number" random points in a vector + std::vector<Point> points; + + // Generates 1000 random 3D points on a sphere of radius 4.0 + CGAL::Random_points_on_sphere_d<Point> rand_outside(3, 4.0, rng); + CGAL::cpp11::copy_n(rand_outside, 1000, std::back_inserter(points)); + // Generates 2000 random 3D points in a sphere of radius 3.0 + CGAL::Random_points_in_ball_d<Point> rand_inside(3, 3.0, rng); + CGAL::cpp11::copy_n(rand_inside, 2000, std::back_inserter(points)); + + return points; +} + +/* + * Compare two intervals by dimension, then by length. + */ +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Alpha_complex * sc) + : sc_(sc) { } + + template<typename Persistent_interval> + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2))) + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + else + return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); + } + Alpha_complex* sc_; +}; + +int main(int argc, char **argv) { + std::vector<Point> points = random_points(); + + // Alpha complex persistence computation from generated points + Alpha_complex alpha_complex_from_points(points, 0.6); + + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Alpha_complex, + Gudhi::persistent_cohomology::Field_Zp >; + Persistent_cohomology pcoh(alpha_complex_from_points); + + // initializes the coefficient field for homology - Z/3Z + pcoh.init_coefficients(3); + pcoh.compute_persistent_cohomology(0.2); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(&alpha_complex_from_points); + auto persistent_pairs = pcoh.get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + for (auto pair : persistent_pairs) { + std::cout << alpha_complex_from_points.dimension(get<0>(pair)) << " " + << alpha_complex_from_points.filtration(get<0>(pair)) << " " + << alpha_complex_from_points.filtration(get<1>(pair)) << std::endl; + } + + // Persistent Betti numbers + std::cout << "The persistent Betti numbers in interval [0.40, 0.41] are : "; + for (int dim = 0; dim < alpha_complex_from_points.dimension(); dim++) + std::cout << "b" << dim << " = " << pcoh.persistent_betti_number(dim, 0.40, 0.41) << " ; "; + std::cout << std::endl; + + // Betti numbers + std::vector<int> betti_numbers = pcoh.betti_numbers(); + std::cout << "The Betti numbers are : "; + for (std::size_t i = 0; i < betti_numbers.size(); i++) + std::cout << "b" << i << " = " << betti_numbers[i] << " ; "; + std::cout << std::endl; + + return 0; +} + diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp index 0a692717..2156ffd9 100644 --- a/src/Persistent_cohomology/example/plain_homology.cpp +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -24,6 +24,7 @@ #include <gudhi/Persistent_cohomology.h> #include <iostream> +#include <vector> using namespace Gudhi; @@ -76,9 +77,16 @@ int main() { // Print the result. The format is, on each line: 2 dim 0 inf // where 2 represents the field, dim the dimension of the feature. - // 2 0 0 inf - // 2 0 0 inf - // 2 1 0 inf + // 2 0 0 inf + // 2 0 0 inf + // 2 1 0 inf // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. pcoh.output_diagram(); + + // Print the Betti numbers are b0=2 and b1=1. + std::cout << std::endl; + std::cout << "The Betti numbers are : "; + for (int i = 0; i < st.dimension(); i++) + std::cout << "b" << i << " = " << pcoh.betti_number(i) << " ; "; + std::cout << std::endl; } diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 1b86f1f9..efa9a367 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -46,136 +46,10 @@ namespace Gudhi { namespace persistent_cohomology { -/** \defgroup persistent_cohomology Persistent Cohomology - * - \author Clément Maria - - Computation of persistent cohomology using the algorithm of - \cite DBLP:journals/dcg/SilvaMV11 and \cite DBLP:journals/corr/abs-1208-5018 - and the Compressed Annotation Matrix - implementation of \cite DBLP:conf/esa/BoissonnatDM13 - - The theory of homology consists in attaching to a topological space a sequence of - (homology) groups, - capturing global topological features - like connected components, holes, cavities, etc. Persistent homology studies the evolution - -- birth, life and death -- of - these features when the topological space is changing. Consequently, the theory is essentially - composed of three elements: - topological spaces, their homology groups and an evolution scheme. - - <DT>Topological Spaces:</DT> - Topological spaces are represented by simplicial complexes. - Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>. - A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices - \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$ - on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$, - \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K} - \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ - is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is - a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq - f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$. - - We define the concept FilteredComplex which enumerates the requirements for a class - to represent a filtered complex from which persistent homology may be computed. - We use the vocabulary of simplicial complexes, but the concept - is valid for any type of cell complex. The main requirements - are the definition of: - \li type <CODE>Indexing_tag</CODE>, which is a model of the concept - <CODE>IndexingTag</CODE>, - describing the nature of the indexing scheme, - \li type Simplex_handle to manipulate simplices, - \li method <CODE>int dimension(Simplex_handle)</CODE> returning - the dimension of a simplex, - \li type and method <CODE>Boundary_simplex_range - boundary_simplex_range(Simplex_handle)</CODE> that returns - a range giving access to the codimension 1 subsimplices of the - input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the - definition of the operator \f$\partial\f$. The iterators have - value type <CODE>Simplex_handle</CODE>, - \li type and method - <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE> - that returns a range giving - access to all the simplices of the complex read in the order - assigned by the indexing scheme, - \li type and method - <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of - the filtration on the simplex represented by the handle. - - <DT>Homology:</DT> - For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>, - denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the - group of formal sums of - n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a - linear operator - \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$ - such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] = - \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$, - where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain - groups form a sequence: - - \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R}) - \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ } - \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ } \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f] - - of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms - \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$. - The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$ - for every \f$n > 0\f$ - and we define the homology groups: - - \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \ \partial_{n+1}\f] - - We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology - theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology. - - <DT>Indexing Scheme:</DT> - "Changing" a simplicial complex consists in applying a simplicial map. - An <EM>indexing scheme</EM> is a directed graph together with a traversal - order, such that two - consecutive nodes in the graph are connected by an arrow (either forward or backward). - The nodes represent simplicial complexes and the directed edges simplicial maps. - - From the computational point of view, there are two types of indexing schemes of - interest - in persistent homology: <EM>linear</EM> ones - \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet - \longrightarrow \bullet\f$ - in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , - and <EM>zigzag</EM> ones - \f$\bullet \longrightarrow \bullet \longleftarrow \cdots - \longrightarrow \bullet - \longleftarrow \bullet \f$ in zigzag persistent - homology \cite DBLP:journals/focm/CarlssonS10. - These indexing schemes have a natural left-to-right traversal order, and we - describe them with ranges and iterators. - In the current release of the Gudhi library, only the linear case is implemented. - - In the following, we consider the case where the indexing scheme is induced - by a filtration. - Ordering the simplices - by increasing filtration values (breaking ties so as a simplex appears after - its subsimplices of same filtration value) provides an indexing scheme. - -\section Examples - We provide several example files: run these examples with -h for details on their use, and read the README file. - -\li <CODE>rips_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram. - -\li <CODE>rips_multifield_persistence.cpp</CODE> computes the Rips complex of a point cloud and its persistence diagram -with a family of field coefficients. - -\li <CODE>performance_rips_persistence.cpp</CODE> provides timings for the construction of the Rips complex on a set of -points sampling a Klein bottle in \f$\mathbb{R}^5\f$ with a simplex tree, its conversion to a -Hasse diagram and the computation of persistent homology and multi-field persistent homology for the -different representations. - - \copyright GNU General Public License v3. - @{ - */ - /** \brief Computes the persistent cohomology of a filtered complex. * + * \ingroup persistent_cohomology + * * The computation is implemented with a Compressed Annotation Matrix * (CAM)\cite DBLP:conf/esa/BoissonnatDM13, * and is adapted to the computation of Multi-Field Persistent Homology (MF) @@ -184,8 +58,8 @@ different representations. * \implements PersistentHomology * */ -// Memory allocation policy: classic, use a mempool, etc.*/ -template<class FilteredComplex, class CoefficientField> // to do mem allocation policy: classic, mempool, etc. +// TODO(CM): Memory allocation policy: classic, use a mempool, etc. +template<class FilteredComplex, class CoefficientField> class Persistent_cohomology { public: typedef FilteredComplex Complex_ds; @@ -194,7 +68,7 @@ class Persistent_cohomology { typedef typename Complex_ds::Simplex_handle Simplex_handle; typedef typename Complex_ds::Filtration_value Filtration_value; typedef typename CoefficientField::Element Arith_element; -// Compressed Annotation Matrix types: + // Compressed Annotation Matrix types: // Column type typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column; // contains 1 set_hook // Cell type @@ -206,15 +80,15 @@ class Persistent_cohomology { typedef boost::intrusive::set<Column, boost::intrusive::constant_time_size<false> > Cam; -// Sparse column type for the annotation of the boundary of an element. + // Sparse column type for the annotation of the boundary of an element. typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type; -// Persistent interval type. The Arith_element field is used for the multi-field framework. + // Persistent interval type. The Arith_element field is used for the multi-field framework. typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval; /** \brief Initializes the Persistent_cohomology class. * * @param[in] cpx Complex for which the persistent homology is computed. - cpx is a model of FilteredComplex + * cpx is a model of FilteredComplex */ explicit Persistent_cohomology(Complex_ds& cpx) : cpx_(&cpx), @@ -243,7 +117,7 @@ class Persistent_cohomology { /** \brief Initializes the Persistent_cohomology class. * * @param[in] cpx Complex for which the persistent homology is compiuted. - cpx is a model of FilteredComplex + * cpx is a model of FilteredComplex * * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the * complex is computed. If false, it is ignored. Default is false. @@ -332,7 +206,7 @@ class Persistent_cohomology { persistent_pairs_.emplace_back( cpx_->simplex(zero_idx.second), cpx_->null_simplex(), coeff_field_.characteristic()); } -// Compute infinite interval of dimension > 0 + // Compute infinite interval of dimension > 0 for (auto cocycle : transverse_idx_) { persistent_pairs_.emplace_back( cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_); @@ -434,14 +308,14 @@ class Persistent_cohomology { // Find its annotation vector curr_col = ds_repr_[dsets_.find_set(key)]; if (curr_col != NULL) { // and insert it in annotations_in_boundary with multyiplicative factor "sign". - annotations_in_boundary.emplace_back(curr_col, sign); + annotations_in_boundary.emplace_back(curr_col, sign); } } sign = -sign; } // Place identical annotations consecutively so we can easily sum their multiplicities. std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(), - [](annotation_t const& a, annotation_t const& b) { return a.first < b.first; }); + [](annotation_t const& a, annotation_t const& b) { return a.first < b.first; }); // Sum the annotations with multiplicity, using a map<key,coeff> // to represent a sparse vector. @@ -451,7 +325,7 @@ class Persistent_cohomology { Column* col = ann_it->first; int mult = ann_it->second; while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) { - mult += ann_it->second; + mult += ann_it->second; } // The following test is just a heuristic, it is not required, and it is fine that is misses p == 0. if (mult != coeff_field_.additive_identity()) { // For all columns in the boundary, @@ -692,7 +566,6 @@ class Persistent_cohomology { * feature exists in homology with Z/piZ coefficients. */ void output_diagram(std::ostream& ostream = std::cout) { - cmp_intervals_by_length cmp(cpx_); std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp); bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity; @@ -720,12 +593,99 @@ class Persistent_cohomology { } } + /** @brief Returns Betti numbers. + * @return A vector of persistent Betti numbers. + */ + std::vector<int> betti_numbers() const { + // Init Betti numbers vector with zeros until Simplicial complex dimension + std::vector<int> betti_numbers(cpx_->dimension(), 0); + + for (auto pair : persistent_pairs_) { + // Count never ended persistence intervals + if (cpx_->null_simplex() == get<1>(pair)) { + // Increment corresponding betti number + betti_numbers[cpx_->dimension(get<0>(pair))] += 1; + } + } + return betti_numbers; + } + + /** @brief Returns the Betti number of the dimension passed by parameter. + * @param[in] dimension The Betti number dimension to get. + * @return Betti number of the given dimension + * + */ + int betti_number(int dimension) const { + int betti_number = 0; + + for (auto pair : persistent_pairs_) { + // Count never ended persistence intervals + if (cpx_->null_simplex() == get<1>(pair)) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + // Increment betti number found + ++betti_number; + } + } + } + return betti_number; + } + + /** @brief Returns the persistent Betti numbers. + * @param[in] from The persistence birth limit to be added in the number \f$(persistent birth \leq from)\f$. + * @param[in] to The persistence death limit to be added in the number \f$(persistent death > from)\f$. + * @return A vector of persistent Betti numbers. + */ + std::vector<int> persistent_betti_numbers(Filtration_value from, Filtration_value to) const { + // Init Betti numbers vector with zeros until Simplicial complex dimension + std::vector<int> betti_numbers(cpx_->dimension(), 0); + + for (auto pair : persistent_pairs_) { + // Count persistence intervals that covers the given interval + if (cpx_->filtration(get<0>(pair)) <= from && cpx_->filtration(get<1>(pair)) > to) { + // Increment corresponding betti number + betti_numbers[cpx_->dimension(get<0>(pair))] += 1; + } + } + return betti_numbers; + } + + /** @brief Returns the persistent Betti number of the dimension passed by parameter. + * @param[in] dimension The Betti number dimension to get. + * @param[in] from The persistence birth limit to be added in the number \f$(persistent birth \leq from)\f$. + * @param[in] to The persistence death limit to be added in the number \f$(persistent death > from)\f$. + * @return Persistent Betti number of the given dimension + */ + int persistent_betti_number(int dimension, Filtration_value from, Filtration_value to) const { + int betti_number = 0; + + for (auto pair : persistent_pairs_) { + // Count persistence intervals that covers the given interval + if (cpx_->filtration(get<0>(pair)) <= from && cpx_->filtration(get<1>(pair)) > to) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + // Increment betti number found + ++betti_number; + } + } + } + return betti_number; + } + + /** @brief Returns the persistent pairs. + * @return Persistent pairs + * + */ + const std::vector<Persistent_interval>& get_persistent_pairs() const { + return persistent_pairs_; + } + private: /* * Structure representing a cocycle. */ struct cocycle { - cocycle() { + cocycle() + : row_(nullptr), + characteristics_() { } cocycle(Arith_element characteristics, Hcell * row) : row_(row), @@ -766,8 +726,6 @@ class Persistent_cohomology { Simple_object_pool<Cell> cell_pool_; }; -/** @} */ // end defgroup persistent_cohomology - } // namespace persistent_cohomology } // namespace Gudhi diff --git a/src/Persistent_cohomology/test/CMakeLists.txt b/src/Persistent_cohomology/test/CMakeLists.txt index 459cc000..a034031a 100644 --- a/src/Persistent_cohomology/test/CMakeLists.txt +++ b/src/Persistent_cohomology/test/CMakeLists.txt @@ -12,8 +12,11 @@ endif() add_executable ( PersistentCohomologyUT persistent_cohomology_unit_test.cpp ) target_link_libraries(PersistentCohomologyUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( BettiNumbersUT betti_numbers_unit_test.cpp ) +target_link_libraries(BettiNumbersUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) if (TBB_FOUND) target_link_libraries(PersistentCohomologyUT ${TBB_RELEASE_LIBRARY}) + target_link_libraries(BettiNumbersUT ${TBB_RELEASE_LIBRARY}) endif() # Unitary tests @@ -23,6 +26,11 @@ add_test(NAME PersistentCohomologyUT # XML format for Jenkins xUnit plugin --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/PersistentCohomologyUT.xml --log_level=test_suite --report_level=no) +add_test(NAME BettiNumbersUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BettiNumbersUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BettiNumbersUT.xml --log_level=test_suite --report_level=no) + if(GMPXX_FOUND AND GMP_FOUND) add_executable ( PersistentCohomologyMultiFieldUT persistent_cohomology_unit_test_multi_field.cpp ) target_link_libraries(PersistentCohomologyMultiFieldUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp new file mode 100644 index 00000000..a4e00b45 --- /dev/null +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -0,0 +1,234 @@ +#include <iostream> +#include <string> +#include <algorithm> +#include <utility> // std::pair, std::make_pair +#include <cmath> // float comparison +#include <limits> + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "betti_numbers" +#include <boost/test/unit_test.hpp> + +#include "gudhi/Simplex_tree.h" +#include "gudhi/Persistent_cohomology.h" + +struct MyOptions : Gudhi::Simplex_tree_options_full_featured { + // Implicitly use 0 as filtration value for all simplices + static const bool store_filtration = false; + // The persistence algorithm needs this + static const bool store_key = true; + // I have few vertices + typedef short Vertex_handle; +}; + +using Mini_simplex_tree = Gudhi::Simplex_tree<MyOptions>; +using Mini_st_persistence = + Gudhi::persistent_cohomology::Persistent_cohomology<Mini_simplex_tree, Gudhi::persistent_cohomology::Field_Zp>; + +/* + * Compare two intervals by dimension, then by length. + */ +template<class Simplicial_complex> +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Simplicial_complex * sc) + : sc_(sc) { } + + template<typename Persistent_interval> + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2))) + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + else + return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); + } + Simplicial_complex* sc_; +}; + +BOOST_AUTO_TEST_CASE( plain_homology_betti_numbers ) +{ + Mini_simplex_tree st; + + /* Complex to build. */ + /* 1 4 */ + /* o---o */ + /* /3\ / */ + /* o---o o */ + /* 2 0 5 */ + const short tetra0123[] = {0, 1, 2, 3}; + const short edge04[] = {0, 4}; + const short edge14[] = {1, 4}; + const short vertex5[] = {5}; + st.insert_simplex_and_subfaces(tetra0123); + st.insert_simplex_and_subfaces(edge04); + st.insert_simplex(edge14); + st.insert_simplex(vertex5); + // FIXME: Remove this line + st.set_dimension(3); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Class for homology computation + Mini_st_persistence pcoh(st); + + // Initialize the coefficient field Z/3Z for homology + pcoh.init_coefficients(3); + + // Compute the persistence diagram of the complex + pcoh.compute_persistent_cohomology(); + + // Print the result. The format is, on each line: 2 dim 0 inf + // where 2 represents the field, dim the dimension of the feature. + // 2 0 0 inf + // 2 0 0 inf + // 2 1 0 inf + // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. + + BOOST_CHECK(pcoh.betti_number(0) == 2); + BOOST_CHECK(pcoh.betti_number(1) == 1); + BOOST_CHECK(pcoh.betti_number(2) == 0); + + std::vector<int> bns = pcoh.betti_numbers(); + BOOST_CHECK(bns.size() == 3); + BOOST_CHECK(bns[0] == 2); + BOOST_CHECK(bns[1] == 1); + BOOST_CHECK(bns[2] == 0); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length<Mini_simplex_tree> cmp(&st); + auto persistent_pairs = pcoh.get_persistent_pairs(); + + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + + BOOST_CHECK(persistent_pairs.size() == 3); + // persistent_pairs[0] = 2 1 0 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[0])) == 1); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[0])) == 0); + BOOST_CHECK(get<1>(persistent_pairs[0]) == st.null_simplex()); + + // persistent_pairs[1] = 2 0 0 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[1])) == 0); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[1])) == 0); + BOOST_CHECK(get<1>(persistent_pairs[1]) == st.null_simplex()); + + // persistent_pairs[2] = 2 0 0 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 0); + BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); +} + +using Simplex_tree = Gudhi::Simplex_tree<>; +using St_persistence = + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp>; + +BOOST_AUTO_TEST_CASE( betti_numbers ) +{ + Simplex_tree st; + + /* Complex to build. */ + /* 1 4 */ + /* o---o */ + /* /3\ / */ + /* o---o o */ + /* 2 0 5 */ + const short tetra0123[] = {0, 1, 2, 3}; + const short edge04[] = {0, 4}; + const short edge14[] = {1, 4}; + const short vertex5[] = {5}; + st.insert_simplex_and_subfaces(tetra0123, 4.0); + st.insert_simplex_and_subfaces(edge04, 2.0); + st.insert_simplex(edge14, 2.0); + st.insert_simplex(vertex5, 1.0); + // FIXME: Remove this line + st.set_dimension(3); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Class for homology computation + St_persistence pcoh(st); + + // Initialize the coefficient field Z/3Z for homology + pcoh.init_coefficients(3); + + // Compute the persistence diagram of the complex + pcoh.compute_persistent_cohomology(); + + // Check the Betti numbers are b0=2, b1=1 and b2=0. + BOOST_CHECK(pcoh.betti_number(0) == 2); + BOOST_CHECK(pcoh.betti_number(1) == 1); + BOOST_CHECK(pcoh.betti_number(2) == 0); + + // Check the Betti numbers are b0=2, b1=1 and b2=0. + std::vector<int> bns = pcoh.betti_numbers(); + BOOST_CHECK(bns.size() == 3); + BOOST_CHECK(bns[0] == 2); + BOOST_CHECK(bns[1] == 1); + BOOST_CHECK(bns[2] == 0); + + // Check the persistent Betti numbers in [4., 10.] are b0=2, b1=1 and b2=0. + BOOST_CHECK(pcoh.persistent_betti_number(0, 4., 10.) == 2); + BOOST_CHECK(pcoh.persistent_betti_number(1, 4., 10.) == 1); + BOOST_CHECK(pcoh.persistent_betti_number(2, 4., 10.) == 0); + + // Check the persistent Betti numbers in [2., 100.] are b0=2, b1=0 and b2=0. + BOOST_CHECK(pcoh.persistent_betti_number(0, 2., 100.) == 2); + BOOST_CHECK(pcoh.persistent_betti_number(1, 2., 100.) == 0); + BOOST_CHECK(pcoh.persistent_betti_number(2, 2., 100.) == 0); + + // Check the persistent Betti numbers in [1., 1000.] are b0=1, b1=0 and b2=0. + BOOST_CHECK(pcoh.persistent_betti_number(0, 1., 1000.) == 1); + BOOST_CHECK(pcoh.persistent_betti_number(1, 1., 1000.) == 0); + BOOST_CHECK(pcoh.persistent_betti_number(2, 1., 1000.) == 0); + + // Check the persistent Betti numbers in [.9, 1000.] are b0=0, b1=0 and b2=0. + BOOST_CHECK(pcoh.persistent_betti_number(0, .9, 1000.) == 0); + BOOST_CHECK(pcoh.persistent_betti_number(1, .9, 1000.) == 0); + BOOST_CHECK(pcoh.persistent_betti_number(2, .9, 1000.) == 0); + + // Check the persistent Betti numbers in [4.1, 10000.] are b0=2, b1=1 and b2=0. + bns = pcoh.persistent_betti_numbers(4.1, 10000.); + BOOST_CHECK(bns[0] == 2); + BOOST_CHECK(bns[1] == 1); + BOOST_CHECK(bns[2] == 0); + + // Check the persistent Betti numbers in [2.1, 100000.] are b0=2, b1=0 and b2=0. + bns = pcoh.persistent_betti_numbers(2.1, 100000.); + BOOST_CHECK(bns[0] == 2); + BOOST_CHECK(bns[1] == 0); + BOOST_CHECK(bns[2] == 0); + + // Check the persistent Betti numbers in [1.1, 1000000.] are b0=1, b1=0 and b2=0. + bns = pcoh.persistent_betti_numbers(1.1, 1000000.); + BOOST_CHECK(bns[0] == 1); + BOOST_CHECK(bns[1] == 0); + BOOST_CHECK(bns[2] == 0); + + // Check the persistent Betti numbers in [.1, 10000000.] are b0=0, b1=0 and b2=0. + bns = pcoh.persistent_betti_numbers(.1, 10000000.); + BOOST_CHECK(bns[0] == 0); + BOOST_CHECK(bns[1] == 0); + BOOST_CHECK(bns[2] == 0); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length<Simplex_tree> cmp(&st); + auto persistent_pairs = pcoh.get_persistent_pairs(); + + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + + BOOST_CHECK(persistent_pairs.size() == 3); + // persistent_pairs[0] = 2 1 4 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[0])) == 1); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[0])) == 4); + BOOST_CHECK(get<1>(persistent_pairs[0]) == st.null_simplex()); + + // persistent_pairs[1] = 2 0 2 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[1])) == 0); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[1])) == 2); + BOOST_CHECK(get<1>(persistent_pairs[1]) == st.null_simplex()); + + // persistent_pairs[2] = 2 0 1 inf + BOOST_CHECK(st.dimension(get<0>(persistent_pairs[2])) == 0); + BOOST_CHECK(st.filtration(get<0>(persistent_pairs[2])) == 1); + BOOST_CHECK(get<1>(persistent_pairs[2]) == st.null_simplex()); +} diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 19dece4b..2b9f2f40 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -209,6 +209,8 @@ * Persistent_cohomology/alpha_complex_persistence.cpp</a> * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html"> * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a> + * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> + * Persistent_cohomology/custom_persistence_sort.cpp</a> * * \subsection eigen3 Eigen3: * <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen3</a> is a C++ template library for linear algebra: @@ -224,6 +226,8 @@ * Persistent_cohomology/alpha_complex_persistence.cpp</a> * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html"> * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a> + * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> + * Persistent_cohomology/custom_persistence_sort.cpp</a> * * \subsection tbb Threading Building Blocks: * <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel® TBB</a> lets you easily write parallel @@ -271,6 +275,8 @@ * Persistent_cohomology/rips_persistence.cpp</a> * \li <a href="_persistent_cohomology_2periodic_alpha_complex_3d_persistence_8cpp-example.html"> * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp</a> + * \li <a href="_persistent_cohomology_2custom_persistence_sort_8cpp-example.html"> + * Persistent_cohomology/custom_persistence_sort.cpp</a> * * \subsection demos Demos and examples * To build the demos and examples, run the following commands in a terminal: @@ -325,6 +331,7 @@ make \endverbatim * @example Persistent_cohomology/plain_homology.cpp * @example Persistent_cohomology/rips_multifield_persistence.cpp * @example Persistent_cohomology/rips_persistence.cpp + * @example Persistent_cohomology/custom_persistence_sort.cpp * @example Simplex_tree/mini_simplex_tree.cpp * @example Simplex_tree/simple_simplex_tree.cpp * @example Simplex_tree/simplex_tree_from_alpha_shapes_3.cpp |