diff options
Diffstat (limited to 'src/Persistent_cohomology')
13 files changed, 2109 insertions, 0 deletions
diff --git a/src/Persistent_cohomology/concept/CoefficientField.h b/src/Persistent_cohomology/concept/CoefficientField.h new file mode 100644 index 00000000..953b06c2 --- /dev/null +++ b/src/Persistent_cohomology/concept/CoefficientField.h @@ -0,0 +1,52 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/** \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..1834903b --- /dev/null +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -0,0 +1,143 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/** \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 +{ +/** Handle to specify a simplex. */ + typedef unspecified Simplex_handle; +/** \brief Key associated to each simplex. + * + * Must be a signed integer type. */ + typedef unspecified Simplex_key; +/** \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; + +/** Returns a Simplex_hanlde 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 a key that is different from the keys associated + * to the simplices. */ + Simplex_key null_key (); +/** \brief Returns the key associated to a simplex. */ + Simplex_key key ( Simplex_handle sh ); +/** \brief Returns the simplex associated to a key. + * + * If key is different from null_key(), there must be a unique + * simplex having this key. */ + Simplex_handle simplex ( Simplex_key key ); +/** \brief Assign a key to a simplex. */ + void assign_key(Simplex_handle sh, Simplex_key key); + +/** \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(); + + +/* \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..111723a5 --- /dev/null +++ b/src/Persistent_cohomology/concept/PersistentHomology.h @@ -0,0 +1,41 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +/** \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/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt new file mode 100644 index 00000000..dca4a04d --- /dev/null +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -0,0 +1,22 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIExPersCohom) + +# NEED TO FIND BOOST NEEDED COMPONENTS TO LINK WITH +find_package(Boost 1.45.0 COMPONENTS system program_options) +message("Boost_lib = ${Boost_LIBRARIES}") + +add_executable ( rips_persistence rips_persistence.cpp ) +target_link_libraries(rips_persistence ${Boost_LIBRARIES}) +add_executable ( rips_persistence_from_alpha_shape_3 rips_persistence_from_alpha_shape_3.cpp ) +target_link_libraries(rips_persistence_from_alpha_shape_3 ${Boost_LIBRARIES}) + +if(GMPXX_FOUND AND GMP_FOUND) + message("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") + message("GMP_LIBRARIES = ${GMP_LIBRARIES}") + + add_executable ( rips_multifield_persistence rips_multifield_persistence.cpp ) + target_link_libraries(rips_multifield_persistence ${Boost_LIBRARIES} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + + add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) + target_link_libraries(performance_rips_persistence ${Boost_LIBRARIES} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) +endif() diff --git a/src/Persistent_cohomology/example/README b/src/Persistent_cohomology/example/README new file mode 100644 index 00000000..5dd34b85 --- /dev/null +++ b/src/Persistent_cohomology/example/README @@ -0,0 +1,47 @@ +To build the example, run in a Terminal: + +cd /path-to-example/ +cmake . +make + + +Example of use : + +Computation of the persistent homology with Z/2Z coefficients of the Rips complex on points +sampling a Klein bottle: + +./rips_persistence ../../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -m 100 + +output: +210 0 0 inf +210 1 0.0702103 inf +2 1 0.0702103 inf +2 2 0.159992 inf + +with Z/3Z coefficients: + +./rips_persistence ../../../data/points/Kl.txt -r 0.25 -d 3 -p 3 -m 100 + +output: +3 0 0 inf +3 1 0.0702103 inf + +and the computation with Z/2Z and Z/3Z coefficients simultaneously: + +./rips_multifield_persistence ../../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 3 -m 100 + +output: +6 0 0 inf +6 1 0.0702103 inf +2 1 0.0702103 inf +2 2 0.159992 inf + +and finally the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers): + + ./rips_multifield_persistence ../../../data/points/Kl.txt -r 0.25 -d 3 -p 2 -q 71 -m 100 + +output: +557940830126698960967415390 0 0 inf +557940830126698960967415390 1 0.0702103 inf +2 1 0.0702103 inf +2 2 0.159992 inf diff --git a/src/Persistent_cohomology/example/performance_rips_persistence.cpp b/src/Persistent_cohomology/example/performance_rips_persistence.cpp new file mode 100644 index 00000000..37d41ee1 --- /dev/null +++ b/src/Persistent_cohomology/example/performance_rips_persistence.cpp @@ -0,0 +1,168 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#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 "gudhi/Persistent_cohomology/Multi_field.h" +#include "gudhi/Hasse_complex.h" + +#include <chrono> + +using namespace Gudhi; + +/* 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 enlapsed_sec; + + std::string filepoints = "../examples/Kl.txt"; + Filtration_value threshold = 0.3; + int dim_max = 3; + int p = 2; + int q = 1223; + +// Extract the points from the file filepoints + typedef std::vector<double> Point_t; + std::vector< Point_t > points; + read_points( filepoints, points ); + +// Compute the proximity graph of the points + start = std::chrono::system_clock::now(); + Graph_t prox_graph = compute_proximity_graph( points, threshold + , euclidean_distance<Point_t> ); + end = std::chrono::system_clock::now(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << "Compute Rips graph in " << enlapsed_sec << " sec.\n"; + +// Construct the Rips complex in a Simplex Tree + Simplex_tree<> st; + start = std::chrono::system_clock::now(); + + st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree + st.expansion( dim_max ); // expand the graph until dimension dim_max + + end = std::chrono::system_clock::now(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << "Compute Rips complex in " << enlapsed_sec << " sec.\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(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << "Order the simplices of the filtration in " << enlapsed_sec << " sec.\n"; + +// Convert the simplex tree into a hasse diagram + start = std::chrono::system_clock::now(); + Hasse_complex<> hcpx(st); + end = std::chrono::system_clock::now(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << "Convert the simplex tree into a Hasse diagram in " << enlapsed_sec << " sec.\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); + + return 0; +} + + +template< typename FilteredComplex> +void +timing_persistence( FilteredComplex & cpx + , int p ) +{ + std::chrono::time_point<std::chrono::system_clock> start, end; + int enlapsed_sec; + + Persistent_cohomology< FilteredComplex, Field_Zp > pcoh (cpx); + pcoh.init_coefficients( p ); //initilizes the coefficient field for homology + + start = std::chrono::system_clock::now(); + + pcoh.compute_persistent_cohomology( INFINITY ); + + end = std::chrono::system_clock::now(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << " Compute persistent homology in Z/"<<p<<"Z in " << enlapsed_sec << " sec.\n"; +} + +template< typename FilteredComplex> +void +timing_persistence( FilteredComplex & cpx + , int p + , int q ) +{ + std::chrono::time_point<std::chrono::system_clock> start, end; + int enlapsed_sec; + + Persistent_cohomology< FilteredComplex, Multi_field > pcoh (cpx); + pcoh.init_coefficients( p, q ); //initilizes the coefficient field for homology + // 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(); + enlapsed_sec = std::chrono::duration_cast<std::chrono::seconds>(end-start).count(); + std::cout << " Compute multi-field persistent homology in all coefficient fields Z/pZ " + << "with p in ["<<p<<";"<<q<<"] in " << enlapsed_sec << " sec.\n"; +} 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..95ef8809 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -0,0 +1,152 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#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 "gudhi/Persistent_cohomology/Multi_field.h" + +#include <boost/program_options.hpp> + +using namespace Gudhi; + +typedef int Vertex_handle; +typedef double Filtration_value; + +void program_options( int argc, char * argv[] + , std::string & filepoints + , 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 filepoints; + std::string filediag; + Filtration_value threshold; + int dim_max; + int min_p; + int max_p; + Filtration_value min_persistence; + + program_options(argc,argv,filepoints,filediag,threshold,dim_max,min_p,max_p,min_persistence); + +// Extract the points from the file filepoints + typedef std::vector<double> Point_t; + std::vector< Point_t > points; + read_points( filepoints, points ); + +// Compute the proximity graph of the points + Graph_t prox_graph = compute_proximity_graph( points, threshold + , euclidean_distance<Point_t> ); + +// Construct the Rips complex in a Simplex Tree + Simplex_tree<> st; + st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree + st.expansion( dim_max ); // expand the graph until dimension dim_max + +// 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_p, max_p ); //initilizes the coefficient field for homology + // compute persistent homology, disgarding persistent features of life shorter than min_persistence + 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 & filepoints + , 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>(&filepoints), + "Name of file containing a point set. Format is one point per line: X1 ... Xd \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; + std::abort(); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp new file mode 100644 index 00000000..d7927223 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -0,0 +1,146 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#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> + +using namespace Gudhi; + +typedef int Vertex_handle; +typedef double Filtration_value; + +void program_options( int argc, char * argv[] + , std::string & filepoints + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence ); + +int main (int argc, char * argv[]) +{ + std::string filepoints; + std::string filediag ; + Filtration_value threshold ; + int dim_max ; + int p ; + Filtration_value min_persistence; + + program_options(argc,argv,filepoints,filediag,threshold,dim_max,p,min_persistence); + +// Extract the points from the file filepoints + typedef std::vector<double> Point_t; + std::vector< Point_t > points; + read_points( filepoints, points ); + +// Compute the proximity graph of the points + Graph_t prox_graph = compute_proximity_graph( points, threshold + , euclidean_distance<Point_t> ); + +// Construct the Rips complex in a Simplex Tree + Simplex_tree<> st; + st.insert_graph(prox_graph); // insert the proximity graph in the simplex tree + st.expansion( dim_max ); // expand the graph until dimension 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< Simplex_tree<>, Field_Zp > pcoh( st ); + pcoh.init_coefficients( p ); //initilizes the coefficient field for homology + + 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 & filepoints + , 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>(&filepoints), + "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"); + + 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; + std::abort(); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp b/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp new file mode 100644 index 00000000..ddaafea2 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp @@ -0,0 +1,140 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 INRIA Saclay (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "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> + +using namespace Gudhi; + +typedef int Vertex_handle; +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; + + // Construct the Rips complex in a Simplex Tree + 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() << " - filtration " << simplex_tree.filtration() << 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 ); + pcoh.init_coefficients( p ); //initilizes the coefficient field for homology + + 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; + std::abort(); + } +} 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..cd4a2bcc --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -0,0 +1,772 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef _PERSISTENCECOMPUTATION_SIMPLEXTREE_ +#define _PERSISTENCECOMPUTATION_SIMPLEXTREE_ + +#include <boost/tuple/tuple.hpp> +#include <boost/intrusive/set.hpp> +#include <boost/pending/disjoint_sets.hpp> +#include <boost/intrusive/list.hpp> +#include <boost/pool/object_pool.hpp> +#include "gudhi/Persistent_cohomology/Persistent_cohomology_column.h" +#include "gudhi/Persistent_cohomology/Field_Zp.h" + +namespace Gudhi{ + +/** \defgroup persistent_cohomology Persistent Cohomology Package + * + * Computation of persistent cohomology using the algorithm of + * \cite DBLP:journals/dcg/SilvaMV11 and \cite DBLP:journals/corr/abs-1208-5018 + * and the Compressed Annotation Matrix + * implementation of \cite DBLP:conf/esa/BoissonnatDM13 + * + * + * + + The theory of homology consists in attaching to a topological space a sequence of + (homology) groups, + capturing global topological features + like connected components, holes, cavities, etc. Persistent homology studies the evolution + -- birth, life and death -- of + these features when the topological space is changing. Consequently, the theory is essentially + composed of three elements: + topological spaces, their homology groups and an evolution scheme. + + <DT>Topological Spaces:</DT> + Topological spaces are represented by simplicial complexes. + Let \f$V = \{1, \cdots ,|V|\}\f$ be a set of <EM>vertices</EM>. + A <EM>simplex</EM> \f$\sigma\f$ is a subset of vertices + \f$\sigma \subseteq V\f$. A <EM>simplicial complex</EM> \f$\mathbf{K}\f$ + on \f$V\f$ is a collection of simplices \f$\{\sigma\}\f$, + \f$\sigma \subseteq V\f$, such that \f$\tau \subseteq \sigma \in \mathbf{K} + \Rightarrow \tau \in \mathbf{K}\f$. The dimension \f$n=|\sigma|-1\f$ of \f$\sigma\f$ + is its number of elements minus 1. A <EM>filtration</EM> of a simplicial complex is + a function \f$f:\mathbf{K} \rightarrow \mathbb{R}\f$ satisfying \f$f(\tau)\leq + f(\sigma)\f$ whenever \f$\tau \subseteq \sigma\f$. + + We define the concept FilteredComplex which enumerates the requirements for a class + to represent a filtered complex from which persistent homology may be computed. + We use the vocabulary of simplicial complexes, but the concept + is valid for any type of cell complex. The main requirements + are the definition of: + \li type <CODE>Indexing_tag</CODE>, which is a model of the concept + <CODE>IndexingTag</CODE>, + describing the nature of the indexing scheme, + \li type Simplex_handle to manipulate simplices, + \li method <CODE>int dimension(Simplex_handle)</CODE> returning + the dimension of a simplex, + \li type and method <CODE>Boundary_simplex_range + boundary_simplex_range(Simplex_handle)</CODE> that returns + a range giving access to the codimension 1 subsimplices of the + input simplex, as-well-as the coefficients \f$(-1)^i\f$ in the + definition of the operator \f$\partial\f$. The iterators have + value type <CODE>Simplex_handle</CODE>, + \li type and method + <CODE>Filtration_simplex_range filtration_simplex_range ()</CODE> + that returns a range giving + access to all the simplices of the complex read in the order + assigned by the indexing scheme, + \li type and method + <CODE>Filtration_value filtration (Simplex_handle)</CODE> that returns the value of + the filtration on the simplex represented by the handle. + + <DT>Homology:</DT> + For a ring \f$\mathcal{R}\f$, the group of <EM>n-chains</EM>, + denoted \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$, of \f$\mathbf{K}\f$ is the + group of formal sums of + n-simplices with \f$\mathcal{R}\f$ coefficients. The <EM>boundary operator</EM> is a + linear operator + \f$\partial_n: \mathbf{C}_n(\mathbf{K},\mathcal{R}) \rightarrow \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R})\f$ + such that \f$\partial_n \sigma = \partial_n [v_0, \cdots , v_n] = + \sum_{i=0}^n (-1)^{i}[v_0,\cdots ,\widehat{v_i}, \cdots,v_n]\f$, + where \f$\widehat{v_i}\f$ means \f$v_i\f$ is omitted from the list. The chain + groups form a sequence: + + \f[\cdots \ \ \mathbf{C}_n(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_n\ } \mathbf{C}_{n-1}(\mathbf{K},\mathcal{R}) + \xrightarrow{\partial_{n-1}} \cdots \xrightarrow{\ \partial_2 \ } + \mathbf{C}_1(\mathbf{K},\mathcal{R}) \xrightarrow{\ \partial_1 \ } \mathbf{C}_0(\mathbf{K},\mathcal{R}) \f] + + of finitely many groups \f$\mathbf{C}_n(\mathbf{K},\mathcal{R})\f$ and homomorphisms + \f$\partial_n\f$, indexed by the dimension \f$n \geq 0\f$. + The boundary operators satisfy the property \f$\partial_n \circ \partial_{n+1}=0\f$ + for every \f$n > 0\f$ + and we define the homology groups: + + \f[\mathbf{H}_n(\mathbf{K},\mathcal{R}) = \ker \partial_n / \mathrm{im} \ \partial_{n+1}\f] + + We refer to \cite Munkres-elementsalgtop1984 for an introduction to homology + theory and to \cite DBLP:books/daglib/0025666 for an introduction to persistent homology. + + <DT>Indexing Scheme:</DT> + "Changing" a simplicial complex consists in applying a simplicial map. + An <EM>indexing scheme</EM> is a directed graph together with a traversal + order, such that two + consecutive nodes in the graph are connected by an arrow (either forward or backward). + The nodes represent simplicial complexes and the directed edges simplicial maps. + + From the computational point of view, there are two types of indexing schemes of + interest + in persistent homology: <EM>linear</EM> ones + \f$\bullet \longrightarrow \bullet \longrightarrow \cdots \longrightarrow \bullet + \longrightarrow \bullet\f$ + in persistent homology \cite DBLP:journals/dcg/ZomorodianC05 , + and <EM>zigzag</EM> ones + \f$\bullet \longrightarrow \bullet \longleftarrow \cdots + \longrightarrow \bullet + \longleftarrow \bullet \f$ in zigzag persistent + homology \cite DBLP:journals/focm/CarlssonS10. + These indexing schemes have a natural left-to-right traversal order, and we + describe them with ranges and iterators. + In the current release of the Gudhi library, only the linear case is implemented. + + In the following, we consider the case where the indexing scheme is induced + by a filtration. + Ordering the simplices + by increasing filtration values (breaking ties so as a simplex appears after + its subsimplices of same filtration value) provides an indexing scheme. + + + + <DT>Implementations:</DT> + We use the <EM>Compressed Annotation Matrix</EM> of \cite DBLP:conf/esa/BoissonnatDM13 to + implement the + persistent cohomology algorithm of \cite DBLP:journals/dcg/SilvaMV11 + and \cite DBLP:conf/compgeom/DeyFW14 for persistence in the class Persistent_cohomology. + + The coefficient fields available as models of CoefficientField are Field_Zp + for \f$\mathbb{Z}_p\f$ (for any prime p) and Multi_field for the multi-field persistence algorithm + -- computing persistence simultaneously in various coefficient fields -- described + in \cite boissonnat:hal-00922572. + + + \author Clément Maria + \version 1.0 + \date 2014 + \copyright GNU General Public License v3. + @{ + */ + +/** \brief Computes the persistent cohomology of a filtered complex. +* +* 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 +* +*/ +//Memory allocation policy: classic, use a mempool, etc.*/ + template < class FilteredComplex + , class CoefficientField + > //to do mem allocation policy: classic, mempool, etc. + class Persistent_cohomology { + public: + + typedef FilteredComplex Complex_ds; + // Data attached to each simplex to interface with a Property Map. + typedef typename Complex_ds::Simplex_key Simplex_key; + typedef typename Complex_ds::Simplex_handle Simplex_handle; + typedef typename Complex_ds::Filtration_value Filtration_value; + typedef typename CoefficientField::Element Arith_element; +// Compressed Annotation Matrix types: + // 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; +// Persistent interval type. The Arith_element field is used for the multi-field framework. + typedef boost::tuple< Simplex_handle + , Simplex_handle + , Arith_element > Persistent_interval; + +/** \brief Initializes the Persistent_cohomology class. + * + * @param[in] cpx Complex for which the persistent homology is compiuted. + 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. + */ +Persistent_cohomology ( Complex_ds & 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. +, ds_rank_ (cpx_->num_simplices()) // union-find +, ds_parent_(cpx_->num_simplices()) // union-find +, ds_repr_ (cpx_->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_(new boost::object_pool< Column > ()) // memory pools for the CAM +, cell_pool_(new boost::object_pool< Cell > ()) +{ + if( persistence_dim_max ) { ++dim_max_; } + 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)); + } +} + +~Persistent_cohomology() +{ +//Clean the remaining columns in the matrix. + for(auto & cam_ref : cam_) { cam_ref.col_.clear(); } +//Clean the transversal lists + for( auto & transverse_ref : transverse_idx_ ) + { transverse_ref.second.row_->clear(); delete transverse_ref.second.row_; } +//Clear the memory pools + delete column_pool_; + delete cell_pool_; +} + +private: +struct length_interval { + length_interval ( Complex_ds * 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; } + + Complex_ds * 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 disgards 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_.push_back( Persistent_interval ( cpx_->simplex(key) + , cpx_->null_simplex() + , coeff_field_.characteristic() ) + ); + } + } + for( auto zero_idx : zero_cocycles_ ) + { + persistent_pairs_.push_back( Persistent_interval ( cpx_->simplex(zero_idx.second) + , cpx_->null_simplex() + , coeff_field_.characteristic() ) + ); + } +// Compute infinite interval of dimension > 0 + for(auto cocycle : transverse_idx_) + { + persistent_pairs_.push_back( Persistent_interval ( 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_.push_back ( Persistent_interval ( 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_.push_back ( Persistent_interval ( 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 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, in a map. + std::map < Column *, int > annotations_in_boundary; + std::pair < typename std::map< Column *, int >::iterator + , bool > result_insert_bound; + 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 + { // vector. + // 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". + result_insert_bound = + annotations_in_boundary.insert(std::pair<Column *,int>(curr_col,sign)); + if( !(result_insert_bound.second) ) { result_insert_bound.first->second += sign; } + } + } + sign = -sign; + } + // 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_ref : annotations_in_boundary ) + { + if(ann_ref.second != coeff_field_.additive_identity()) // For all columns in the boundary, + { + for( auto cell_ref : ann_ref.first->col_ ) // insert every cell in map_a_ds with multiplicity + { + Arith_element w_y = + coeff_field_.times(cell_ref.coefficient_ , ann_ref.second); //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 + { + 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(Column(key)); + Cell * new_cell = cell_pool_->construct(Cell (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 + if(interval_length_policy(cpx_->simplex(death_key),sigma)) { + persistent_pairs_.push_back ( Persistent_interval ( cpx_->simplex(death_key) //creator + , sigma //destructor + , charac ) //fields + ); // for which the interval exists + } + + 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_it = curr_col->col_.begin(); + col_cell_it != curr_col->col_.end(); ++col_cell_it ) + { col_cell_it->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_->free(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_it = curr_col->col_.begin(); + col_cell_it != curr_col->col_.end(); ++col_cell_it ) //re-establish the row links + { transverse_idx_[ col_cell_it->key_ ].row_->push_front(*col_cell_it); } + } + 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; + column_pool_->free(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)); + + 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 + 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 + + coeff_field_.clear_coefficient(tmp_cell_ptr->coefficient_); + cell_pool_->free(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 //key + , coeff_field_.additive_identity() + , &target)); + + 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 { + cmp_intervals_by_length( Complex_ds * sc ) : sc_ (sc) {} + bool operator() ( Persistent_interval & p1 + , Persistent_interval & p2 ) + { + return ( sc_->filtration( get<1>(p1) ) - sc_->filtration( get<0>(p1) ) + > sc_->filtration( get<1>(p2) ) - sc_->filtration( get<0>(p2) ) ); + } + Complex_ds * 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_ ); + persistent_pairs_.sort( cmp ); + for(auto pair : persistent_pairs_) + { + ostream << get<2>(pair) << " " + << cpx_->dimension(get<0>(pair)) << " " + << cpx_->filtration(get<0>(pair)) << " " + << cpx_->filtration(get<1>(pair)) << " " + << std::endl; + } +} + +private: +/* + * Structure representing a cocycle. + */ +struct cocycle { + cocycle() {} + 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: + Complex_ds * cpx_; + int dim_max_; + CoefficientField coeff_field_; + +/* 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::list< Persistent_interval > persistent_pairs_; + length_interval interval_length_policy; + + boost::object_pool< Column > * column_pool_; + boost::object_pool< Cell > * cell_pool_; +}; + +/** @} */ //end defgroup persistent_cohomology + +} // namespace GUDHI + +#endif // _PERSISTENCECOMPUTATION_SIMPLEXTREE_ 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..af0d6605 --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -0,0 +1,106 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_FIELD_ZP_H +#define GUDHI_FIELD_ZP_H + +namespace Gudhi{ + +/** \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(-1) +, inverse_() {} + +void init(int charac ) { + assert(charac <= 32768); + 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*/ +void plus_times_equal ( Element & x, Element y, Element w ) +{ x = (x + w * y) % Prime; } + +// operator= defined on Element + +/** Returns y * w */ +Element times ( Element y, int w ) { + Element res = (y * w) % Prime; + if(res < 0) return res+Prime; + else return res; +} + +void clear_coefficient(Element x) {} + +void plus_equal(Element & x, Element y) { x = ((x+y)%Prime); } + +/** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ +Element additive_identity () { return 0; } +/** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ +Element multiplicative_identity ( Element P = 0) { 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 ) +{ + Element out = (-x * y) % Prime; + return (out < 0) ? out + Prime : out; +} + + +bool is_one ( Element x ) { return x == 1; } +bool is_zero ( Element x ) { return x == 0; } + +//bool is_null() + +/** \brief Returns the characteristic \f$p\f$ of the field.*/ +Element characteristic() { return Prime; } + +private: + Element Prime; +/** Property map Element -> Element, which associate to an element its inverse in the field.*/ + std::vector< Element > inverse_; +}; + +} // namespace GUDHI + +#endif // GUDHI_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..9dd0998c --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -0,0 +1,167 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_MULTI_FIELD_H +#define GUDHI_MULTI_FIELD_H + +#include <iostream> +#include <vector> +#include <gmpxx.h> + +namespace Gudhi{ + +/** \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 () {} + +/* 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 + unsigned 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); + } + //set m to primorial(bound_prime) + prod_characteristics_ = 1; + for(auto p : primes_) + { mpz_mul_ui(prod_characteristics_.get_mpz_t(), + prod_characteristics_.get_mpz_t(), + p); + } + + num_primes_ = primes_.size(); + + //Uvect_ + Element Ui; Element tmp_elem; + for(auto p : primes_) + { + 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(int idx = 0; idx < num_primes_; ++idx) + { mult_id_all = (mult_id_all + Uvect_[idx]) % prod_characteristics_; } + + } + + void clear_coefficient(Element & x) { mpz_clear(x.get_mpz_t()); } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + Element additive_identity () { return 0; } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + Element multiplicative_identity () { return mult_id_all; }// 1 everywhere + + Element multiplicative_identity (Element Q) + { + if(Q == prod_characteristics_) { return multiplicative_identity(); } + + Element mult_id = 0; + for(int idx = 0; idx < num_primes_; ++idx) { + if( (Q % primes_[idx]) == 0 ) + { mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; } + } + return mult_id; + } + + /** Returns y * w */ + Element times ( Element y, int w ) { + Element tmp = (y*w) % prod_characteristics_; + if(tmp < 0) return prod_characteristics_ + tmp; + return tmp; + } + + void plus_equal(Element & x, Element y) + { x += y; x %= prod_characteristics_; } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + Element characteristic() { 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()); + + return std::pair<Element,Element>( + (inv_qt * multiplicative_identity(QT)) % prod_characteristics_ + , QT ); + } + /** Returns -x * y.*/ + Element times_minus ( Element x, Element y ) + { return prod_characteristics_ - ((x*y)%prod_characteristics_); } + + /** Set x <- x + w * y*/ + void plus_times_equal ( Element & x, Element y, Element w ) + { x = (x + w * y) % prod_characteristics_; } + + 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_; + size_t num_primes_; //number of fields + Element mult_id_all; + +}; + +} // namespace GUDHI + +#endif // GUDHI_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..0702c58e --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h @@ -0,0 +1,153 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_COLUMN_LIST_H +#define GUDHI_COLUMN_LIST_H + +#include "boost/tuple/tuple.hpp" +#include "boost/intrusive/set.hpp" +#include "boost/intrusive/list.hpp" + +namespace Gudhi{ + +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. + */ +template < typename SimplexKey + , typename ArithmeticElement > +class Persistent_cohomology_column +: public boost::intrusive::set_base_hook + < boost::intrusive::link_mode< boost::intrusive::normal_link > > +{ +private: + template < class T1, class T2 > friend class Persistent_cohomology; + + 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.*/ + Persistent_cohomology_column (SimplexKey key) + { + class_key_ = key; + col_ = Col_type(); + } +public: + /** Copy constructor.*/ + Persistent_cohomology_column( Persistent_cohomology_column const &other ) + : col_() + , class_key_(other.class_key_) + { if(!other.col_.empty()) std::cerr << "Copying a non-empty column.\n"; } + +/** \brief Returns true iff the column is null.*/ + bool is_null() { 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 () { 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()); + } + + // void display() + // { + // for(auto cell : col_) + // { std::cout << "(" << cell.key_ <<":"<<cell.coefficient_<<") "; } + // } + + Col_type col_; + SimplexKey class_key_; +}; + +} // namespace GUDHI + +#endif // GUDHI_COLUMN_LIST_H |