diff options
Diffstat (limited to 'src/Persistent_cohomology')
28 files changed, 3786 insertions, 0 deletions
diff --git a/src/Persistent_cohomology/benchmark/CMakeLists.txt b/src/Persistent_cohomology/benchmark/CMakeLists.txt new file mode 100644 index 00000000..2bb3b0c7 --- /dev/null +++ b/src/Persistent_cohomology/benchmark/CMakeLists.txt @@ -0,0 +1,12 @@ +project(Persistent_cohomology_benchmark) + +if(GMP_FOUND) + if(GMPXX_FOUND) + add_executable ( performance_rips_persistence EXCLUDE_FROM_ALL performance_rips_persistence.cpp ) + target_link_libraries(performance_rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + file(COPY "${CMAKE_SOURCE_DIR}/data/points/Kl.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + endif(GMPXX_FOUND) +endif(GMP_FOUND) diff --git a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp new file mode 100644 index 00000000..45757002 --- /dev/null +++ b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp @@ -0,0 +1,204 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Persistent_cohomology/Multi_field.h> +#include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> + +#include <chrono> +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +/* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ +template< typename FilteredComplex> +void timing_persistence(FilteredComplex & cpx + , int p); + +/* Compute multi-field persistent homology of the complex cpx with coefficients in + * Z/rZ for all prime number r in [p;q].*/ +template< typename FilteredComplex> +void timing_persistence(FilteredComplex & cpx + , int p + , int q); + +/* Timings for the computation of persistent homology with different + * representations of a Rips complex and different coefficient fields. The + * Rips complex is built on a set of 10000 points sampling a Klein bottle embedded + * in dimension 5. + * We represent complexes with a simplex tree and + * with a Hasse diagram. The Hasse diagram represents explicitly all + * codimension 1 incidence relations in the complex, and hence leads to + * a faster computation of persistence because boundaries are precomputed. + * Hovewer, the simplex tree may be constructed directly from a point cloud and + * is more compact. + * We compute persistent homology with coefficient fields Z/2Z and Z/1223Z. + * We present also timings for the computation of multi-field persistent + * homology in all fields Z/rZ for r prime between 2 and 1223. + */ +int main(int argc, char * argv[]) { + std::chrono::time_point<std::chrono::system_clock> start, end; + int elapsed_sec; + { + + std::string off_file_points = "Kl.off"; + Filtration_value threshold = 0.27; + int dim_max = 3; + int p = 2; + int q = 1223; + + // Extract the points from the file off_file_points + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + start = std::chrono::system_clock::now(); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; + + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + start = std::chrono::system_clock::now(); + + // insert the proximity graph in the simplex tree + // expand the graph until dimension dim_max + rips_complex_from_file.create_complex(st, dim_max); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n"; + std::cout << " - dimension = " << st.dimension() << std::endl; + std::cout << " - number of simplices = " << st.num_simplices() << std::endl; + + // Sort the simplices in the order of the filtration + start = std::chrono::system_clock::now(); + st.initialize_filtration(); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n"; + + // Copy the keys inside the simplices + start = std::chrono::system_clock::now(); + { + int count = 0; + for (auto sh : st.filtration_simplex_range()) + st.assign_key(sh, count++); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n"; + + // Convert the simplex tree into a hasse diagram + start = std::chrono::system_clock::now(); + Gudhi::Hasse_complex<> hcpx(st); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n"; + + + std::cout << "Timings when using a simplex tree: \n"; + timing_persistence(st, p); + timing_persistence(st, q); + timing_persistence(st, p, q); + + std::cout << "Timings when using a Hasse complex: \n"; + timing_persistence(hcpx, p); + timing_persistence(hcpx, q); + timing_persistence(hcpx, p, q); + + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n"; + return 0; +} + +template< typename FilteredComplex> +void +timing_persistence(FilteredComplex & cpx + , int p) { + std::chrono::time_point<std::chrono::system_clock> start, end; + int elapsed_sec; + { + start = std::chrono::system_clock::now(); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; + // initializes the coefficient field for homology + start = std::chrono::system_clock::now(); + pcoh.init_coefficients(p); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; + + start = std::chrono::system_clock::now(); + + pcoh.compute_persistent_cohomology(INFINITY); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n"; + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; +} + +template< typename FilteredComplex> +void +timing_persistence(FilteredComplex & cpx + , int p + , int q) { + std::chrono::time_point<std::chrono::system_clock> start, end; + int elapsed_sec; + { + start = std::chrono::system_clock::now(); + Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; + // initializes the coefficient field for homology + start = std::chrono::system_clock::now(); + pcoh.init_coefficients(p, q); + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; + // compute persistent homology, disgarding persistent features of life shorter than min_persistence + + start = std::chrono::system_clock::now(); + + pcoh.compute_persistent_cohomology(INFINITY); + + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Compute multi-field persistent homology in all coefficient fields Z/pZ " + << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n"; + start = std::chrono::system_clock::now(); + } + end = std::chrono::system_clock::now(); + elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; +} diff --git a/src/Persistent_cohomology/concept/CoefficientField.h b/src/Persistent_cohomology/concept/CoefficientField.h new file mode 100644 index 00000000..916f49e2 --- /dev/null +++ b/src/Persistent_cohomology/concept/CoefficientField.h @@ -0,0 +1,40 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** \brief Concept describing the requirements for a class to represent + * a field of coefficients to compute persistent homology. + */ +struct CoefficientField { + +/** \brief Type of element of the field. + * + * Must be Assignable. */ + typedef unspecified Element; + +/** Default constructible. */ + CoefficientField(); + + void init(Element charac); + void init(Element charac_min, Element charac_max); + +/** Return the characteristic of the field. */ + Element characteristic(); +/** Return the element 1 of the field. */ + Element multiplicative_identity(); +/** Return the element 0 of the field. */ + Element additive_identity(); + +/** Assign: x <- x + y */ + void plus_equal(Element x, Element y); + +/** */ +//... inverse() + + };
\ No newline at end of file diff --git a/src/Persistent_cohomology/concept/FilteredComplex.h b/src/Persistent_cohomology/concept/FilteredComplex.h new file mode 100644 index 00000000..26ac7ac8 --- /dev/null +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -0,0 +1,131 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** \brief The concept FilteredComplex describes the requirements + * for a type to implement a filtered cell complex, from which + * one can compute persistent homology via a model of the concept + * PersistentHomology. + */ +struct FilteredComplex +{ +/** \brief Handle to specify a simplex. */ + typedef unspecified Simplex_handle; +/** \brief Type for the value of the filtration function. + * + * Must be comparable with <. */ + typedef unspecified Filtration_value; + +/** \brief Specifies the nature of the indexing scheme. + * + * is model of IndexingTag. */ + typedef unspecified Indexing_tag; + +/** \brief Returns a Simplex_handle that is different from all simplex handles of the simplices. */ + Simplex_handle null_simplex(); +/** \brief Returns the number of simplices in the complex. + * + * Does not count the empty simplex. */ + size_t num_simplices(); +/** \brief Returns the dimension of a simplex. */ + int dimension(Simplex_handle sh); +/** \brief Returns the filtration value of a simplex. + * + * If sh is null_simplex(), returns the maximal value of the + * filtration function on the complex. */ + Filtration_value filtration(Simplex_handle sh); + +/** \brief Returns the simplex that has index idx in the filtration. + * + * This is only called on valid indices. */ + Simplex_handle simplex ( size_t idx ); +/** \brief Iterator on the simplices belonging to the boundary of a simplex. + * + * <CODE>value_type</CODE> must be 'Simplex_handle'. + */ +typedef unspecified Boundary_simplex_iterator; +/** \brief Range giving access to the simplices in the boundary of a simplex. + * + * .begin() and .end() return type Boundary_simplex_iterator. + */ +typedef unspecified Boundary_simplex_range; + +/** \brief Returns a range giving access to all simplices of the boundary of a simplex, i.e. the set of codimension 1 + * subsimplices of the Simplex. + * + * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation + * induced by \f$ v_0 < \cdots < v_d \f$, the iterator enumerates the + * simplices of the boundary in the order: + * \f$[v_0,\cdots,\widehat{v_i},\cdots,v_d]\f$ for \f$i\f$ from 0 to d + * + * We note that the alternate sum of the simplices given by the iterator + * gives the chains corresponding to the boundary of the simplex.*/ +Boundary_simplex_range boundary_simplex_range(Simplex_handle sh); + +/** \brief Iterator over all simplices of the complex in the order of the indexing scheme. + * + * 'value_type' must be 'Simplex_handle'. + */ +typedef unspecified Filtration_simplex_iterator; +/** \brief Range over the simplices of the complex in the order of the filtration. + * + * .begin() and .end() return type Filtration_simplex_iterator.*/ +typedef unspecified Filtration_simplex_range; +/** \brief Returns a range over the simplices of the complex in the order of the filtration. + * + * .begin() and .end() return type Filtration_simplex_iterator.*/ +Filtration_simplex_range filtration_simplex_range(); + +/** \name Map interface + * Conceptually a `std::unordered_map<Simplex_handle,std::size_t>`. + * @{ */ +/** \brief Data stored for each simplex. + * + * Must be an integer type. */ + typedef unspecified Simplex_key; +/** \brief Returns a constant dummy number that is either negative, + * or at least as large as `num_simplices()`. Suggested value: -1. */ + Simplex_key null_key (); +/** \brief Returns the number stored for a simplex by `assign_key`. + * + * This is never called on null_simplex(). */ + Simplex_key key ( Simplex_handle sh ); +/** \brief Store a number for a simplex, which can later be retrieved with `key(sh)`. + * + * This is never called on null_simplex(). */ + void assign_key(Simplex_handle sh, Simplex_key n); +/** @} */ + + +/* \brief Iterator over the simplices of the complex, + * in an arbitrary order. + * + * 'value_type' must be 'Simplex_handle'.*/ +//typedef unspecified Complex_simplex_iterator; +//typedef unspecified Complex_simplex_range; + +/* +* Returns a range over all the simplices of a +* complex. +*/ +//Complex_simplex_range complex_simplex_range(); + +/*************************************************/ /** +* @details Compares the filtration values of simplices s and t +* +* @return -1 if s comes before t in the filtration, +1 +* if it comes later, and 0 if they come at the same time +* +* @note OPTIONAL +* @todo use an enum? Just a bool? +*/ +//int is_before_in_filtration(Simplex_handle s, Simplex_handle t); +/*************************************************/ + +}; diff --git a/src/Persistent_cohomology/concept/PersistentHomology.h b/src/Persistent_cohomology/concept/PersistentHomology.h new file mode 100644 index 00000000..373832af --- /dev/null +++ b/src/Persistent_cohomology/concept/PersistentHomology.h @@ -0,0 +1,29 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +/** \brief Concept describing the requirements for a class to compute + * persistent homology. */ +struct PersistentHomology { + +/** \brief Type of filtered cell complex on which persistent homology + * is computed. + * + * Must be a model of concept FilteredComplex. + */ + typedef unspecified Filtered_complex; + +/** \brief Type of coefficients to be used for computing persistent + * homology. + * + * Must be a model of concept CoefficientField. + */ + typedef unspecified Coefficient_field; + + }; diff --git a/src/Persistent_cohomology/doc/3DTorus_poch.png b/src/Persistent_cohomology/doc/3DTorus_poch.png Binary files differnew file mode 100644 index 00000000..1c9d8328 --- /dev/null +++ b/src/Persistent_cohomology/doc/3DTorus_poch.png diff --git a/src/Persistent_cohomology/doc/COPYRIGHT b/src/Persistent_cohomology/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Persistent_cohomology/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + +Author(s): Vincent Rouvreau + +Copyright (C) 2015 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT.
\ No newline at end of file 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..46b784d8 --- /dev/null +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -0,0 +1,261 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#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. + + \section persistencetopolocalspaces Topological Spaces + 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. + + \section persistencehomology Homology + 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. + + \section persistenceindexingscheme Indexing Scheme + "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 pcohexamples Examples + +We provide several example files: run these examples with -h for details on their use, and read the README file. + +\li <a href="_rips_complex_2rips_persistence_8cpp-example.html"> +Rips_complex/rips_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its persistence +diagram. +\code $> ./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3 \endcode +\code The complex contains 177838 simplices + and has dimension 3 +3 0 0 inf +3 1 0.0983494 inf +3 1 0.104347 inf +3 2 0.138335 inf \endcode + +More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. + +\li <a href="_persistent_cohomology_2rips_multifield_persistence_8cpp-example.html"> +Persistent_cohomology/rips_multifield_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its +persistence diagram with a family of field coefficients. + +\li <a href="_rips_complex_2rips_distance_matrix_persistence_8cpp-example.html"> +Rips_complex/rips_distance_matrix_persistence.cpp</a> computes the Rips complex of a distance matrix and +outputs its persistence diagram. + +The file should contain square or lower triangular distance matrix with semicolons as separators. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. + +More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. + +\li <a href="_rips_complex_2rips_correlation_matrix_persistence_8cpp-example.html"> +Rips_complex/rips_correlation_matrix_persistence.cpp</a> +computes the Rips complex of a correlation matrix and outputs its persistence diagram. + +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. The input is to be given either as a square or a lower +triangular matrix. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + +More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. + +\li <a href="_alpha_complex_2alpha_complex_3d_persistence_8cpp-example.html"> +Alpha_complex/alpha_complex_3d_persistence.cpp</a> computes the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. +\code $> ./alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45 \endcode +\code Simplex_tree dim: 3 +2 0 0 inf +2 1 0.0682162 1.0001 +2 1 0.0934117 1.00003 +2 2 0.56444 1.03938 \endcode + +More details on the <a href="../../alphacomplex/">Alpha complex utilities</a> dedicated page. + +CGAL can be forced to compute the exact values, it is slower, but it is necessary when points are on a grid +for instance (the fast version `--fast` would give incorrect values). +\code $> ./alpha_complex_3d_persistence ../../data/points/sphere3D_pts_on_grid.off --exact -p 2 -m 0.1 \endcode +\code Simplex_tree dim: 3 +2 0 0 inf +2 2 0.0002 0.2028 \endcode + +It can also compute the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the weighted alpha complex on points sampling from an OFF file +and a weights file. +\code $> ./alpha_complex_3d_persistence ../../data/points/tore3D_300.off +--weight-file ../../data/points/tore3D_300.weights -p 2 -m 0.45 \endcode +\code Simplex_tree dim: 3 +2 0 -1 inf +2 1 -0.931784 0.000103311 +2 1 -0.906588 2.60165e-05 +2 2 -0.43556 0.0393798 \endcode + +One can also compute the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the periodic alpha complex on points sampling from an OFF file. +The second parameter is a \ref FileFormatsIsoCuboid file with coordinates of the periodic cuboid. +Note that the lengths of the sides of the periodic cuboid have to be the same. +\code $> ./alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off +--cuboid-file ../../data/points/iso_cuboid_3_in_0_1.txt -p 3 -m 1.0 \endcode +\code Simplex_tree dim: 3 +3 0 0 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 2 0.005 inf +3 2 0.005 inf +3 2 0.005 inf +3 3 0.0075 inf \endcode + +In order to compute the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the periodic alpha complex on weighted points from an OFF file. The +additional parameters of this program are:<br> +(a) The file with the weights of points. The file consist of a sequence of numbers (as many as points). +Note that the weight of each single point have to be bounded by 1/64 times the square of the cuboid edge length.<br> +(b) A \ref FileFormatsIsoCuboid file with coordinates of the periodic cuboid. +Note that the lengths of the sides of the periodic cuboid have to be the same.<br> +\code $> ./alpha_complex_3d_persistence ../../data/points/shifted_sphere.off +--weight-file ../../data/points/shifted_sphere.weights +--cuboid-file ../../data/points/iso_cuboid_3_in_0_10.txt -p 3 -m 1.0 \endcode +\code Simplex_tree dim: 3 +3 0 -0.0001 inf +3 1 16.0264 inf +3 1 16.0273 inf +3 1 16.0303 inf +3 2 36.8635 inf +3 2 36.8704 inf +3 2 36.8838 inf +3 3 58.6783 inf \endcode + +\li <a href="_alpha_complex_2alpha_complex_persistence_8cpp-example.html"> +Alpha_complex/alpha_complex_persistence.cpp</a> computes the persistent homology with +\f$\mathbb{Z}/p\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. +\code $> ./alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off \endcode +\code Alpha complex is of dimension 3 - 9273 simplices - 300 vertices. +Simplex_tree dim: 3 +2 0 0 inf +2 1 0.0682162 1.0001 +2 1 0.0934117 1.00003 +2 2 0.56444 1.03938 \endcode + +More details on the <a href="../../alphacomplex/">Alpha complex utilities</a> dedicated page. + +\li <a href="_persistent_cohomology_2plain_homology_8cpp-example.html"> +Persistent_cohomology/plain_homology.cpp</a> computes the plain homology of a simple simplicial complex without +filtration values. + + */ + +} // 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 new file mode 100644 index 00000000..94ec13c5 --- /dev/null +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -0,0 +1,68 @@ +project(Persistent_cohomology_examples) + +add_executable(plain_homology plain_homology.cpp) + +add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp) + +add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp) +target_link_libraries(rips_persistence_step_by_step ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(rips_persistence_via_boundary_matrix rips_persistence_via_boundary_matrix.cpp) +target_link_libraries(rips_persistence_via_boundary_matrix ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(persistence_from_file persistence_from_file.cpp) +target_link_libraries(persistence_from_file ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(plain_homology ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_step_by_step ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_file ${TBB_LIBRARIES}) +endif() + +add_test(NAME Persistent_cohomology_example_plain_homology COMMAND $<TARGET_FILE:plain_homology>) +add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $<TARGET_FILE:persistence_from_simple_simplex_tree> + "1" "0") +add_test(NAME Persistent_cohomology_example_from_rips_step_by_step_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence_step_by_step> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +add_test(NAME Persistent_cohomology_example_via_boundary_matrix COMMAND $<TARGET_FILE:rips_persistence_via_boundary_matrix> + "${CMAKE_SOURCE_DIR}/data/points/Kl.off" "-r" "0.16" "-d" "3" "-p" "3" "-m" "100") +add_test(NAME Persistent_cohomology_example_from_file_3_2_0 COMMAND $<TARGET_FILE:persistence_from_file> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc" "-p" "2" "-m" "0") +add_test(NAME Persistent_cohomology_example_from_file_3_3_100 COMMAND $<TARGET_FILE:persistence_from_file> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc" "-p" "3" "-m" "100") + +install(TARGETS plain_homology DESTINATION bin) +install(TARGETS persistence_from_simple_simplex_tree DESTINATION bin) +install(TARGETS rips_persistence_step_by_step DESTINATION bin) +install(TARGETS rips_persistence_via_boundary_matrix DESTINATION bin) +install(TARGETS persistence_from_file DESTINATION bin) + +if(GMP_FOUND) + if(GMPXX_FOUND) + add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) + target_link_libraries(rips_multifield_persistence + ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistent_cohomology_example_multifield_2_71 COMMAND $<TARGET_FILE:rips_multifield_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "2" "-q" "71") + install(TARGETS rips_multifield_persistence DESTINATION bin) + endif(GMPXX_FOUND) +endif(GMP_FOUND) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + + add_executable(custom_persistence_sort custom_persistence_sort.cpp) + target_link_libraries(custom_persistence_sort ${CGAL_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistent_cohomology_example_custom_persistence_sort COMMAND $<TARGET_FILE:custom_persistence_sort>) + + install(TARGETS custom_persistence_sort DESTINATION bin) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Persistent_cohomology/example/README b/src/Persistent_cohomology/example/README new file mode 100644 index 00000000..f39d9584 --- /dev/null +++ b/src/Persistent_cohomology/example/README @@ -0,0 +1,67 @@ +To build the examples, run in a Terminal: + +cd /path-to-examples/ +cmake . +make + +*********************************************************************************************************************** +Example of use of RIPS: + +Computation of the persistent homology with Z/2Z and Z/3Z coefficients simultaneously of the Rips complex +on points sampling a 3D torus: + +./rips_multifield_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.12 -d 3 -p 2 -q 3 + +output: +6 0 0 inf +6 1 0.0983494 inf +6 1 0.104347 inf +6 2 0.138335 inf +6 0 0 0.122545 +6 0 0 0.121171 +6 0 0 0.120964 +6 0 0 0.12057 +6 0 0 0.12047 +6 0 0 0.120414 + +Every line is of this format: p1*...*pr dim b d +where + p1*...*pr is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients. + dim is the dimension of the homological feature, + b and d are respectively the birth and death of the feature and + +and the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers): + + ./rips_multifield_persistence ../../data/points/Kl.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71 + +output: +557940830126698960967415390 0 0 inf +557940830126698960967415390 1 0.0983494 inf +557940830126698960967415390 1 0.104347 inf +557940830126698960967415390 2 0.138335 inf +557940830126698960967415390 0 0 0.122545 +557940830126698960967415390 0 0 0.121171 +557940830126698960967415390 0 0 0.120964 +557940830126698960967415390 0 0 0.12057 +557940830126698960967415390 0 0 0.12047 +557940830126698960967415390 0 0 0.120414 + +*********************************************************************************************************************** +Example of use of PLAIN HOMOLOGY: + +This example computes the plain homology of the following simplicial complex without filtration values: + /* Complex to build. */ + /* 1 3 */ + /* o---o */ + /* /X\ / */ + /* o---o o */ + /* 2 0 4 */ + +./plain_homology + +output: +2 0 0 inf +2 0 0 inf +2 1 0 inf + +Here we retrieve the 2 entities {0,1,2,3} and {4} (Betti numbers[0] = 2) and the hole in {0,1,3} (Betti numbers[1] = 1) 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..be74cf50 --- /dev/null +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -0,0 +1,125 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#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> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.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::alpha_complex::Alpha_complex<Kernel>; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, + Gudhi::persistent_cohomology::Field_Zp >; + +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(Simplex_tree * 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))); + } + Simplex_tree* sc_; +}; + +int main(int argc, char **argv) { + std::vector<Point> points = random_points(); + + std::cout << "Points size=" << points.size() << std::endl; + // Alpha complex persistence computation from generated points + Alpha_complex alpha_complex_from_points(points); + std::cout << "alpha_complex_from_points" << std::endl; + + Simplex_tree simplex; + std::cout << "simplex" << std::endl; + if (alpha_complex_from_points.create_complex(simplex, 0.6)) { + std::cout << "simplex" << std::endl; + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Simplicial complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + // Sort the simplices in the order of the filtration + simplex.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; + + Persistent_cohomology pcoh(simplex); + + // 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(&simplex); + 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 << simplex.dimension(get<0>(pair)) << " " + << simplex.filtration(get<0>(pair)) << " " + << simplex.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 < simplex.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/persistence_from_file.cpp b/src/Persistent_cohomology/example/persistence_from_file.cpp new file mode 100644 index 00000000..d169cc63 --- /dev/null +++ b/src/Persistent_cohomology/example/persistence_from_file.cpp @@ -0,0 +1,130 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/reader_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <boost/program_options.hpp> + +#include <string> + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; + +typedef double Filtration_value; + +void program_options(int argc, char * argv[] + , std::string & simplex_tree_file + , std::string & output_file + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string simplex_tree_file; + std::string output_file; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, simplex_tree_file, output_file, p, min_persistence); + + std::cout << "Simplex_tree from file=" << simplex_tree_file.c_str() << " - output_file=" << output_file.c_str() + << std::endl; + std::cout << " - p=" << p << " - min_persistence=" << min_persistence << std::endl; + + // Read the list of simplices from a file. + Simplex_tree<> simplex_tree; + + std::ifstream simplex_tree_stream(simplex_tree_file); + simplex_tree_stream >> simplex_tree; + + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices" << std::endl; + std::cout << " - dimension " << simplex_tree.dimension() << std::endl; + + /* + std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for( auto f_simplex : simplex_tree.filtration_simplex_range() ) + { std::cout << " " << "[" << simplex_tree.filtration(f_simplex) << "] "; + for( auto vertex : simplex_tree.simplex_vertex_range(f_simplex) ) + { std::cout << vertex << " "; } + std::cout << std::endl; + }*/ + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in output_file + if (output_file.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(output_file); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & simplex_tree_file + , std::string & output_file + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&simplex_tree_file), + "Name of file containing a simplex set. Format is one simplex per line (cf. reader_utils.h - read_simplex): Dim1 X11 X12 ... X1d Fil1 "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&output_file)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp new file mode 100644 index 00000000..3c91662f --- /dev/null +++ b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp @@ -0,0 +1,163 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <ctime> +#include <utility> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using typeVectorVertex = std::vector< Simplex_tree::Vertex_handle >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; + exit(-1); +} + +int main(int argc, char * const argv[]) { + // program args management + if (argc != 3) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = 0; + int returnedScanValue = sscanf(argv[1], "%d", &coeff_field_characteristic); + if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) { + std::cerr << "Error: " << argv[1] << " is not correct\n"; + usage(argv[0]); + } + + Filtration_value min_persistence = 0.0; + returnedScanValue = sscanf(argv[2], "%lf", &min_persistence); + if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { + std::cerr << "Error: " << argv[2] << " is not correct\n"; + usage(argv[0]); + } + + // TEST OF INSERTION + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF INSERTION" << std::endl; + Simplex_tree st; + + // ++ FIRST + std::cout << " - INSERT (0,1,2)" << std::endl; + typeVectorVertex SimplexVector = {0, 1, 2}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ SECOND + std::cout << " - INSERT 3" << std::endl; + SimplexVector = {3}; + st.insert_simplex_and_subfaces(SimplexVector, 0.1); + + // ++ THIRD + std::cout << " - INSERT (0,3)" << std::endl; + SimplexVector = {0, 3}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ FOURTH + std::cout << " - INSERT (0,1) (already inserted)" << std::endl; + SimplexVector = {0, 1}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ FIFTH + std::cout << " - INSERT (3,4,5)" << std::endl; + SimplexVector = {3, 4, 5}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ SIXTH + std::cout << " - INSERT (0,1,6,7)" << std::endl; + SimplexVector = {0, 1, 6, 7}; + st.insert_simplex_and_subfaces(SimplexVector, 0.4); + + // ++ SEVENTH + std::cout << " - INSERT (4,5,8,9)" << std::endl; + SimplexVector = {4, 5, 8, 9}; + st.insert_simplex_and_subfaces(SimplexVector, 0.4); + + // ++ EIGHTH + std::cout << " - INSERT (9,10,11)" << std::endl; + SimplexVector = {9, 10, 11}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ NINETH + std::cout << " - INSERT (2,10,12)" << std::endl; + SimplexVector = {2, 10, 12}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ TENTH + std::cout << " - INSERT (11,6)" << std::endl; + SimplexVector = {6, 11}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ ELEVENTH + std::cout << " - INSERT (13,14,15)" << std::endl; + SimplexVector = {13, 14, 15}; + st.insert_simplex_and_subfaces(SimplexVector, 0.25); + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ 4 2 */ + /* o---o---o---o o */ + /* 2 0 3\X/8\ 10 /X\ */ + /* o---o---o---o */ + /* 5 9\X/ 12 */ + /* o---o */ + /* 11 6 */ + /* In other words: */ + /* A facet [2,1,0] */ + /* An edge [0,3] */ + /* A facet [3,4,5] */ + /* A cell [0,1,6,7] */ + /* A cell [4,5,8,9] */ + /* A facet [9,10,11] */ + /* An edge [11,6] */ + /* An edge [10,12,2] */ + + + std::cout << "The complex contains " << st.num_simplices() << " simplices - " << st.num_vertices() << " vertices " + << std::endl; + std::cout << " - dimension " << st.dimension() << std::endl; + std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" + << std::endl; + std::cout << "**************************************************************" << std::endl; + std::cout << "strict graph G { " << std::endl; + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << " " << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << static_cast<int>(vertex) << " -- "; + } + std::cout << ";" << std::endl; + } + + std::cout << "}" << std::endl; + std::cout << "**************************************************************" << std::endl; + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + pcoh.output_diagram(); + return 0; +} diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp new file mode 100644 index 00000000..84333e46 --- /dev/null +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -0,0 +1,91 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <vector> +#include <cstdint> // for std::uint8_t + +/* We could perfectly well use the default Simplex_tree<> (which uses + * Simplex_tree_options_full_featured), the following simply demonstrates + * how to save on storage by not storing a filtration value. */ + +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; + // Maximum number of simplices to compute persistence is 2^8 - 1 = 255. One is reserved for null_key + typedef std::uint8_t Simplex_key; +}; + +using ST = Gudhi::Simplex_tree<MyOptions>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<ST, Field_Zp>; + +int main() { + ST st; + + /* Complex to build. */ + /* 1 3 5 */ + /* o---o---o */ + /* / \ / */ + /* o---o o */ + /* 2 0 4 */ + + const short edge01[] = {0, 1}; + const short edge02[] = {0, 2}; + const short edge12[] = {1, 2}; + const short edge03[] = {0, 3}; + const short edge13[] = {1, 3}; + const short edge35[] = {3, 5}; + const short vertex4[] = {4}; + st.insert_simplex_and_subfaces(edge01); + st.insert_simplex_and_subfaces(edge02); + st.insert_simplex_and_subfaces(edge12); + st.insert_simplex_and_subfaces(edge03); + st.insert_simplex(edge13); + st.insert_simplex_and_subfaces(edge35); + st.insert_simplex(vertex4); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Class for homology computation + // By default, since the complex has dimension 1, only 0-dimensional homology would be computed. + // Here we also want persistent homology to be computed for the maximal dimension in the complex (persistence_dim_max = true) + Persistent_cohomology pcoh(st, true); + + // Initialize the coefficient field Z/2Z for homology + pcoh.init_coefficients(2); + + // 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 + // 2 1 0 inf + // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=2. + pcoh.output_diagram(); + + // Print the Betti numbers are b0=2 and b1=2. + std::cout << std::endl; + std::cout << "The Betti numbers are : "; + for (int i = 0; i < 3; i++) + std::cout << "b" << i << " = " << pcoh.betti_number(i) << " ; "; + std::cout << std::endl; +} diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp new file mode 100644 index 00000000..9eb5ccfc --- /dev/null +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -0,0 +1,142 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Persistent_cohomology/Multi_field.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Multi_field >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & min_p + , int & max_p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int min_p; + int max_p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(min_p, max_p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & min_p + , int & max_p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options"); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("min-field-charac,p", po::value<int>(&min_p)->default_value(2), + "Minimal characteristic p of the coefficient field Z/pZ.") + ("max-field-charac,q", po::value<int>(&max_p)->default_value(1223), + "Minimial characteristic q of the coefficient field Z/pZ.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with various coefficient fields \n"; + std::cout << "of a Rips complex defined on a set of input points. The coefficient \n"; + std::cout << "fields are all the Z/rZ for a prime number r contained in the \n"; + std::cout << "specified range [p,q]\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p1*...*pr dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p1*...*pr is the product of prime numbers pi such that the homology \n"; + std::cout << "feature exists in homology with Z/piZ coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp new file mode 100644 index 00000000..02db05ec --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -0,0 +1,154 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <map> + +// ---------------------------------------------------------------------------- +// rips_persistence_step_by_step is an example of each step that is required to +// build a Rips over a Simplex_tree. Please refer to rips_persistence to see +// how to do the same thing with the Rips_complex wrapper for less detailed +// steps. +// ---------------------------------------------------------------------------- + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + // insert the proximity graph in the simplex tree + st.insert_graph(prox_graph); + // expand the graph until dimension dim_max + st.expansion(dim_max); + + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << st.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp new file mode 100644 index 00000000..37fa5e93 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -0,0 +1,160 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria, Marc Glisse + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> + +#include <boost/program_options.hpp> + +#ifdef GUDHI_USE_TBB +#include <tbb/task_scheduler_init.h> +#endif + +#include <string> +#include <vector> + +//////////////////////////////////////////////////////////////// +// // +// WARNING: persistence computation itself is not parallel, // +// and this uses more memory than rips_persistence. // +// // +//////////////////////////////////////////////////////////////// + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree& st = *new Simplex_tree; + rips_complex_from_file.create_complex(st, dim_max); + + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << st.dimension() << " \n"; + +#ifdef GUDHI_USE_TBB + // Unnecessary, but clarifies which operations are parallel. + tbb::task_scheduler_init ts; +#endif + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + int count = 0; + for (auto sh : st.filtration_simplex_range()) + st.assign_key(sh, count++); + + // Convert to a more convenient representation. + Gudhi::Hasse_complex<> hcpx(st); + +#ifdef GUDHI_USE_TBB + ts.terminate(); +#endif + + // Free some space. + delete &st; + + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Hasse_complex<>, Field_Zp > pcoh(hcpx); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h new file mode 100644 index 00000000..944b6d35 --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -0,0 +1,756 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_H_ +#define PERSISTENT_COHOMOLOGY_H_ + +#include <gudhi/Persistent_cohomology/Persistent_cohomology_column.h> +#include <gudhi/Persistent_cohomology/Field_Zp.h> +#include <gudhi/Simple_object_pool.h> + +#include <boost/intrusive/set.hpp> +#include <boost/pending/disjoint_sets.hpp> +#include <boost/intrusive/list.hpp> + +#include <map> +#include <utility> +#include <list> +#include <vector> +#include <set> +#include <fstream> // std::ofstream +#include <limits> // for numeric_limits<> +#include <tuple> +#include <algorithm> +#include <string> +#include <stdexcept> // for std::out_of_range + +namespace Gudhi { + +namespace persistent_cohomology { + +/** \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) + * \cite boissonnat:hal-00922572 . + * + * \implements PersistentHomology + * + */ +// TODO(CM): Memory allocation policy: classic, use a mempool, etc. +template<class FilteredComplex, class CoefficientField> +class Persistent_cohomology { + public: + // Data attached to each simplex to interface with a Property Map. + + /** \brief Data stored for each simplex. */ + typedef typename FilteredComplex::Simplex_key Simplex_key; + /** \brief Handle to specify a simplex. */ + typedef typename FilteredComplex::Simplex_handle Simplex_handle; + /** \brief Type for the value of the filtration function. */ + typedef typename FilteredComplex::Filtration_value Filtration_value; + /** \brief Type of element of the field. */ + typedef typename CoefficientField::Element Arith_element; + /** \brief Type for birth and death FilteredComplex::Simplex_handle. + * The Arith_element field is used for the multi-field framework. */ + typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval; + + private: + // Compressed Annotation Matrix types: + // Column type + typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column; // contains 1 set_hook + // Cell type + typedef typename Column::Cell Cell; // contains 2 list_hooks + // Remark: constant_time_size must be false because base_hook_cam_h has auto_unlink link_mode + typedef boost::intrusive::list<Cell, + boost::intrusive::constant_time_size<false>, + boost::intrusive::base_hook<base_hook_cam_h> > Hcell; + + typedef boost::intrusive::set<Column, + boost::intrusive::constant_time_size<false> > Cam; + // Sparse column type for the annotation of the boundary of an element. + typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type; + + public: + /** \brief Initializes the Persistent_cohomology class. + * + * @param[in] cpx Complex for which the persistent homology is computed. + * 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. + * + * @exception std::out_of_range In case the number of simplices is more than Simplex_key type numeric limit. + */ + explicit Persistent_cohomology(FilteredComplex& cpx, bool persistence_dim_max = false) + : cpx_(&cpx), + dim_max_(cpx.dimension()), // upper bound on the dimension of the simplices + coeff_field_(), // initialize the field coefficient structure. + num_simplices_(cpx_->num_simplices()), // num_simplices save to avoid to call thrice the function + ds_rank_(num_simplices_), // union-find + ds_parent_(num_simplices_), // union-find + ds_repr_(num_simplices_, NULL), // union-find -> annotation vectors + dsets_(&ds_rank_[0], &ds_parent_[0]), // union-find + cam_(), // collection of annotation vectors + zero_cocycles_(), // union-find -> Simplex_key of creator for 0-homology + transverse_idx_(), // key -> row + persistent_pairs_(), + interval_length_policy(&cpx, 0), + column_pool_(), // memory pools for the CAM + cell_pool_() { + if (cpx_->num_simplices() > std::numeric_limits<Simplex_key>::max()) { + // num_simplices must be strictly lower than the limit, because a value is reserved for null_key. + throw std::out_of_range("The number of simplices is more than Simplex_key type numeric limit."); + } + Simplex_key idx_fil = 0; + for (auto sh : cpx_->filtration_simplex_range()) { + cpx_->assign_key(sh, idx_fil); + ++idx_fil; + dsets_.make_set(cpx_->key(sh)); + } + if (persistence_dim_max) { + ++dim_max_; + } + } + + ~Persistent_cohomology() { + // Clean the transversal lists + for (auto & transverse_ref : transverse_idx_) { + // Destruct all the cells + transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();}); + delete transverse_ref.second.row_; + } + } + + private: + struct length_interval { + length_interval(FilteredComplex * cpx, Filtration_value min_length) + : cpx_(cpx), + min_length_(min_length) { + } + + bool operator()(Simplex_handle sh1, Simplex_handle sh2) { + return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_; + } + + void set_length(Filtration_value new_length) { + min_length_ = new_length; + } + + FilteredComplex * cpx_; + Filtration_value min_length_; + }; + + public: + /** \brief Initializes the coefficient field.*/ + void init_coefficients(int charac) { + coeff_field_.init(charac); + } + /** \brief Initializes the coefficient field for multi-field persistent homology.*/ + void init_coefficients(int charac_min, int charac_max) { + coeff_field_.init(charac_min, charac_max); + } + + /** \brief Compute the persistent homology of the filtered simplicial + * complex. + * + * @param[in] min_interval_length the computation discards all intervals of length + * less or equal than min_interval_length + * + * Assumes that the filtration provided by the simplicial complex is + * valid. Undefined behavior otherwise. */ + void compute_persistent_cohomology(Filtration_value min_interval_length = 0) { + interval_length_policy.set_length(min_interval_length); + // Compute all finite intervals + for (auto sh : cpx_->filtration_simplex_range()) { + int dim_simplex = cpx_->dimension(sh); + switch (dim_simplex) { + case 0: + break; + case 1: + update_cohomology_groups_edge(sh); + break; + default: + update_cohomology_groups(sh, dim_simplex); + break; + } + } + // Compute infinite intervals of dimension 0 + Simplex_key key; + for (auto v_sh : cpx_->skeleton_simplex_range(0)) { // for all 0-dimensional simplices + key = cpx_->key(v_sh); + + 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()); + } + } + for (auto zero_idx : zero_cocycles_) { + persistent_pairs_.emplace_back( + 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_); + } + } + + private: + /** \brief Update the cohomology groups under the insertion of an edge. + * + * The 0-homology is maintained with a simple Union-Find data structure, which + * explains the existance of a specific function of edge insertions. */ + void update_cohomology_groups_edge(Simplex_handle sigma) { + Simplex_handle u, v; + boost::tie(u, v) = cpx_->endpoints(sigma); + + Simplex_key ku = dsets_.find_set(cpx_->key(u)); + Simplex_key kv = dsets_.find_set(cpx_->key(v)); + + if (ku != kv) { // Destroy a connected component + dsets_.link(ku, kv); + // Keys of the simplices which created the connected components containing + // respectively u and v. + Simplex_key idx_coc_u, idx_coc_v; + auto map_it_u = zero_cocycles_.find(ku); + // If the index of the cocycle representing the class is already ku. + if (map_it_u == zero_cocycles_.end()) { + idx_coc_u = ku; + } else { + idx_coc_u = map_it_u->second; + } + + auto map_it_v = zero_cocycles_.find(kv); + // If the index of the cocycle representing the class is already kv. + if (map_it_v == zero_cocycles_.end()) { + idx_coc_v = kv; + } else { + idx_coc_v = map_it_v->second; + } + + if (cpx_->filtration(cpx_->simplex(idx_coc_u)) + < 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()); + } + // Maintain the index of the 0-cocycle alive. + if (kv != idx_coc_v) { + zero_cocycles_.erase(map_it_v); + } + if (kv == dsets_.find_set(kv)) { + if (ku != idx_coc_u) { + zero_cocycles_.erase(map_it_u); + } + zero_cocycles_[kv] = idx_coc_u; + } + } 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()); + } + // Maintain the index of the 0-cocycle alive. + if (ku != idx_coc_u) { + zero_cocycles_.erase(map_it_u); + } + if (ku == dsets_.find_set(ku)) { + if (kv != idx_coc_v) { + zero_cocycles_.erase(map_it_v); + } + zero_cocycles_[ku] = idx_coc_v; + } + } + cpx_->assign_key(sigma, cpx_->null_key()); + } else if (dim_max_ > 1) { // If ku == kv, same connected component: create a 1-cocycle class. + create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic()); + } + } + + /* + * Compute the annotation of the boundary of a simplex. + */ + void annotation_of_the_boundary( + std::map<Simplex_key, Arith_element> & map_a_ds, Simplex_handle sigma, + int dim_sigma) { + // traverses the boundary of sigma, keeps track of the annotation vectors, + // with multiplicity. We used to sum the coefficients directly in + // annotations_in_boundary by using a map, we now do it later. + typedef std::pair<Column *, int> annotation_t; +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL + thread_local +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::vector<annotation_t> annotations_in_boundary; + annotations_in_boundary.clear(); + int sign = 1 - 2 * (dim_sigma % 2); // \in {-1,1} provides the sign in the + // alternate sum in the boundary. + Simplex_key key; + Column * curr_col; + + for (auto sh : cpx_->boundary_simplex_range(sigma)) { + key = cpx_->key(sh); + if (key != cpx_->null_key()) { // A simplex with null_key is a killer, and have null annotation + // 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); + } + } + 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; }); + + // Sum the annotations with multiplicity, using a map<key,coeff> + // to represent a sparse vector. + std::pair<typename std::map<Simplex_key, Arith_element>::iterator, bool> result_insert_a_ds; + + for (auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); /**/) { + 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; + } + // 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, + for (auto cell_ref : col->col_) { // insert every cell in map_a_ds with multiplicity + Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult); // coefficient * multiplicity + + if (w_y != coeff_field_.additive_identity()) { // if != 0 + result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y)); + if (!(result_insert_a_ds.second)) { // if cell_ref.key_ already a Key in map_a_ds + result_insert_a_ds.first->second = coeff_field_.plus_equal(result_insert_a_ds.first->second, w_y); + if (result_insert_a_ds.first->second == coeff_field_.additive_identity()) { + map_a_ds.erase(result_insert_a_ds.first); + } + } + } + } + } + } + } + + /* + * Update the cohomology groups under the insertion of a simplex. + */ + void update_cohomology_groups(Simplex_handle sigma, int dim_sigma) { +// Compute the annotation of the boundary of sigma: + std::map<Simplex_key, Arith_element> map_a_ds; + annotation_of_the_boundary(map_a_ds, sigma, dim_sigma); +// Update the cohomology groups: + if (map_a_ds.empty()) { // sigma is a creator in all fields represented in coeff_field_ + if (dim_sigma < dim_max_) { + create_cocycle(sigma, coeff_field_.multiplicative_identity(), + coeff_field_.characteristic()); + } + } else { // sigma is a destructor in at least a field in coeff_field_ + // Convert map_a_ds to a vector + A_ds_type a_ds; // admits reverse iterators + for (auto map_a_ds_ref : map_a_ds) { + a_ds.push_back( + std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first, + map_a_ds_ref.second)); + } + + Arith_element inv_x, charac; + Arith_element prod = coeff_field_.characteristic(); // Product of characteristic of the fields + for (auto a_ds_rit = a_ds.rbegin(); + (a_ds_rit != a_ds.rend()) + && (prod != coeff_field_.multiplicative_identity()); ++a_ds_rit) { + std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod); + + if (inv_x != coeff_field_.additive_identity()) { + destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac); + prod /= charac; + } + } + if (prod != coeff_field_.multiplicative_identity() + && dim_sigma < dim_max_) { + create_cocycle(sigma, coeff_field_.multiplicative_identity(prod), prod); + } + } + } + + /* \brief Create a new cocycle class. + * + * The class is created by the insertion of the simplex sigma. + * The methods adds a cocycle, representing the new cocycle class, + * to the matrix representing the cohomology groups. + * The new cocycle has value 0 on every simplex except on sigma + * where it worths 1.*/ + void create_cocycle(Simplex_handle sigma, Arith_element x, + Arith_element charac) { + Simplex_key key = cpx_->key(sigma); + // Create a column containing only one cell, + Column * new_col = column_pool_.construct(key); + Cell * new_cell = cell_pool_.construct(key, x, new_col); + new_col->col_.push_back(*new_cell); + // and insert it in the matrix, in constant time thanks to the hint cam_.end(). + // Indeed *new_col has the biggest lexicographic value because key is the + // biggest key used so far. + cam_.insert(cam_.end(), *new_col); + // Update the disjoint sets data structure. + Hcell * new_hcell = new Hcell; + new_hcell->push_back(*new_cell); + transverse_idx_[key] = cocycle(charac, new_hcell); // insert the new row + ds_repr_[key] = new_col; + } + + /* \brief Destroy a cocycle class. + * + * The cocycle class is destroyed by the insertion of sigma. + * The methods proceeds to a reduction of the matrix representing + * the cohomology groups using Gauss pivoting. The reduction zeros-out + * the row containing the cell with highest key in + * a_ds, the annotation of the boundary of simplex sigma. This key + * is "death_key".*/ + void destroy_cocycle(Simplex_handle sigma, A_ds_type const& a_ds, + Simplex_key death_key, Arith_element inv_x, + Arith_element charac) { + // Create a finite persistent interval for which the interval exists + if (interval_length_policy(cpx_->simplex(death_key), sigma)) { + persistent_pairs_.emplace_back(cpx_->simplex(death_key) // creator + , sigma // destructor + , charac); // fields + } + + auto death_key_row = transverse_idx_.find(death_key); // Find the beginning of the row. + std::pair<typename Cam::iterator, bool> result_insert_cam; + + auto row_cell_it = death_key_row->second.row_->begin(); + + while (row_cell_it != death_key_row->second.row_->end()) { // Traverse all cells in + // the row at index death_key. + Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_); + + if (w != coeff_field_.additive_identity()) { + Column * curr_col = row_cell_it->self_col_; + ++row_cell_it; + // Disconnect the column from the rows in the CAM. + for (auto& col_cell : curr_col->col_) { + col_cell.base_hook_cam_h::unlink(); + } + + // Remove the column from the CAM before modifying its value + cam_.erase(cam_.iterator_to(*curr_col)); + // Proceed to the reduction of the column + plus_equal_column(*curr_col, a_ds, w); + + if (curr_col->col_.empty()) { // If the column is null + ds_repr_[curr_col->class_key_] = NULL; + column_pool_.destroy(curr_col); // delete curr_col; + } else { + // Find whether the column obtained is already in the CAM + result_insert_cam = cam_.insert(*curr_col); + if (result_insert_cam.second) { // If it was not in the CAM before: insertion has succeeded + for (auto& col_cell : curr_col->col_) { + // re-establish the row links + transverse_idx_[col_cell.key_].row_->push_front(col_cell); + } + } else { // There is already an identical column in the CAM: + // merge two disjoint sets. + dsets_.link(curr_col->class_key_, + result_insert_cam.first->class_key_); + + Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_); + ds_repr_[key_tmp] = &(*(result_insert_cam.first)); + result_insert_cam.first->class_key_ = key_tmp; + // intrusive containers don't own their elements, we have to release them manually + curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);}); + column_pool_.destroy(curr_col); // delete curr_col; + } + } + } else { + ++row_cell_it; + } // If w == 0, pass. + } + + // Because it is a killer simplex, set the data of sigma to null_key(). + if (charac == coeff_field_.characteristic()) { + cpx_->assign_key(sigma, cpx_->null_key()); + } + if (death_key_row->second.characteristics_ == charac) { + delete death_key_row->second.row_; + transverse_idx_.erase(death_key_row); + } else { + death_key_row->second.characteristics_ /= charac; + } + } + + /* + * Assign: target <- target + w * other. + */ + void plus_equal_column(Column & target, A_ds_type const& other // value_type is pair<Simplex_key,Arith_element> + , Arith_element w) { + auto target_it = target.col_.begin(); + auto other_it = other.begin(); + while (target_it != target.col_.end() && other_it != other.end()) { + if (target_it->key_ < other_it->first) { + ++target_it; + } else { + if (target_it->key_ > other_it->first) { + Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first // key + , coeff_field_.additive_identity(), &target)); + + cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); + + target.col_.insert(target_it, *cell_tmp); + + ++other_it; + } else { // it1->key == it2->key + // target_it->coefficient_ <- target_it->coefficient_ + other_it->second * w + target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w); + if (target_it->coefficient_ == coeff_field_.additive_identity()) { + auto tmp_it = target_it; + ++target_it; + ++other_it; // iterators remain valid + Cell * tmp_cell_ptr = &(*tmp_it); + target.col_.erase(tmp_it); // removed from column + + cell_pool_.destroy(tmp_cell_ptr); // delete from memory + } else { + ++target_it; + ++other_it; + } + } + } + } + while (other_it != other.end()) { + Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.additive_identity(), &target)); + cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w); + target.col_.insert(target.col_.end(), *cell_tmp); + + ++other_it; + } + } + + /* + * Compare two intervals by length. + */ + struct cmp_intervals_by_length { + explicit cmp_intervals_by_length(FilteredComplex * sc) + : sc_(sc) { + } + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + } + FilteredComplex * sc_; + }; + + public: + /** \brief Output the persistence diagram in ostream. + * + * The file format is the following: + * p1*...*pr dim b d + * + * where "dim" is the dimension of the homological feature, + * b and d are respectively the birth and death of the feature and + * p1*...*pr is the product of prime numbers pi such that the homology + * 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; + 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<Filtration_value>::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)) << " " + << cpx_->filtration(get<0>(pair)) << " " + << cpx_->filtration(get<1>(pair)) << " " << std::endl; + } + } + } + + void write_output_diagram(std::string diagram_name) { + std::ofstream diagram_out(diagram_name.c_str()); + 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; + for (auto pair : persistent_pairs_) { + // Special case on windows, inf is "1.#INF" + if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) { + diagram_out << cpx_->dimension(get<0>(pair)) << " " + << cpx_->filtration(get<0>(pair)) << " inf" << std::endl; + } else { + diagram_out << cpx_->dimension(get<0>(pair)) << " " + << cpx_->filtration(get<0>(pair)) << " " + << cpx_->filtration(get<1>(pair)) << std::endl; + } + } + } + + /** @brief Returns Betti numbers. + * @return A vector of Betti numbers. + */ + std::vector<int> betti_numbers() const { + // Init Betti numbers vector with zeros until Simplicial complex dimension + std::vector<int> betti_numbers(dim_max_, 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 > to)\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(dim_max_, 0); + for (auto pair : persistent_pairs_) { + // Count persistence intervals that covers the given interval + // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will + // still work if we change the complex filtration function to reject null simplices. + if (cpx_->filtration(get<0>(pair)) <= from && + (get<1>(pair) == cpx_->null_simplex() || 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 > to)\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 + // null_simplex test : if the function is called with to=+infinity, we still get something useful. And it will + // still work if we change the complex filtration function to reject null simplices. + if (cpx_->filtration(get<0>(pair)) <= from && + (get<1>(pair) == cpx_->null_simplex() || cpx_->filtration(get<1>(pair)) > to)) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + // Increment betti number found + ++betti_number; + } + } + } + return betti_number; + } + + /** @brief Returns a list of persistence birth and death FilteredComplex::Simplex_handle pairs. + * @return A list of Persistent_cohomology::Persistent_interval + */ + const std::vector<Persistent_interval>& get_persistent_pairs() const { + return persistent_pairs_; + } + + /** @brief Returns persistence intervals for a given dimension. + * @param[in] dimension Dimension to get the birth and death pairs from. + * @return A vector of persistence intervals (birth and death) on a fixed dimension. + */ + std::vector< std::pair< Filtration_value , Filtration_value > > + intervals_in_dimension(int dimension) { + std::vector< std::pair< Filtration_value , Filtration_value > > result; + // auto && pair, to avoid unnecessary copying + for (auto && pair : persistent_pairs_) { + if (cpx_->dimension(get<0>(pair)) == dimension) { + result.emplace_back(cpx_->filtration(get<0>(pair)), cpx_->filtration(get<1>(pair))); + } + } + return result; + } + + private: + /* + * Structure representing a cocycle. + */ + struct cocycle { + cocycle() + : row_(nullptr), + characteristics_() { + } + cocycle(Arith_element characteristics, Hcell * row) + : row_(row), + characteristics_(characteristics) { + } + + Hcell * row_; // points to the corresponding row in the CAM + Arith_element characteristics_; // product of field characteristics for which the cocycle exist + }; + + public: + FilteredComplex * cpx_; + int dim_max_; + CoefficientField coeff_field_; + size_t num_simplices_; + + /* Disjoint sets data structure to link the model of FilteredComplex + * with the compressed annotation matrix. + * ds_rank_ is a property map Simplex_key -> int, ds_parent_ is a property map + * Simplex_key -> simplex_key_t */ + std::vector<int> ds_rank_; + std::vector<Simplex_key> ds_parent_; + std::vector<Column *> ds_repr_; + boost::disjoint_sets<int *, Simplex_key *> dsets_; + /* The compressed annotation matrix fields.*/ + Cam cam_; + /* Dictionary establishing the correspondance between the Simplex_key of + * the root vertex in the union-find ds and the Simplex_key of the vertex which + * created the connected component as a 0-dimension homology feature.*/ + std::map<Simplex_key, Simplex_key> zero_cocycles_; + /* Key -> row. */ + std::map<Simplex_key, cocycle> transverse_idx_; + /* Persistent intervals. */ + std::vector<Persistent_interval> persistent_pairs_; + length_interval interval_length_policy; + + Simple_object_pool<Column> column_pool_; + Simple_object_pool<Cell> cell_pool_; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h new file mode 100644 index 00000000..0673625c --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -0,0 +1,104 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ +#define PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ + +#include <utility> +#include <vector> + +namespace Gudhi { + +namespace persistent_cohomology { + +/** \brief Structure representing the coefficient field \f$\mathbb{Z}/p\mathbb{Z}\f$ + * + * \implements CoefficientField + * \ingroup persistent_cohomology + */ +class Field_Zp { + public: + typedef int Element; + + Field_Zp() + : Prime(0), + inverse_() { + } + + void init(int charac) { + assert(charac > 0); // division by zero + non negative values + Prime = charac; + inverse_.clear(); + inverse_.reserve(charac); + inverse_.push_back(0); + for (int i = 1; i < Prime; ++i) { + int inv = 1; + while (((inv * i) % Prime) != 1) + ++inv; + inverse_.push_back(inv); + } + } + + /** Set x <- x + w * y*/ + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(Prime > 0); // division by zero + non negative values + Element result = (x + w * y) % Prime; + if (result < 0) + result += Prime; + return result; + } + +// operator= defined on Element + + /** Returns y * w */ + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, (Element)w); + } + + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + Element additive_identity() const { + return 0; + } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + Element multiplicative_identity(Element = 0) const { + return 1; + } + /** Returns the inverse in the field. Modifies P. ??? */ + std::pair<Element, Element> inverse(Element x, Element P) { + return std::pair<Element, Element>(inverse_[x], P); + } // <------ return the product of field characteristic for which x is invertible + + /** Returns -x * y.*/ + Element times_minus(Element x, Element y) { + assert(Prime > 0); // division by zero + non negative values + Element out = (-x * y) % Prime; + return (out < 0) ? out + Prime : out; + } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + int characteristic() const { + return Prime; + } + + private: + int Prime; + /** Property map Element -> Element, which associate to an element its inverse in the field.*/ + std::vector<Element> inverse_; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h new file mode 100644 index 00000000..1754a2ec --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -0,0 +1,173 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ +#define PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ + +#include <gmpxx.h> + +#include <vector> +#include <utility> + +namespace Gudhi { + +namespace persistent_cohomology { + +/** \brief Structure representing coefficients in a set of finite fields simultaneously + * using the chinese remainder theorem. + * + * \implements CoefficientField + * \ingroup persistent_cohomology + + * Details on the algorithms may be found in \cite boissonnat:hal-00922572 + */ +class Multi_field { + public: + typedef mpz_class Element; + + Multi_field() + : prod_characteristics_(0), + mult_id_all(0), + add_id_all(0) { + } + + /* Initialize the multi-field. The generation of prime numbers might fail with + * a very small probability.*/ + void init(int min_prime, int max_prime) { + if (max_prime < 2) { + std::cerr << "There is no prime less than " << max_prime << std::endl; + } + if (min_prime > max_prime) { + std::cerr << "No prime in [" << min_prime << ":" << max_prime << "]" + << std::endl; + } + // fill the list of prime numbers + int curr_prime = min_prime; + mpz_t tmp_prime; + mpz_init_set_ui(tmp_prime, min_prime); + // test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime, 25); // probabilistic primality test + + if (is_prime == 0) { // min_prime is composite + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + + while (curr_prime <= max_prime) { + primes_.push_back(curr_prime); + mpz_nextprime(tmp_prime, tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + mpz_clear(tmp_prime); + // set m to primorial(bound_prime) + prod_characteristics_ = 1; + for (auto p : primes_) { + prod_characteristics_ *= p; + } + + // Uvect_ + Element Ui; + Element tmp_elem; + for (auto p : primes_) { + assert(p > 0); // division by zero + non negative values + tmp_elem = prod_characteristics_ / p; + // Element tmp_elem_bis = 10; + mpz_powm_ui(tmp_elem.get_mpz_t(), tmp_elem.get_mpz_t(), p - 1, + prod_characteristics_.get_mpz_t()); + Uvect_.push_back(tmp_elem); + } + mult_id_all = 0; + for (auto uvect : Uvect_) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + mult_id_all = (mult_id_all + uvect) % prod_characteristics_; + } + } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + const Element& additive_identity() const { + return add_id_all; + } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + const Element& multiplicative_identity() const { + return mult_id_all; + } // 1 everywhere + + Element multiplicative_identity(Element Q) { + if (Q == prod_characteristics_) { + return multiplicative_identity(); + } + + assert(prod_characteristics_ > 0); // division by zero + non negative values + Element mult_id = 0; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + assert(primes_[idx] > 0); // division by zero + non negative values + if ((Q % primes_[idx]) == 0) { + mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; + } + } + return mult_id; + } + + /** Returns y * w */ + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, w); + } + + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + const Element& characteristic() const { + return prod_characteristics_; + } + + /** Returns the inverse in the field. Modifies P. ??? */ + std::pair<Element, Element> inverse(Element x, Element QS) { + Element QR; + mpz_gcd(QR.get_mpz_t(), x.get_mpz_t(), QS.get_mpz_t()); // QR <- gcd(x,QS) + if (QR == QS) + return std::pair<Element, Element>(additive_identity(), multiplicative_identity()); // partial inverse is 0 + Element QT = QS / QR; + Element inv_qt; + mpz_invert(inv_qt.get_mpz_t(), x.get_mpz_t(), QT.get_mpz_t()); + + assert(prod_characteristics_ > 0); // division by zero + non negative values + return { (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT }; + } + /** Returns -x * y.*/ + Element times_minus(const Element& x, const Element& y) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + /* This assumes that (x*y)%pc cannot be zero, but Field_Zp has specific code for the 0 case ??? */ + return prod_characteristics_ - ((x * y) % prod_characteristics_); + } + + /** Set x <- x + w * y*/ + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + Element result = (x + w * y) % prod_characteristics_; + if (result < 0) + result += prod_characteristics_; + return result; + } + + Element prod_characteristics_; // product of characteristics of the fields + // represented by the multi-field class + std::vector<int> primes_; // all the characteristics of the fields + std::vector<Element> Uvect_; + Element mult_id_all; + const Element add_id_all; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h new file mode 100644 index 00000000..480be389 --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h @@ -0,0 +1,128 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ +#define PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ + +#include <boost/intrusive/set.hpp> +#include <boost/intrusive/list.hpp> + +#include <list> + +namespace Gudhi { + +namespace persistent_cohomology { + +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_column; + +struct cam_h_tag; +// for horizontal traversal in the CAM +struct cam_v_tag; +// for vertical traversal in the CAM + +typedef boost::intrusive::list_base_hook<boost::intrusive::tag<cam_h_tag>, + boost::intrusive::link_mode<boost::intrusive::auto_unlink> // allows .unlink() +> base_hook_cam_h; + +typedef boost::intrusive::list_base_hook<boost::intrusive::tag<cam_v_tag>, + boost::intrusive::link_mode<boost::intrusive::normal_link> // faster hook, less safe +> base_hook_cam_v; + +/** \internal + * \brief + * + */ +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_cell : public base_hook_cam_h, + public base_hook_cam_v { + public: + template<class T1, class T2> friend class Persistent_cohomology; + friend class Persistent_cohomology_column<SimplexKey, ArithmeticElement>; + + typedef Persistent_cohomology_column<SimplexKey, ArithmeticElement> Column; + + Persistent_cohomology_cell(SimplexKey key, ArithmeticElement x, + Column * self_col) + : key_(key), + coefficient_(x), + self_col_(self_col) { + } + + SimplexKey key_; + ArithmeticElement coefficient_; + Column * self_col_; +}; + +/* + * \brief Sparse column for the Compressed Annotation Matrix. + * + * The non-zero coefficients of the column are stored in a + * boost::intrusive::list. Contains a hook to be stored in a + * boost::intrusive::set. + * + * Movable but not Copyable. + */ +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_column : public boost::intrusive::set_base_hook< + boost::intrusive::link_mode<boost::intrusive::normal_link> > { + template<class T1, class T2> friend class Persistent_cohomology; + + public: + typedef Persistent_cohomology_cell<SimplexKey, ArithmeticElement> Cell; + typedef boost::intrusive::list<Cell, + boost::intrusive::constant_time_size<false>, + boost::intrusive::base_hook<base_hook_cam_v> > Col_type; + + /** \brief Creates an empty column.*/ + explicit Persistent_cohomology_column(SimplexKey key) + : col_(), + class_key_(key) {} + + /** \brief Returns true iff the column is null.*/ + bool is_null() const { + return col_.empty(); + } + /** \brief Returns the key of the representative simplex of + * the set of simplices having this column as annotation vector + * in the compressed annotation matrix.*/ + SimplexKey class_key() const { + return class_key_; + } + + /** \brief Lexicographic comparison of two columns.*/ + friend bool operator<(const Persistent_cohomology_column& c1, + const Persistent_cohomology_column& c2) { + typename Col_type::const_iterator it1 = c1.col_.begin(); + typename Col_type::const_iterator it2 = c2.col_.begin(); + while (it1 != c1.col_.end() && it2 != c2.col_.end()) { + if (it1->key_ == it2->key_) { + if (it1->coefficient_ == it2->coefficient_) { + ++it1; + ++it2; + } else { + return it1->coefficient_ < it2->coefficient_; + } + } else { + return it1->key_ < it2->key_; + } + } + return (it2 != c2.col_.end()); + } + + Col_type col_; + SimplexKey class_key_; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ diff --git a/src/Persistent_cohomology/test/CMakeLists.txt b/src/Persistent_cohomology/test/CMakeLists.txt new file mode 100644 index 00000000..f8baf861 --- /dev/null +++ b/src/Persistent_cohomology/test/CMakeLists.txt @@ -0,0 +1,37 @@ +project(Persistent_cohomology_tests) + +include(GUDHI_test_coverage) + +add_executable ( Persistent_cohomology_test_unit persistent_cohomology_unit_test.cpp ) +target_link_libraries(Persistent_cohomology_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( Persistent_cohomology_test_betti_numbers betti_numbers_unit_test.cpp ) +target_link_libraries(Persistent_cohomology_test_betti_numbers ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(Persistent_cohomology_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Persistent_cohomology_test_betti_numbers ${TBB_LIBRARIES}) +endif(TBB_FOUND) + +# Do not forget to copy test results files in current binary dir +file(COPY "${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt" + DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +# Unitary tests +gudhi_add_coverage_test(Persistent_cohomology_test_unit) +gudhi_add_coverage_test(Persistent_cohomology_test_betti_numbers) + +if(GMPXX_FOUND AND GMP_FOUND) + add_executable ( Persistent_cohomology_test_unit_multi_field persistent_cohomology_unit_test_multi_field.cpp ) + target_link_libraries(Persistent_cohomology_test_unit_multi_field + ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(Persistent_cohomology_test_unit_multi_field ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + # Do not forget to copy test results files in current binary dir + file(COPY "${CMAKE_SOURCE_DIR}/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt" + DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + # Unitary tests + gudhi_add_coverage_test(Persistent_cohomology_test_unit_multi_field) + +endif(GMPXX_FOUND AND GMP_FOUND) + diff --git a/src/Persistent_cohomology/test/README b/src/Persistent_cohomology/test/README new file mode 100644 index 00000000..0c41feed --- /dev/null +++ b/src/Persistent_cohomology/test/README @@ -0,0 +1,29 @@ +To compile: +*********** + +cd /path-to-gudhi/ +cmake . +cd /path-to-test/ +make + +To launch with details: +*********************** + +SINGLE FIELD +------------ +./Persistent_cohomology_test_unit --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) + +MULTI FIELD +----------- +./Persistent_cohomology_test_unit_multi_field --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) + +BETTI NUMBERS +------------- +./Persistent_cohomology_test_betti_numbers --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) + 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..0a08d200 --- /dev/null +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -0,0 +1,287 @@ +#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 MiniSTOptions : 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<MiniSTOptions>; +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); + + // 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. + + std::cout << "BETTI NUMBERS" << std::endl; + + 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); + + std::cout << "GET PERSISTENT PAIRS" << std::endl; + + // 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()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 0); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 0); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 0); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); +} + +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); + + // 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()); + + std::cout << "INTERVALS IN DIMENSION" << std::endl; + + auto intervals_in_dimension_0 = pcoh.intervals_in_dimension(0); + std::cout << "intervals_in_dimension_0.size() = " << intervals_in_dimension_0.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_0.size(); i++) + std::cout << "intervals_in_dimension_0[" << i << "] = [" << intervals_in_dimension_0[i].first << "," << + intervals_in_dimension_0[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_0.size() == 2); + BOOST_CHECK(intervals_in_dimension_0[0].first == 2); + BOOST_CHECK(intervals_in_dimension_0[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + BOOST_CHECK(intervals_in_dimension_0[1].first == 1); + BOOST_CHECK(intervals_in_dimension_0[1].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_1 = pcoh.intervals_in_dimension(1); + std::cout << "intervals_in_dimension_1.size() = " << intervals_in_dimension_1.size() << std::endl; + for (std::size_t i = 0; i < intervals_in_dimension_1.size(); i++) + std::cout << "intervals_in_dimension_1[" << i << "] = [" << intervals_in_dimension_1[i].first << "," << + intervals_in_dimension_1[i].second << "]" << std::endl; + BOOST_CHECK(intervals_in_dimension_1.size() == 1); + BOOST_CHECK(intervals_in_dimension_1[0].first == 4); + BOOST_CHECK(intervals_in_dimension_1[0].second == std::numeric_limits<Mini_simplex_tree::Filtration_value>::infinity()); + + auto intervals_in_dimension_2 = pcoh.intervals_in_dimension(2); + std::cout << "intervals_in_dimension_2.size() = " << intervals_in_dimension_2.size() << std::endl; + BOOST_CHECK(intervals_in_dimension_2.size() == 0); +} diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp new file mode 100644 index 00000000..a1c106d5 --- /dev/null +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -0,0 +1,212 @@ +#include <iostream> +#include <string> +#include <algorithm> +#include <utility> // std::pair, std::make_pair +#include <cmath> // float comparison +#include <limits> +#include <cstdint> // for std::uint8_t + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "persistent_cohomology" +#include <boost/test/unit_test.hpp> + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; +using namespace boost::unit_test; + +typedef Simplex_tree<> typeST; + +std::string test_rips_persistence(int coefficient, int min_persistence) { + // file is copied in CMakeLists.txt + std::ifstream simplex_tree_stream; + simplex_tree_stream.open("simplex_tree_file_for_unit_test.txt"); + typeST st; + simplex_tree_stream >> st; + simplex_tree_stream.close(); + + // Display the Simplex_tree + std::cout << "The complex contains " << st.num_simplices() << " simplices" << " - dimension= " << st.dimension() + << std::endl; + + // Check + BOOST_CHECK(st.num_simplices() == 98); + BOOST_CHECK(st.dimension() == 3); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology<Simplex_tree<>, Field_Zp> pcoh(st); + + pcoh.init_coefficients( coefficient ); // initializes the coefficient field for homology + // Check infinite rips + pcoh.compute_persistent_cohomology( min_persistence ); // Minimal lifetime of homology feature to be recorded. + std::ostringstream ossInfinite; + + pcoh.output_diagram(ossInfinite); + std::string strInfinite = ossInfinite.str(); + return strInfinite; +} + +void test_rips_persistence_in_dimension(int dimension) { + std::string value0(" 0 0.02 1.12"); + std::string value1(" 0 0.03 1.13"); + std::string value2(" 0 0.04 1.14"); + std::string value3(" 0 0.05 1.15"); + std::string value4(" 0 0.06 1.16"); + std::string value5(" 0 0.07 1.17"); + std::string value6(" 0 0.08 1.18"); + std::string value7(" 0 0.09 1.19"); + std::string value8(" 0 0 inf" ); + std::string value9(" 0 0.01 inf" ); + + value0.insert(0,std::to_string(dimension)); + value1.insert(0,std::to_string(dimension)); + value2.insert(0,std::to_string(dimension)); + value3.insert(0,std::to_string(dimension)); + value4.insert(0,std::to_string(dimension)); + value5.insert(0,std::to_string(dimension)); + value6.insert(0,std::to_string(dimension)); + value7.insert(0,std::to_string(dimension)); + value8.insert(0,std::to_string(dimension)); + value9.insert(0,std::to_string(dimension)); + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=0" << std::endl; + + std::string str_rips_persistence = test_rips_persistence(dimension, 0); + std::cout << str_rips_persistence << std::endl; + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=1" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, 1); + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=2" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, 2); + + BOOST_CHECK(str_rips_persistence.find(value0) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value1) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value2) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_SINGLE_FIELD DIM=" << dimension << " MIN_PERS=Inf" << std::endl; + + str_rips_persistence = test_rips_persistence(dimension, (std::numeric_limits<int>::max)()); + + BOOST_CHECK(str_rips_persistence.find(value0) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value1) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value2) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value8) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value9) != std::string::npos); // Check found + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_1 ) +{ + test_rips_persistence_in_dimension(1); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_2 ) +{ + test_rips_persistence_in_dimension(2); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_3 ) +{ + test_rips_persistence_in_dimension(3); +} + +BOOST_AUTO_TEST_CASE( rips_persistent_cohomology_single_field_dim_5 ) +{ + test_rips_persistence_in_dimension(5); +} + +// TODO(VR): not working from 6 +// std::string str_rips_persistence = test_rips_persistence(6, 0); +// TODO(VR): division by zero +// std::string str_rips_persistence = test_rips_persistence(0, 0); + +/** SimplexTree minimal options to test the limits. + * + * Maximum number of simplices to compute persistence is <CODE>std::numeric_limits<std::uint8_t>::max()<\CODE> = 256.*/ +struct MiniSTOptions { + typedef linear_indexing_tag Indexing_tag; + typedef short Vertex_handle; + typedef double Filtration_value; + // Maximum number of simplices to compute persistence is 2^8 - 1 = 255. One is reserved for null_key + typedef std::uint8_t Simplex_key; + static const bool store_key = true; + static const bool store_filtration = false; + static const bool contiguous_vertices = false; +}; + +using Mini_simplex_tree = Gudhi::Simplex_tree<MiniSTOptions>; +using Mini_st_persistence = + Gudhi::persistent_cohomology::Persistent_cohomology<Mini_simplex_tree, Gudhi::persistent_cohomology::Field_Zp>; + +BOOST_AUTO_TEST_CASE( persistence_constructor_exception ) +{ + Mini_simplex_tree st; + + // To make number of simplices = 255 + const short simplex_0[] = {0, 1, 2, 3, 4, 5, 6, 7}; + st.insert_simplex_and_subfaces(simplex_0); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + BOOST_CHECK(st.num_simplices() <= std::numeric_limits<MiniSTOptions::Simplex_key>::max()); + // Class for homology computation + BOOST_CHECK_NO_THROW(Mini_st_persistence pcoh(st)); + + st.insert_simplex({8}); + BOOST_CHECK(st.num_simplices() > std::numeric_limits<MiniSTOptions::Simplex_key>::max()); + // Class for homology computation + BOOST_CHECK_THROW(Mini_st_persistence pcoh2(st), std::out_of_range); + +} diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp new file mode 100644 index 00000000..9e767943 --- /dev/null +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp @@ -0,0 +1,115 @@ +#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 "persistent_cohomology_multi_field" +#include <boost/test/unit_test.hpp> + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Persistent_cohomology/Multi_field.h> + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; +using namespace boost::unit_test; + +typedef Simplex_tree<> typeST; + +std::string test_rips_persistence(int min_coefficient, int max_coefficient, double min_persistence) { + // file is copied in CMakeLists.txt + std::ifstream simplex_tree_stream; + simplex_tree_stream.open("simplex_tree_file_for_multi_field_unit_test.txt"); + typeST st; + simplex_tree_stream >> st; + simplex_tree_stream.close(); + + // Display the Simplex_tree + std::cout << "The complex contains " << st.num_simplices() << " simplices" << " - dimension= " << st.dimension() + << std::endl; + + // Check + BOOST_CHECK(st.num_simplices() == 58); + BOOST_CHECK(st.dimension() == 3); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology<Simplex_tree<>, Multi_field> pcoh(st); + + pcoh.init_coefficients(min_coefficient, max_coefficient); // initializes the coefficient field for homology + // Check infinite rips + pcoh.compute_persistent_cohomology(min_persistence); // Minimal lifetime of homology feature to be recorded. + + std::ostringstream ossRips; + pcoh.output_diagram(ossRips); + + std::string strRips = ossRips.str(); + return strRips; +} + +void test_rips_persistence_in_dimension(int min_dimension, int max_dimension) { + // there are 2 discontinued ensembles + std::string value0(" 0 0.25 inf"); + std::string value1(" 1 0.4 inf"); + // And a big hole - cut in 2 pieces after 0.3 + std::string value2(" 0 0.2 0.3"); + + // For dim <= 1 => + std::string value3(" 1 0.25 inf"); + std::string value4(" 2 0.25 inf"); + std::string value5(" 1 0.3 inf"); + std::string value6(" 2 0.3 inf"); + std::string value7(" 2 0.4 inf"); + + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF RIPS_PERSISTENT_COHOMOLOGY_MULTI_FIELD MIN_DIM=" << min_dimension << " MAX_DIM=" << max_dimension << " MIN_PERS=0" << std::endl; + + std::string str_rips_persistence = test_rips_persistence(min_dimension, max_dimension, 0.0); + std::cout << "str_rips_persistence=" << str_rips_persistence << std::endl; + + BOOST_CHECK(str_rips_persistence.find(value0) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value1) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value2) != std::string::npos); // Check found + + if ((min_dimension < 2) && (max_dimension < 2)) { + BOOST_CHECK(str_rips_persistence.find(value3) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value4) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value5) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value6) != std::string::npos); // Check found + BOOST_CHECK(str_rips_persistence.find(value7) != std::string::npos); // Check found + } else { + BOOST_CHECK(str_rips_persistence.find(value3) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value4) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value5) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value6) == std::string::npos); // Check not found + BOOST_CHECK(str_rips_persistence.find(value7) == std::string::npos); // Check not found + } + +} + +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_1_2) { + test_rips_persistence_in_dimension(0, 1); +} + +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_2_3) { + test_rips_persistence_in_dimension(1, 3); +} + +BOOST_AUTO_TEST_CASE(rips_persistent_cohomology_multi_field_dim_1_5) { + test_rips_persistence_in_dimension(1, 5); +} + +// TODO(VR): not working from 6 +// std::string str_rips_persistence = test_rips_persistence(6, 0); +// TODO(VR): division by zero +// std::string str_rips_persistence = test_rips_persistence(0, 0); +// TODO(VR): is result OK of : +// test_rips_persistence_in_dimension(3, 4); + diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt new file mode 100644 index 00000000..ed2c0c3d --- /dev/null +++ b/src/Persistent_cohomology/test/simplex_tree_file_for_multi_field_unit_test.txt @@ -0,0 +1,58 @@ +0 0 0.2 +0 3 0.2 +1 3 0 0.2 +0 6 0.2 +0 11 0.2 +1 11 6 0.2 +0 13 0.25 +0 14 0.25 +1 14 13 0.25 +0 15 0.25 +1 15 13 0.25 +1 15 14 0.25 +2 15 14 13 0.25 +0 1 0.3 +1 1 0 0.3 +0 2 0.3 +1 2 0 0.3 +1 2 1 0.3 +2 2 1 0 0.3 +0 4 0.3 +1 4 3 0.3 +0 5 0.3 +1 5 3 0.3 +1 5 4 0.3 +2 5 4 3 0.3 +0 9 0.3 +0 10 0.3 +1 10 2 0.3 +1 10 9 0.3 +1 11 9 0.3 +1 11 10 0.3 +2 11 10 9 0.3 +0 12 0.3 +1 12 2 0.3 +1 12 10 0.3 +2 12 10 2 0.3 +1 6 0 0.4 +1 6 1 0.4 +2 6 1 0 0.4 +0 7 0.4 +1 7 0 0.4 +1 7 1 0.4 +2 7 1 0 0.4 +1 7 6 0.4 +2 7 6 0 0.4 +2 7 6 1 0.4 +3 7 6 1 0 0.4 +0 8 0.4 +1 8 4 0.4 +1 8 5 0.4 +2 8 5 4 0.4 +1 9 4 0.4 +1 9 5 0.4 +2 9 5 4 0.4 +1 9 8 0.4 +2 9 8 4 0.4 +2 9 8 5 0.4 +3 9 8 5 4 0.4 diff --git a/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt b/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt new file mode 100644 index 00000000..90756ce0 --- /dev/null +++ b/src/Persistent_cohomology/test/simplex_tree_file_for_unit_test.txt @@ -0,0 +1,98 @@ +0 0 0 +0 1 0.01 +0 2 0.02 +0 3 0.03 +0 4 0.04 +0 5 0.05 +0 6 0.06 +0 7 0.07 +0 8 0.08 +0 9 0.09 +1 2 1 1.12 +1 3 1 1.13 +1 4 1 1.14 +1 5 1 1.15 +1 6 1 1.16 +1 7 1 1.17 +1 8 1 1.18 +1 9 1 1.19 +1 3 2 1.23 +2 3 2 1 1.23 +1 6 2 1.26 +2 6 2 1 1.26 +1 8 2 1.28 +2 8 2 1 1.28 +1 9 2 1.29 +2 9 2 1 1.29 +1 6 3 1.36 +2 6 3 1 1.36 +2 6 3 2 1.36 +3 6 3 2 1 1.36 +1 7 3 1.37 +2 7 3 1 1.37 +1 8 3 1.38 +2 8 3 1 1.38 +2 8 3 2 1.38 +3 8 3 2 1 1.38 +1 5 4 1.45 +2 5 4 1 1.45 +1 6 4 1.46 +2 6 4 1 1.46 +1 8 4 1.48 +2 8 4 1 1.48 +1 6 5 1.56 +2 6 5 1 1.56 +2 6 5 4 1.56 +3 6 5 4 1 1.56 +1 7 5 1.57 +2 7 5 1 1.57 +1 8 5 1.58 +2 8 5 1 1.58 +2 8 5 4 1.58 +3 8 5 4 1 1.58 +1 9 5 1.59 +2 9 5 1 1.59 +1 7 6 1.67 +2 7 6 1 1.67 +2 7 6 3 1.67 +3 7 6 3 1 1.67 +2 7 6 5 1.67 +3 7 6 5 1 1.67 +1 8 6 1.68 +2 8 6 1 1.68 +2 8 6 2 1.68 +3 8 6 2 1 1.68 +2 8 6 3 1.68 +3 8 6 3 1 1.68 +3 8 6 3 2 1.68 +2 8 6 4 1.68 +3 8 6 4 1 1.68 +2 8 6 5 1.68 +3 8 6 5 1 1.68 +3 8 6 5 4 1.68 +1 9 6 1.69 +2 9 6 1 1.69 +2 9 6 2 1.69 +3 9 6 2 1 1.69 +2 9 6 5 1.69 +3 9 6 5 1 1.69 +1 8 7 1.78 +2 8 7 1 1.78 +2 8 7 3 1.78 +3 8 7 3 1 1.78 +2 8 7 5 1.78 +3 8 7 5 1 1.78 +2 8 7 6 1.78 +3 8 7 6 1 1.78 +3 8 7 6 3 1.78 +3 8 7 6 5 1.78 +1 9 8 1.89 +2 9 8 1 1.89 +2 9 8 2 1.89 +3 9 8 2 1 1.89 +2 9 8 5 1.89 +3 9 8 5 1 1.89 +2 9 8 6 1.89 +3 9 8 6 1 1.89 +3 9 8 6 2 1.89 +3 9 8 6 5 1.89 |
