From b16152118574b9f2127df3b7f4495f75f3b079c1 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 19 May 2016 15:27:00 +0000 Subject: Seperate doc from code in Persistent_cohomology module. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1184 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 90b406b11cfa6356243bb4e881c8c7db06f14d6e --- .../doc/Intro_persistent_cohomology.h | 162 +++++++++++++++++++++ .../include/gudhi/Persistent_cohomology.h | 136 +---------------- 2 files changed, 166 insertions(+), 132 deletions(-) create mode 100644 src/Persistent_cohomology/doc/Intro_persistent_cohomology.h (limited to 'src') 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 . + */ + +#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. + +
Topological Spaces:
+ Topological spaces are represented by simplicial complexes. + Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of vertices. + A simplex \f$\sigma\f$ is a subset of vertices + \f$\sigma \subseteq V\f$. A simplicial complex \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 filtration 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 Indexing_tag, which is a model of the concept + IndexingTag, + describing the nature of the indexing scheme, + \li type Simplex_handle to manipulate simplices, + \li method int dimension(Simplex_handle) returning + the dimension of a simplex, + \li type and method Boundary_simplex_range + boundary_simplex_range(Simplex_handle) 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 Simplex_handle, + \li type and method + Filtration_simplex_range filtration_simplex_range () + 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 + Filtration_value filtration (Simplex_handle) that returns the value of + the filtration on the simplex represented by the handle. + +
Homology:
+ For a ring \f$\mathcal{R}\f$, the group of n-chains, + 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 boundary operator 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. + +
Indexing Scheme:
+ "Changing" a simplicial complex consists in applying a simplicial map. + An indexing scheme 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: linear ones + \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet + \longrightarrow \bullet\f$ + in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , + and zigzag 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 rips_persistence.cpp computes the Rips complex of a point cloud and its persistence diagram. + +\li rips_multifield_persistence.cpp computes the Rips complex of a point cloud and its persistence diagram +with a family of field coefficients. + +\li performance_rips_persistence.cpp 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/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 1b86f1f9..cd1bd438 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. - -
Topological Spaces:
- Topological spaces are represented by simplicial complexes. - Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of vertices. - A simplex \f$\sigma\f$ is a subset of vertices - \f$\sigma \subseteq V\f$. A simplicial complex \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 filtration 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 Indexing_tag, which is a model of the concept - IndexingTag, - describing the nature of the indexing scheme, - \li type Simplex_handle to manipulate simplices, - \li method int dimension(Simplex_handle) returning - the dimension of a simplex, - \li type and method Boundary_simplex_range - boundary_simplex_range(Simplex_handle) 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 Simplex_handle, - \li type and method - Filtration_simplex_range filtration_simplex_range () - 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 - Filtration_value filtration (Simplex_handle) that returns the value of - the filtration on the simplex represented by the handle. - -
Homology:
- For a ring \f$\mathcal{R}\f$, the group of n-chains, - 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 boundary operator 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. - -
Indexing Scheme:
- "Changing" a simplicial complex consists in applying a simplicial map. - An indexing scheme 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: linear ones - \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet - \longrightarrow \bullet\f$ - in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , - and zigzag 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 rips_persistence.cpp computes the Rips complex of a point cloud and its persistence diagram. - -\li rips_multifield_persistence.cpp computes the Rips complex of a point cloud and its persistence diagram -with a family of field coefficients. - -\li performance_rips_persistence.cpp 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 // to do mem allocation policy: classic, mempool, etc. +// TODO(CM): Memory allocation policy: classic, use a mempool, etc. +template class Persistent_cohomology { public: typedef FilteredComplex Complex_ds; @@ -766,8 +640,6 @@ class Persistent_cohomology { Simple_object_pool cell_pool_; }; -/** @} */ // end defgroup persistent_cohomology - } // namespace persistent_cohomology } // namespace Gudhi -- cgit v1.2.3 From 6dae848dd4672fae2c9a6b5ae9059b44855e26a3 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 24 May 2016 15:17:14 +0000 Subject: get_persistence and Betti numbers git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1190 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 6da18e680111a9a1ad7753779a94a1e8ce49702c --- .../example/plain_homology.cpp | 51 ++++++++++++++++ .../include/gudhi/Persistent_cohomology.h | 67 +++++++++++++++++----- 2 files changed, 104 insertions(+), 14 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp index 0a692717..9e5adb5d 100644 --- a/src/Persistent_cohomology/example/plain_homology.cpp +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -24,6 +24,7 @@ #include #include +#include using namespace Gudhi; @@ -41,6 +42,24 @@ struct MyOptions : Simplex_tree_options_full_featured { }; typedef Simplex_tree ST; + /* + * Compare two intervals by dimension, then by length. + */ + struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(ST * sc) + : sc_(sc) { + } + template + 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))); + } + ST* sc_; + }; + int main() { ST st; @@ -81,4 +100,36 @@ int main() { // 2 1 0 inf // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. pcoh.output_diagram(); + + + // ******************************************************** + // get_persistence + // ******************************************************** + std::cout << std::endl; + std::cout << std::endl; + + // First version + std::vector 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; + + // Second version + 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; + + // Get persistence + cmp_intervals_by_dim_then_length cmp(&st); + 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 << st.dimension(get<0>(pair)) << " " + << st.filtration(get<0>(pair)) << " " + << st.filtration(get<1>(pair)) << std::endl; + } + + } diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index cd1bd438..9fc3bf2e 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -83,7 +83,7 @@ class Persistent_cohomology { // Sparse column type for the annotation of the boundary of an element. typedef std::vector > A_ds_type; // Persistent interval type. The Arith_element field is used for the multi-field framework. - typedef std::tuple Persistent_interval; + typedef std::tuple Persistent_interval; /** \brief Initializes the Persistent_cohomology class. * @@ -199,17 +199,17 @@ class Persistent_cohomology { if (ds_parent_[key] == key // root of its tree && zero_cocycles_.find(key) == zero_cocycles_.end()) { persistent_pairs_.emplace_back( - cpx_->simplex(key), cpx_->null_simplex(), coeff_field_.characteristic()); + cpx_->simplex(key), cpx_->null_simplex()/*, coeff_field_.characteristic()*/); } } for (auto zero_idx : zero_cocycles_) { persistent_pairs_.emplace_back( - cpx_->simplex(zero_idx.second), cpx_->null_simplex(), coeff_field_.characteristic()); + cpx_->simplex(zero_idx.second), cpx_->null_simplex()/*, coeff_field_.characteristic()*/); } // 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_); + cpx_->simplex(cocycle.first), cpx_->null_simplex()/*, cocycle.second.characteristics_*/); } } @@ -250,7 +250,7 @@ class Persistent_cohomology { < cpx_->filtration(cpx_->simplex(idx_coc_v))) { // Kill cocycle [idx_coc_v], which is younger. if (interval_length_policy(cpx_->simplex(idx_coc_v), sigma)) { persistent_pairs_.emplace_back( - cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic()); + cpx_->simplex(idx_coc_v), sigma/*, coeff_field_.characteristic()*/); } // Maintain the index of the 0-cocycle alive. if (kv != idx_coc_v) { @@ -265,7 +265,7 @@ class Persistent_cohomology { } else { // Kill cocycle [idx_coc_u], which is younger. if (interval_length_policy(cpx_->simplex(idx_coc_u), sigma)) { persistent_pairs_.emplace_back( - cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic()); + cpx_->simplex(idx_coc_u), sigma/*, coeff_field_.characteristic()*/); } // Maintain the index of the 0-cocycle alive. if (ku != idx_coc_u) { @@ -308,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 // to represent a sparse vector. @@ -325,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, @@ -427,7 +427,7 @@ class Persistent_cohomology { if (interval_length_policy(cpx_->simplex(death_key), sigma)) { persistent_pairs_.emplace_back(cpx_->simplex(death_key) // creator , sigma // destructor - , charac); // fields + /*, charac*/); // fields } auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row. @@ -566,17 +566,16 @@ 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::has_infinity; for (auto pair : persistent_pairs_) { // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA) if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits::infinity()) { - ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " " + ostream << /*get<2>(pair) <<*/ " " << cpx_->dimension(get<0>(pair)) << " " << cpx_->filtration(get<0>(pair)) << " inf " << std::endl; } else { - ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " " + ostream << /*get<2>(pair) <<*/ " " << cpx_->dimension(get<0>(pair)) << " " << cpx_->filtration(get<0>(pair)) << " " << cpx_->filtration(get<1>(pair)) << " " << std::endl; } @@ -594,12 +593,52 @@ class Persistent_cohomology { } } + /** @brief Returns Betti numbers. + * + */ + const std::vector betti_numbers() const { + // Init Betti numbers vector with zeros until Simplicial complex dimension + std::vector betti_numbers(cpx_->dimension(), 0); + + for (auto pair : persistent_pairs_) { + // Increment corresponding betti number + betti_numbers[cpx_->dimension(get<0>(pair))] += 1; + } + return betti_numbers; + } + + /** @brief Returns the Betti number passed by parameter. + * @param[in] dimension The Betti number dimension to get. + * @return Betti number + * + */ + const int betti_number(int dimension) const { + int betti_number = 0; + + for (auto pair : persistent_pairs_) { + if (cpx_->dimension(get<0>(pair)) == dimension) + // Increment betti number found + ++betti_number; + } + return betti_number; + } + + /** @brief Returns the persistent pairs. + * @return Persistent pairs + * + */ + std::vector 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), -- cgit v1.2.3 From 0741cec5900f05b57f1aacffcef248fab911d27d Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 26 May 2016 08:17:15 +0000 Subject: Betti number fix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1207 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5dacf1565eee2fe83781b0bdfb05898ffa36fe66 --- .../include/gudhi/Persistent_cohomology.h | 54 +++++++++++++--------- 1 file changed, 33 insertions(+), 21 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 9fc3bf2e..f0cb159e 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -68,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 Column; // contains 1 set_hook // Cell type @@ -80,10 +80,10 @@ class Persistent_cohomology { typedef boost::intrusive::set > 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 > A_ds_type; -// Persistent interval type. The Arith_element field is used for the multi-field framework. - typedef std::tuple Persistent_interval; + // Persistent interval type. The Arith_element field is used for the multi-field framework. + typedef std::tuple Persistent_interval; /** \brief Initializes the Persistent_cohomology class. * @@ -199,17 +199,17 @@ class Persistent_cohomology { if (ds_parent_[key] == key // root of its tree && zero_cocycles_.find(key) == zero_cocycles_.end()) { persistent_pairs_.emplace_back( - cpx_->simplex(key), cpx_->null_simplex()/*, coeff_field_.characteristic()*/); + cpx_->simplex(key), cpx_->null_simplex(), coeff_field_.characteristic()); } } for (auto zero_idx : zero_cocycles_) { persistent_pairs_.emplace_back( - cpx_->simplex(zero_idx.second), cpx_->null_simplex()/*, coeff_field_.characteristic()*/); + 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_*/); + cpx_->simplex(cocycle.first), cpx_->null_simplex(), cocycle.second.characteristics_); } } @@ -250,7 +250,7 @@ class Persistent_cohomology { < cpx_->filtration(cpx_->simplex(idx_coc_v))) { // Kill cocycle [idx_coc_v], which is younger. if (interval_length_policy(cpx_->simplex(idx_coc_v), sigma)) { persistent_pairs_.emplace_back( - cpx_->simplex(idx_coc_v), sigma/*, coeff_field_.characteristic()*/); + cpx_->simplex(idx_coc_v), sigma, coeff_field_.characteristic()); } // Maintain the index of the 0-cocycle alive. if (kv != idx_coc_v) { @@ -265,7 +265,7 @@ class Persistent_cohomology { } else { // Kill cocycle [idx_coc_u], which is younger. if (interval_length_policy(cpx_->simplex(idx_coc_u), sigma)) { persistent_pairs_.emplace_back( - cpx_->simplex(idx_coc_u), sigma/*, coeff_field_.characteristic()*/); + cpx_->simplex(idx_coc_u), sigma, coeff_field_.characteristic()); } // Maintain the index of the 0-cocycle alive. if (ku != idx_coc_u) { @@ -427,7 +427,7 @@ class Persistent_cohomology { if (interval_length_policy(cpx_->simplex(death_key), sigma)) { persistent_pairs_.emplace_back(cpx_->simplex(death_key) // creator , sigma // destructor - /*, charac*/); // fields + , charac); // fields } auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row. @@ -554,6 +554,14 @@ class Persistent_cohomology { Complex_ds * sc_; }; + /* + * Returns true when Filtration_value type accepts infinity values and the given value is equal to infinity. + */ + bool is_infinity(Filtration_value value) const { + bool has_infinity = std::numeric_limits::has_infinity; + return (has_infinity && value == std::numeric_limits::infinity()); + } + public: /** \brief Output the persistence diagram in ostream. * @@ -568,10 +576,9 @@ class Persistent_cohomology { 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::has_infinity; for (auto pair : persistent_pairs_) { // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA) - if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits::infinity()) { + if (is_infinity(cpx_->filtration(get<1>(pair)))) { ostream << /*get<2>(pair) <<*/ " " << cpx_->dimension(get<0>(pair)) << " " << cpx_->filtration(get<0>(pair)) << " inf " << std::endl; } else { @@ -596,13 +603,15 @@ class Persistent_cohomology { /** @brief Returns Betti numbers. * */ - const std::vector betti_numbers() const { + std::vector betti_numbers() const { // Init Betti numbers vector with zeros until Simplicial complex dimension std::vector betti_numbers(cpx_->dimension(), 0); for (auto pair : persistent_pairs_) { - // Increment corresponding betti number - betti_numbers[cpx_->dimension(get<0>(pair))] += 1; + if (is_infinity(cpx_->filtration(get<1>(pair)))) { + // Increment corresponding betti number + betti_numbers[cpx_->dimension(get<0>(pair))] += 1; + } } return betti_numbers; } @@ -612,13 +621,16 @@ class Persistent_cohomology { * @return Betti number * */ - const int betti_number(int dimension) const { + int betti_number(int dimension) const { int betti_number = 0; for (auto pair : persistent_pairs_) { - if (cpx_->dimension(get<0>(pair)) == dimension) - // Increment betti number found - ++betti_number; + if (is_infinity(cpx_->filtration(get<1>(pair)))) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + // Increment betti number found + ++betti_number; + } + } } return betti_number; } @@ -627,7 +639,7 @@ class Persistent_cohomology { * @return Persistent pairs * */ - std::vector get_persistent_pairs() const { + const std::vector& get_persistent_pairs() const { return persistent_pairs_; } -- cgit v1.2.3 From 24fddda9ca1d057c57b1ec9a3d24443c3e1abbf9 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 26 May 2016 21:03:43 +0000 Subject: Add doc, examples and unitary tests git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1213 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: dc7980ef3f1baef17f2e42770da146d6c961e672 --- src/Persistent_cohomology/example/CMakeLists.txt | 5 + .../example/custom_persistence_sort.cpp | 115 ++++++++++ .../example/plain_homology.cpp | 44 +--- .../include/gudhi/Persistent_cohomology.h | 70 ++++-- src/Persistent_cohomology/test/CMakeLists.txt | 8 + .../test/betti_numbers_unit_test.cpp | 234 +++++++++++++++++++++ src/common/doc/main_page.h | 7 + 7 files changed, 422 insertions(+), 61 deletions(-) create mode 100644 src/Persistent_cohomology/example/custom_persistence_sort.cpp create mode 100644 src/Persistent_cohomology/test/betti_numbers_unit_test.cpp (limited to 'src') 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..4dab4560 --- /dev/null +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -0,0 +1,115 @@ +/* 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 . + */ + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include // for std::ofstream + + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Point = Kernel::Point_d; +using Alpha_complex = Gudhi::alphacomplex::Alpha_complex; + +std::vector random_points() { + // Instanciate a random point generator + CGAL::Random rng(0); + + // Generate "points_number" random points in a vector + std::vector points; + + // Generates 1000 random 3D points on a sphere of radius 4.0 + CGAL::Random_points_on_sphere_d 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 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 + 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 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 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 9e5adb5d..5afce9e2 100644 --- a/src/Persistent_cohomology/example/plain_homology.cpp +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -42,24 +42,6 @@ struct MyOptions : Simplex_tree_options_full_featured { }; typedef Simplex_tree ST; - /* - * Compare two intervals by dimension, then by length. - */ - struct cmp_intervals_by_dim_then_length { - explicit cmp_intervals_by_dim_then_length(ST * sc) - : sc_(sc) { - } - template - 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))); - } - ST* sc_; - }; - int main() { ST st; @@ -102,34 +84,10 @@ int main() { pcoh.output_diagram(); - // ******************************************************** - // get_persistence - // ******************************************************** - std::cout << std::endl; + // Print the Betti numbers are b0=2 and b1=1. std::cout << std::endl; - - // First version - std::vector 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; - - // Second version 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; - - // Get persistence - cmp_intervals_by_dim_then_length cmp(&st); - 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 << st.dimension(get<0>(pair)) << " " - << st.filtration(get<0>(pair)) << " " - << st.filtration(get<1>(pair)) << std::endl; - } - - } diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index f0cb159e..5b4c215a 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -88,7 +88,7 @@ class Persistent_cohomology { /** \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), @@ -117,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. @@ -554,14 +554,6 @@ class Persistent_cohomology { Complex_ds * sc_; }; - /* - * Returns true when Filtration_value type accepts infinity values and the given value is equal to infinity. - */ - bool is_infinity(Filtration_value value) const { - bool has_infinity = std::numeric_limits::has_infinity; - return (has_infinity && value == std::numeric_limits::infinity()); - } - public: /** \brief Output the persistence diagram in ostream. * @@ -576,13 +568,14 @@ class Persistent_cohomology { 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::has_infinity; for (auto pair : persistent_pairs_) { // Special case on windows, inf is "1.#INF" (cf. unitary tests and R package TDA) - if (is_infinity(cpx_->filtration(get<1>(pair)))) { - ostream << /*get<2>(pair) <<*/ " " << cpx_->dimension(get<0>(pair)) << " " + if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits::infinity()) { + ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " " << cpx_->filtration(get<0>(pair)) << " inf " << std::endl; } else { - ostream << /*get<2>(pair) <<*/ " " << cpx_->dimension(get<0>(pair)) << " " + ostream << get<2>(pair) << " " << cpx_->dimension(get<0>(pair)) << " " << cpx_->filtration(get<0>(pair)) << " " << cpx_->filtration(get<1>(pair)) << " " << std::endl; } @@ -601,14 +594,15 @@ class Persistent_cohomology { } /** @brief Returns Betti numbers. - * + * @return A vector of persistent Betti numbers. */ std::vector betti_numbers() const { // Init Betti numbers vector with zeros until Simplicial complex dimension std::vector betti_numbers(cpx_->dimension(), 0); for (auto pair : persistent_pairs_) { - if (is_infinity(cpx_->filtration(get<1>(pair)))) { + // Count never ended persistence intervals + if (cpx_->null_simplex() == get<1>(pair)) { // Increment corresponding betti number betti_numbers[cpx_->dimension(get<0>(pair))] += 1; } @@ -616,16 +610,56 @@ class Persistent_cohomology { return betti_numbers; } - /** @brief Returns the Betti number passed by parameter. + /** @brief Returns the Betti number of the dimension passed by parameter. * @param[in] dimension The Betti number dimension to get. - * @return Betti number + * @return Betti number of the given dimension * */ int betti_number(int dimension) const { int betti_number = 0; for (auto pair : persistent_pairs_) { - if (is_infinity(cpx_->filtration(get<1>(pair)))) { + // 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 persistent_betti_numbers(Filtration_value from, Filtration_value to) const { + // Init Betti numbers vector with zeros until Simplicial complex dimension + std::vector 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 persistentBetti number of the dimension passed by parameter. + * @param[in] dimension The Betti number dimension to get. + * @return 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; 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 +#include +#include +#include // std::pair, std::make_pair +#include // float comparison +#include + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "betti_numbers" +#include + +#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; +using Mini_st_persistence = + Gudhi::persistent_cohomology::Persistent_cohomology; + +/* + * Compare two intervals by dimension, then by length. + */ +template +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Simplicial_complex * sc) + : sc_(sc) { } + + template + 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 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 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; + +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 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 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 * \li * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp + * \li + * Persistent_cohomology/custom_persistence_sort.cpp * * \subsection eigen3 Eigen3: * Eigen3 is a C++ template library for linear algebra: @@ -224,6 +226,8 @@ * Persistent_cohomology/alpha_complex_persistence.cpp * \li * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp + * \li + * Persistent_cohomology/custom_persistence_sort.cpp * * \subsection tbb Threading Building Blocks: * Intel® TBB lets you easily write parallel @@ -271,6 +275,8 @@ * Persistent_cohomology/rips_persistence.cpp * \li * Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp + * \li + * Persistent_cohomology/custom_persistence_sort.cpp * * \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 -- cgit v1.2.3 From 270c69e3e96931b6f355af421fee414340c642c2 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 30 May 2016 11:12:25 +0000 Subject: cpplint fixes git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1219 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 81a6208f769c5237d4f81bc97fc1f3d581e733be --- .../example/custom_persistence_sort.cpp | 1 + src/Persistent_cohomology/example/plain_homology.cpp | 9 ++++----- .../include/gudhi/Persistent_cohomology.h | 14 +++++++------- 3 files changed, 12 insertions(+), 12 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/example/custom_persistence_sort.cpp b/src/Persistent_cohomology/example/custom_persistence_sort.cpp index 4dab4560..d591a0d0 100644 --- a/src/Persistent_cohomology/example/custom_persistence_sort.cpp +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -32,6 +32,7 @@ #include #include #include // for std::ofstream +#include // for std::sort using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp index 5afce9e2..2156ffd9 100644 --- a/src/Persistent_cohomology/example/plain_homology.cpp +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -77,13 +77,12 @@ 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 : "; diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 5b4c215a..23beaa96 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -599,7 +599,7 @@ class Persistent_cohomology { std::vector betti_numbers() const { // Init Betti numbers vector with zeros until Simplicial complex dimension std::vector betti_numbers(cpx_->dimension(), 0); - + for (auto pair : persistent_pairs_) { // Count never ended persistence intervals if (cpx_->null_simplex() == get<1>(pair)) { @@ -617,7 +617,7 @@ class Persistent_cohomology { */ 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)) { @@ -629,7 +629,7 @@ class Persistent_cohomology { } 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$. @@ -638,7 +638,7 @@ class Persistent_cohomology { std::vector persistent_betti_numbers(Filtration_value from, Filtration_value to) const { // Init Betti numbers vector with zeros until Simplicial complex dimension std::vector 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) { @@ -656,7 +656,7 @@ class Persistent_cohomology { */ 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) { @@ -668,7 +668,7 @@ class Persistent_cohomology { } return betti_number; } - + /** @brief Returns the persistent pairs. * @return Persistent pairs * @@ -676,7 +676,7 @@ class Persistent_cohomology { const std::vector& get_persistent_pairs() const { return persistent_pairs_; } - + private: /* * Structure representing a cocycle. -- cgit v1.2.3 From ca93a90b895e74cc6bdbeed9273e8cb5de195fc9 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 30 May 2016 12:49:54 +0000 Subject: Missing documentation for parameters. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/get_persistence@1220 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 190cfccdca4d4fd15066ad04885d5dc19428d59c --- src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index 23beaa96..efa9a367 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -649,10 +649,11 @@ class Persistent_cohomology { return betti_numbers; } - /** @brief Returns the persistentBetti number of the dimension passed by parameter. + /** @brief Returns the persistent Betti number of the dimension passed by parameter. * @param[in] dimension The Betti number dimension to get. - * @return Betti number of the given dimension - * + * @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; -- cgit v1.2.3