summaryrefslogtreecommitdiff
path: root/src/Persistent_cohomology/example
diff options
context:
space:
mode:
Diffstat (limited to 'src/Persistent_cohomology/example')
-rw-r--r--src/Persistent_cohomology/example/CMakeLists.txt22
-rw-r--r--src/Persistent_cohomology/example/README47
-rw-r--r--src/Persistent_cohomology/example/performance_rips_persistence.cpp168
-rw-r--r--src/Persistent_cohomology/example/rips_multifield_persistence.cpp152
-rw-r--r--src/Persistent_cohomology/example/rips_persistence.cpp146
-rw-r--r--src/Persistent_cohomology/example/rips_persistence_from_alpha_shape_3.cpp140
6 files changed, 675 insertions, 0 deletions
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();
+ }
+}