diff options
Diffstat (limited to 'src/Persistent_cohomology/example')
9 files changed, 1100 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..94ec13c5 --- /dev/null +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -0,0 +1,68 @@ +project(Persistent_cohomology_examples) + +add_executable(plain_homology plain_homology.cpp) + +add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp) + +add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp) +target_link_libraries(rips_persistence_step_by_step ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(rips_persistence_via_boundary_matrix rips_persistence_via_boundary_matrix.cpp) +target_link_libraries(rips_persistence_via_boundary_matrix ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(persistence_from_file persistence_from_file.cpp) +target_link_libraries(persistence_from_file ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(plain_homology ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_step_by_step ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_file ${TBB_LIBRARIES}) +endif() + +add_test(NAME Persistent_cohomology_example_plain_homology COMMAND $<TARGET_FILE:plain_homology>) +add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $<TARGET_FILE:persistence_from_simple_simplex_tree> + "1" "0") +add_test(NAME Persistent_cohomology_example_from_rips_step_by_step_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence_step_by_step> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +add_test(NAME Persistent_cohomology_example_via_boundary_matrix COMMAND $<TARGET_FILE:rips_persistence_via_boundary_matrix> + "${CMAKE_SOURCE_DIR}/data/points/Kl.off" "-r" "0.16" "-d" "3" "-p" "3" "-m" "100") +add_test(NAME Persistent_cohomology_example_from_file_3_2_0 COMMAND $<TARGET_FILE:persistence_from_file> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc" "-p" "2" "-m" "0") +add_test(NAME Persistent_cohomology_example_from_file_3_3_100 COMMAND $<TARGET_FILE:persistence_from_file> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc" "-p" "3" "-m" "100") + +install(TARGETS plain_homology DESTINATION bin) +install(TARGETS persistence_from_simple_simplex_tree DESTINATION bin) +install(TARGETS rips_persistence_step_by_step DESTINATION bin) +install(TARGETS rips_persistence_via_boundary_matrix DESTINATION bin) +install(TARGETS persistence_from_file DESTINATION bin) + +if(GMP_FOUND) + if(GMPXX_FOUND) + add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) + target_link_libraries(rips_multifield_persistence + ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistent_cohomology_example_multifield_2_71 COMMAND $<TARGET_FILE:rips_multifield_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "2" "-q" "71") + install(TARGETS rips_multifield_persistence DESTINATION bin) + endif(GMPXX_FOUND) +endif(GMP_FOUND) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + + add_executable(custom_persistence_sort custom_persistence_sort.cpp) + target_link_libraries(custom_persistence_sort ${CGAL_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistent_cohomology_example_custom_persistence_sort COMMAND $<TARGET_FILE:custom_persistence_sort>) + + install(TARGETS custom_persistence_sort DESTINATION bin) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Persistent_cohomology/example/README b/src/Persistent_cohomology/example/README new file mode 100644 index 00000000..f39d9584 --- /dev/null +++ b/src/Persistent_cohomology/example/README @@ -0,0 +1,67 @@ +To build the examples, run in a Terminal: + +cd /path-to-examples/ +cmake . +make + +*********************************************************************************************************************** +Example of use of RIPS: + +Computation of the persistent homology with Z/2Z and Z/3Z coefficients simultaneously of the Rips complex +on points sampling a 3D torus: + +./rips_multifield_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.12 -d 3 -p 2 -q 3 + +output: +6 0 0 inf +6 1 0.0983494 inf +6 1 0.104347 inf +6 2 0.138335 inf +6 0 0 0.122545 +6 0 0 0.121171 +6 0 0 0.120964 +6 0 0 0.12057 +6 0 0 0.12047 +6 0 0 0.120414 + +Every line is of this format: p1*...*pr dim b d +where + p1*...*pr is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients. + dim is the dimension of the homological feature, + b and d are respectively the birth and death of the feature and + +and the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers): + + ./rips_multifield_persistence ../../data/points/Kl.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71 + +output: +557940830126698960967415390 0 0 inf +557940830126698960967415390 1 0.0983494 inf +557940830126698960967415390 1 0.104347 inf +557940830126698960967415390 2 0.138335 inf +557940830126698960967415390 0 0 0.122545 +557940830126698960967415390 0 0 0.121171 +557940830126698960967415390 0 0 0.120964 +557940830126698960967415390 0 0 0.12057 +557940830126698960967415390 0 0 0.12047 +557940830126698960967415390 0 0 0.120414 + +*********************************************************************************************************************** +Example of use of PLAIN HOMOLOGY: + +This example computes the plain homology of the following simplicial complex without filtration values: + /* Complex to build. */ + /* 1 3 */ + /* o---o */ + /* /X\ / */ + /* o---o o */ + /* 2 0 4 */ + +./plain_homology + +output: +2 0 0 inf +2 0 0 inf +2 1 0 inf + +Here we retrieve the 2 entities {0,1,2,3} and {4} (Betti numbers[0] = 2) and the hole in {0,1,3} (Betti numbers[1] = 1) diff --git a/src/Persistent_cohomology/example/custom_persistence_sort.cpp b/src/Persistent_cohomology/example/custom_persistence_sort.cpp new file mode 100644 index 00000000..be74cf50 --- /dev/null +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -0,0 +1,125 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Persistent_cohomology.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <fstream> // for std::ofstream +#include <algorithm> // for std::sort + + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Point = Kernel::Point_d; +using Alpha_complex = Gudhi::alpha_complex::Alpha_complex<Kernel>; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, + Gudhi::persistent_cohomology::Field_Zp >; + +std::vector<Point> random_points() { + // Instanciate a random point generator + CGAL::Random rng(0); + + // Generate "points_number" random points in a vector + std::vector<Point> points; + + // Generates 1000 random 3D points on a sphere of radius 4.0 + CGAL::Random_points_on_sphere_d<Point> rand_outside(3, 4.0, rng); + CGAL::cpp11::copy_n(rand_outside, 1000, std::back_inserter(points)); + // Generates 2000 random 3D points in a sphere of radius 3.0 + CGAL::Random_points_in_ball_d<Point> rand_inside(3, 3.0, rng); + CGAL::cpp11::copy_n(rand_inside, 2000, std::back_inserter(points)); + + return points; +} + +/* + * Compare two intervals by dimension, then by length. + */ +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Simplex_tree * sc) + : sc_(sc) { } + + template<typename Persistent_interval> + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2))) + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + else + return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); + } + Simplex_tree* sc_; +}; + +int main(int argc, char **argv) { + std::vector<Point> points = random_points(); + + std::cout << "Points size=" << points.size() << std::endl; + // Alpha complex persistence computation from generated points + Alpha_complex alpha_complex_from_points(points); + std::cout << "alpha_complex_from_points" << std::endl; + + Simplex_tree simplex; + std::cout << "simplex" << std::endl; + if (alpha_complex_from_points.create_complex(simplex, 0.6)) { + std::cout << "simplex" << std::endl; + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Simplicial complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + // Sort the simplices in the order of the filtration + simplex.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; + + Persistent_cohomology pcoh(simplex); + + // initializes the coefficient field for homology - Z/3Z + pcoh.init_coefficients(3); + pcoh.compute_persistent_cohomology(0.2); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(&simplex); + auto persistent_pairs = pcoh.get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + for (auto pair : persistent_pairs) { + std::cout << simplex.dimension(get<0>(pair)) << " " + << simplex.filtration(get<0>(pair)) << " " + << simplex.filtration(get<1>(pair)) << std::endl; + } + + // Persistent Betti numbers + std::cout << "The persistent Betti numbers in interval [0.40, 0.41] are : "; + for (int dim = 0; dim < simplex.dimension(); dim++) + std::cout << "b" << dim << " = " << pcoh.persistent_betti_number(dim, 0.40, 0.41) << " ; "; + std::cout << std::endl; + + // Betti numbers + std::vector<int> betti_numbers = pcoh.betti_numbers(); + std::cout << "The Betti numbers are : "; + for (std::size_t i = 0; i < betti_numbers.size(); i++) + std::cout << "b" << i << " = " << betti_numbers[i] << " ; "; + std::cout << std::endl; + } + return 0; +} + diff --git a/src/Persistent_cohomology/example/persistence_from_file.cpp b/src/Persistent_cohomology/example/persistence_from_file.cpp new file mode 100644 index 00000000..d169cc63 --- /dev/null +++ b/src/Persistent_cohomology/example/persistence_from_file.cpp @@ -0,0 +1,130 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/reader_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <boost/program_options.hpp> + +#include <string> + +using namespace Gudhi; +using namespace Gudhi::persistent_cohomology; + +typedef double Filtration_value; + +void program_options(int argc, char * argv[] + , std::string & simplex_tree_file + , std::string & output_file + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string simplex_tree_file; + std::string output_file; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, simplex_tree_file, output_file, p, min_persistence); + + std::cout << "Simplex_tree from file=" << simplex_tree_file.c_str() << " - output_file=" << output_file.c_str() + << std::endl; + std::cout << " - p=" << p << " - min_persistence=" << min_persistence << std::endl; + + // Read the list of simplices from a file. + Simplex_tree<> simplex_tree; + + std::ifstream simplex_tree_stream(simplex_tree_file); + simplex_tree_stream >> simplex_tree; + + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices" << std::endl; + std::cout << " - dimension " << simplex_tree.dimension() << std::endl; + + /* + std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for( auto f_simplex : simplex_tree.filtration_simplex_range() ) + { std::cout << " " << "[" << simplex_tree.filtration(f_simplex) << "] "; + for( auto vertex : simplex_tree.simplex_vertex_range(f_simplex) ) + { std::cout << vertex << " "; } + std::cout << std::endl; + }*/ + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in output_file + if (output_file.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(output_file); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & simplex_tree_file + , std::string & output_file + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&simplex_tree_file), + "Name of file containing a simplex set. Format is one simplex per line (cf. reader_utils.h - read_simplex): Dim1 X11 X12 ... X1d Fil1 "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&output_file)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp new file mode 100644 index 00000000..3c91662f --- /dev/null +++ b/src/Persistent_cohomology/example/persistence_from_simple_simplex_tree.cpp @@ -0,0 +1,163 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <ctime> +#include <utility> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using typeVectorVertex = std::vector< Simplex_tree::Vertex_handle >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; + exit(-1); +} + +int main(int argc, char * const argv[]) { + // program args management + if (argc != 3) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = 0; + int returnedScanValue = sscanf(argv[1], "%d", &coeff_field_characteristic); + if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) { + std::cerr << "Error: " << argv[1] << " is not correct\n"; + usage(argv[0]); + } + + Filtration_value min_persistence = 0.0; + returnedScanValue = sscanf(argv[2], "%lf", &min_persistence); + if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { + std::cerr << "Error: " << argv[2] << " is not correct\n"; + usage(argv[0]); + } + + // TEST OF INSERTION + std::cout << "********************************************************************" << std::endl; + std::cout << "TEST OF INSERTION" << std::endl; + Simplex_tree st; + + // ++ FIRST + std::cout << " - INSERT (0,1,2)" << std::endl; + typeVectorVertex SimplexVector = {0, 1, 2}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ SECOND + std::cout << " - INSERT 3" << std::endl; + SimplexVector = {3}; + st.insert_simplex_and_subfaces(SimplexVector, 0.1); + + // ++ THIRD + std::cout << " - INSERT (0,3)" << std::endl; + SimplexVector = {0, 3}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ FOURTH + std::cout << " - INSERT (0,1) (already inserted)" << std::endl; + SimplexVector = {0, 1}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ FIFTH + std::cout << " - INSERT (3,4,5)" << std::endl; + SimplexVector = {3, 4, 5}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ SIXTH + std::cout << " - INSERT (0,1,6,7)" << std::endl; + SimplexVector = {0, 1, 6, 7}; + st.insert_simplex_and_subfaces(SimplexVector, 0.4); + + // ++ SEVENTH + std::cout << " - INSERT (4,5,8,9)" << std::endl; + SimplexVector = {4, 5, 8, 9}; + st.insert_simplex_and_subfaces(SimplexVector, 0.4); + + // ++ EIGHTH + std::cout << " - INSERT (9,10,11)" << std::endl; + SimplexVector = {9, 10, 11}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ NINETH + std::cout << " - INSERT (2,10,12)" << std::endl; + SimplexVector = {2, 10, 12}; + st.insert_simplex_and_subfaces(SimplexVector, 0.3); + + // ++ TENTH + std::cout << " - INSERT (11,6)" << std::endl; + SimplexVector = {6, 11}; + st.insert_simplex_and_subfaces(SimplexVector, 0.2); + + // ++ ELEVENTH + std::cout << " - INSERT (13,14,15)" << std::endl; + SimplexVector = {13, 14, 15}; + st.insert_simplex_and_subfaces(SimplexVector, 0.25); + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ 4 2 */ + /* o---o---o---o o */ + /* 2 0 3\X/8\ 10 /X\ */ + /* o---o---o---o */ + /* 5 9\X/ 12 */ + /* o---o */ + /* 11 6 */ + /* In other words: */ + /* A facet [2,1,0] */ + /* An edge [0,3] */ + /* A facet [3,4,5] */ + /* A cell [0,1,6,7] */ + /* A cell [4,5,8,9] */ + /* A facet [9,10,11] */ + /* An edge [11,6] */ + /* An edge [10,12,2] */ + + + std::cout << "The complex contains " << st.num_simplices() << " simplices - " << st.num_vertices() << " vertices " + << std::endl; + std::cout << " - dimension " << st.dimension() << std::endl; + std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" + << std::endl; + std::cout << "**************************************************************" << std::endl; + std::cout << "strict graph G { " << std::endl; + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << " " << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << static_cast<int>(vertex) << " -- "; + } + std::cout << ";" << std::endl; + } + + std::cout << "}" << std::endl; + std::cout << "**************************************************************" << std::endl; + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + pcoh.output_diagram(); + return 0; +} diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp new file mode 100644 index 00000000..84333e46 --- /dev/null +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -0,0 +1,91 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <vector> +#include <cstdint> // for std::uint8_t + +/* We could perfectly well use the default Simplex_tree<> (which uses + * Simplex_tree_options_full_featured), the following simply demonstrates + * how to save on storage by not storing a filtration value. */ + +struct MyOptions : Gudhi::Simplex_tree_options_full_featured { + // Implicitly use 0 as filtration value for all simplices + static const bool store_filtration = false; + // The persistence algorithm needs this + static const bool store_key = true; + // I have few vertices + typedef short Vertex_handle; + // Maximum number of simplices to compute persistence is 2^8 - 1 = 255. One is reserved for null_key + typedef std::uint8_t Simplex_key; +}; + +using ST = Gudhi::Simplex_tree<MyOptions>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<ST, Field_Zp>; + +int main() { + ST st; + + /* Complex to build. */ + /* 1 3 5 */ + /* o---o---o */ + /* / \ / */ + /* o---o o */ + /* 2 0 4 */ + + const short edge01[] = {0, 1}; + const short edge02[] = {0, 2}; + const short edge12[] = {1, 2}; + const short edge03[] = {0, 3}; + const short edge13[] = {1, 3}; + const short edge35[] = {3, 5}; + const short vertex4[] = {4}; + st.insert_simplex_and_subfaces(edge01); + st.insert_simplex_and_subfaces(edge02); + st.insert_simplex_and_subfaces(edge12); + st.insert_simplex_and_subfaces(edge03); + st.insert_simplex(edge13); + st.insert_simplex_and_subfaces(edge35); + st.insert_simplex(vertex4); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Class for homology computation + // By default, since the complex has dimension 1, only 0-dimensional homology would be computed. + // Here we also want persistent homology to be computed for the maximal dimension in the complex (persistence_dim_max = true) + Persistent_cohomology pcoh(st, true); + + // Initialize the coefficient field Z/2Z for homology + pcoh.init_coefficients(2); + + // Compute the persistence diagram of the complex + pcoh.compute_persistent_cohomology(); + + // Print the result. The format is, on each line: 2 dim 0 inf + // where 2 represents the field, dim the dimension of the feature. + // 2 0 0 inf + // 2 0 0 inf + // 2 1 0 inf + // 2 1 0 inf + // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=2. + pcoh.output_diagram(); + + // Print the Betti numbers are b0=2 and b1=2. + std::cout << std::endl; + std::cout << "The Betti numbers are : "; + for (int i = 0; i < 3; i++) + std::cout << "b" << i << " = " << pcoh.betti_number(i) << " ; "; + std::cout << std::endl; +} diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp new file mode 100644 index 00000000..9eb5ccfc --- /dev/null +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -0,0 +1,142 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Persistent_cohomology/Multi_field.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Multi_field = Gudhi::persistent_cohomology::Multi_field; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Multi_field >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & min_p + , int & max_p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int min_p; + int max_p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(min_p, max_p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & min_p + , int & max_p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options"); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("min-field-charac,p", po::value<int>(&min_p)->default_value(2), + "Minimal characteristic p of the coefficient field Z/pZ.") + ("max-field-charac,q", po::value<int>(&max_p)->default_value(1223), + "Minimial characteristic q of the coefficient field Z/pZ.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with various coefficient fields \n"; + std::cout << "of a Rips complex defined on a set of input points. The coefficient \n"; + std::cout << "fields are all the Z/rZ for a prime number r contained in the \n"; + std::cout << "specified range [p,q]\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p1*...*pr dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p1*...*pr is the product of prime numbers pi such that the homology \n"; + std::cout << "feature exists in homology with Z/piZ coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp new file mode 100644 index 00000000..02db05ec --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -0,0 +1,154 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <map> + +// ---------------------------------------------------------------------------- +// rips_persistence_step_by_step is an example of each step that is required to +// build a Rips over a Simplex_tree. Please refer to rips_persistence to see +// how to do the same thing with the Rips_complex wrapper for less detailed +// steps. +// ---------------------------------------------------------------------------- + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + // insert the proximity graph in the simplex tree + st.insert_graph(prox_graph); + // expand the graph until dimension dim_max + st.expansion(dim_max); + + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << st.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp new file mode 100644 index 00000000..37fa5e93 --- /dev/null +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -0,0 +1,160 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clément Maria, Marc Glisse + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> + +#include <boost/program_options.hpp> + +#ifdef GUDHI_USE_TBB +#include <tbb/task_scheduler_init.h> +#endif + +#include <string> +#include <vector> + +//////////////////////////////////////////////////////////////// +// // +// WARNING: persistence computation itself is not parallel, // +// and this uses more memory than rips_persistence. // +// // +//////////////////////////////////////////////////////////////// + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree& st = *new Simplex_tree; + rips_complex_from_file.create_complex(st, dim_max); + + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << st.dimension() << " \n"; + +#ifdef GUDHI_USE_TBB + // Unnecessary, but clarifies which operations are parallel. + tbb::task_scheduler_init ts; +#endif + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + int count = 0; + for (auto sh : st.filtration_simplex_range()) + st.assign_key(sh, count++); + + // Convert to a more convenient representation. + Gudhi::Hasse_complex<> hcpx(st); + +#ifdef GUDHI_USE_TBB + ts.terminate(); +#endif + + // Free some space. + delete &st; + + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Hasse_complex<>, Field_Zp > pcoh(hcpx); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} |