diff options
Diffstat (limited to 'example')
69 files changed, 3204 insertions, 923 deletions
diff --git a/example/Alpha_complex/Alpha_complex_from_off.cpp b/example/Alpha_complex/Alpha_complex_from_off.cpp index 7836d59a..d411e90a 100644 --- a/example/Alpha_complex/Alpha_complex_from_off.cpp +++ b/example/Alpha_complex/Alpha_complex_from_off.cpp @@ -1,4 +1,7 @@ #include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + #include <CGAL/Epick_d.h> #include <iostream> @@ -14,14 +17,14 @@ void usage(int nbArgs, char * const progName) { int main(int argc, char **argv) { if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1)); - std::string off_file_name(argv[1]); - double alpha_square_max_value = atof(argv[2]); + std::string off_file_name {argv[1]}; + double alpha_square_max_value {atof(argv[2])}; // ---------------------------------------------------------------------------- // Init of an alpha complex from an OFF file // ---------------------------------------------------------------------------- - typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel; - Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name, alpha_square_max_value); + using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name); std::streambuf* streambufffer; std::ofstream ouput_file_stream; @@ -33,23 +36,27 @@ int main(int argc, char **argv) { streambufffer = std::cout.rdbuf(); } - std::ostream output_stream(streambufffer); - - // ---------------------------------------------------------------------------- - // Display information about the alpha complex - // ---------------------------------------------------------------------------- - output_stream << "Alpha complex is of dimension " << alpha_complex_from_file.dimension() << - " - " << alpha_complex_from_file.num_simplices() << " simplices - " << - alpha_complex_from_file.num_vertices() << " vertices." << std::endl; - - output_stream << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; - for (auto f_simplex : alpha_complex_from_file.filtration_simplex_range()) { - output_stream << " ( "; - for (auto vertex : alpha_complex_from_file.simplex_vertex_range(f_simplex)) { - output_stream << vertex << " "; + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + output_stream << "Alpha complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + output_stream << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + output_stream << std::endl; } - output_stream << ") -> " << "[" << alpha_complex_from_file.filtration(f_simplex) << "] "; - output_stream << std::endl; } ouput_file_stream.close(); return 0; diff --git a/example/Alpha_complex/Alpha_complex_from_points.cpp b/example/Alpha_complex/Alpha_complex_from_points.cpp index 49f77276..c19f7cc8 100644 --- a/example/Alpha_complex/Alpha_complex_from_points.cpp +++ b/example/Alpha_complex/Alpha_complex_from_points.cpp @@ -1,14 +1,17 @@ -#include <CGAL/Epick_d.h> #include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <CGAL/Epick_d.h> #include <iostream> #include <string> #include <vector> #include <limits> // for numeric limits -typedef CGAL::Epick_d< CGAL::Dimension_tag<2> > Kernel; -typedef Kernel::Point_d Point; -typedef std::vector<Point> Vector_of_points; +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >; +using Point = Kernel::Point_d; +using Vector_of_points = std::vector<Point>; void usage(int nbArgs, char * const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; @@ -21,7 +24,7 @@ int main(int argc, char **argv) { if ((argc != 1) && (argc != 2)) usage(argc, (argv[0] - 1)); // Delaunay complex if alpha_square_max_value is not given by the user. - double alpha_square_max_value = std::numeric_limits<double>::infinity(); + double alpha_square_max_value {std::numeric_limits<double>::infinity()}; if (argc == 2) alpha_square_max_value = atof(argv[1]); @@ -40,23 +43,26 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points // ---------------------------------------------------------------------------- - Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points, alpha_square_max_value); + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points); - // ---------------------------------------------------------------------------- - // Display information about the alpha complex - // ---------------------------------------------------------------------------- - std::cout << "Alpha complex is of dimension " << alpha_complex_from_points.dimension() << - " - " << alpha_complex_from_points.num_simplices() << " simplices - " << - alpha_complex_from_points.num_vertices() << " vertices." << std::endl; - - std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; - for (auto f_simplex : alpha_complex_from_points.filtration_simplex_range()) { - std::cout << " ( "; - for (auto vertex : alpha_complex_from_points.simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_points.create_complex(simplex, alpha_square_max_value)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Alpha complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + std::cout << std::endl; } - std::cout << ") -> " << "[" << alpha_complex_from_points.filtration(f_simplex) << "] "; - std::cout << std::endl; } return 0; } diff --git a/example/Alpha_complex/CMakeLists.txt b/example/Alpha_complex/CMakeLists.txt index 71a95d61..a4853d78 100644 --- a/example/Alpha_complex/CMakeLists.txt +++ b/example/Alpha_complex/CMakeLists.txt @@ -2,32 +2,31 @@ cmake_minimum_required(VERSION 2.6) project(Alpha_complex_examples) # need CGAL 4.7 -# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../.. -if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.7.0) - if (EIGEN3_FOUND) - add_executable ( alphapoints Alpha_complex_from_points.cpp ) - target_link_libraries(alphapoints ${Boost_SYSTEM_LIBRARY} ${Boost_THREAD_LIBRARY} ${CGAL_LIBRARY}) - add_executable ( alphaoffreader Alpha_complex_from_off.cpp ) - target_link_libraries(alphaoffreader ${Boost_SYSTEM_LIBRARY} ${Boost_THREAD_LIBRARY} ${CGAL_LIBRARY}) - if (TBB_FOUND) - target_link_libraries(alphapoints ${TBB_LIBRARIES}) - target_link_libraries(alphaoffreader ${TBB_LIBRARIES}) - endif() +# cmake -DCGAL_DIR=~/workspace/CGAL-4.7 .. +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + add_executable ( Alpha_complex_example_from_points Alpha_complex_from_points.cpp ) + target_link_libraries(Alpha_complex_example_from_points ${Boost_SYSTEM_LIBRARY} ${Boost_THREAD_LIBRARY} ${CGAL_LIBRARY}) + add_executable ( Alpha_complex_example_from_off Alpha_complex_from_off.cpp ) + target_link_libraries(Alpha_complex_example_from_off ${Boost_SYSTEM_LIBRARY} ${Boost_THREAD_LIBRARY} ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES}) + target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES}) + endif() - add_test(alphapoints ${CMAKE_CURRENT_BINARY_DIR}/alphapoints) - # Do not forget to copy test files in current binary dir - file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) - add_test(alphaoffreader_doc_60 ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader alphacomplexdoc.off 60.0 ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt) - add_test(alphaoffreader_doc_32 ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader alphacomplexdoc.off 32.0 ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt) - if (DIFF_PATH) - # Do not forget to copy test results files in current binary dir - file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) - file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>) - add_test(alphaoffreader_doc_60_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) - add_test(alphaoffreader_doc_32_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) - endif() - endif(EIGEN3_FOUND) - endif(NOT CGAL_VERSION VERSION_LESS 4.7.0) -endif(CGAL_FOUND) + add_test(NAME Alpha_complex_example_from_off_60 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt") + add_test(NAME Alpha_complex_example_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt") + if (DIFF_PATH) + # Do not forget to copy test results files in current binary dir + file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(Alpha_complex_example_from_off_60_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) + add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) + endif() +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) diff --git a/example/Bitmap_cubical_complex/CMakeLists.txt b/example/Bitmap_cubical_complex/CMakeLists.txt index 2fddc514..241a11e5 100644 --- a/example/Bitmap_cubical_complex/CMakeLists.txt +++ b/example/Bitmap_cubical_complex/CMakeLists.txt @@ -6,21 +6,31 @@ target_link_libraries(Bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY}) if (TBB_FOUND) target_link_libraries(Bitmap_cubical_complex ${TBB_LIBRARIES}) endif() -add_test(Bitmap_cubical_complex_one_sphere ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt) -add_test(Bitmap_cubical_complex_two_sphere ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt) + +add_test(NAME Bitmap_cubical_complex_example_persistence_one_sphere COMMAND $<TARGET_FILE:Bitmap_cubical_complex> + "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt") + +add_test(NAME Bitmap_cubical_complex_example_persistence_two_sphere COMMAND $<TARGET_FILE:Bitmap_cubical_complex> + "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt") add_executable ( Random_bitmap_cubical_complex Random_bitmap_cubical_complex.cpp ) target_link_libraries(Random_bitmap_cubical_complex ${Boost_SYSTEM_LIBRARY}) if (TBB_FOUND) target_link_libraries(Random_bitmap_cubical_complex ${TBB_LIBRARIES}) endif() -add_test(Random_bitmap_cubical_complex ${CMAKE_CURRENT_BINARY_DIR}/Random_bitmap_cubical_complex 2 100 100) +add_test(NAME Bitmap_cubical_complex_example_random COMMAND $<TARGET_FILE:Random_bitmap_cubical_complex> + "2" "100" "100") add_executable ( Bitmap_cubical_complex_periodic_boundary_conditions Bitmap_cubical_complex_periodic_boundary_conditions.cpp ) target_link_libraries(Bitmap_cubical_complex_periodic_boundary_conditions ${Boost_SYSTEM_LIBRARY}) if (TBB_FOUND) target_link_libraries(Bitmap_cubical_complex_periodic_boundary_conditions ${TBB_LIBRARIES}) endif() -add_test(Bitmap_cubical_complex_periodic_2d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt) -add_test(Bitmap_cubical_complex_periodic_3d_torus ${CMAKE_CURRENT_BINARY_DIR}/Bitmap_cubical_complex_periodic_boundary_conditions ${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt) +add_test(NAME Bitmap_cubical_complex_example_periodic_boundary_conditions_2d_torus + COMMAND $<TARGET_FILE:Bitmap_cubical_complex_periodic_boundary_conditions> + "${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt") + +add_test(NAME Bitmap_cubical_complex_example_periodic_boundary_conditions_3d_torus + COMMAND $<TARGET_FILE:Bitmap_cubical_complex_periodic_boundary_conditions> + "${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt") diff --git a/example/Bottleneck_distance/CMakeLists.txt b/example/Bottleneck_distance/CMakeLists.txt new file mode 100644 index 00000000..0d0bff45 --- /dev/null +++ b/example/Bottleneck_distance/CMakeLists.txt @@ -0,0 +1,22 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_examples) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) + add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + + add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(bottleneck_read_file_example ${TBB_LIBRARIES}) + target_link_libraries(bottleneck_basic_example ${TBB_LIBRARIES}) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Bottleneck_distance_example_basic COMMAND $<TARGET_FILE:bottleneck_basic_example>) + + add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck + COMMAND $<TARGET_FILE:alpha_rips_persistence_bottleneck_distance> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3") + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/example/Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp b/example/Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp new file mode 100644 index 00000000..fd164b22 --- /dev/null +++ b/example/Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp @@ -0,0 +1,190 @@ +/* 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) 2017 INRIA + * + * 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/Alpha_complex.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/Bottleneck.h> + +#include <CGAL/Epick_d.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <algorithm> // for transform + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; +using Point_d = Kernel::Point_d; +using Points_off_reader = Gudhi::Points_off_reader<Point_d>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +static inline std::pair<double, double> compute_root_square(std::pair<double, double> input) { + return std::make_pair(std::sqrt(input.first), std::sqrt(input.second)); +} + +int main(int argc, char * argv[]) { + std::string off_file_points; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + + // -------------------------------------------- + // Rips persistence + // -------------------------------------------- + Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree rips_stree; + + rips_complex.create_complex(rips_stree, dim_max); + std::cout << "The Rips complex contains " << rips_stree.num_simplices() << " simplices and has dimension " + << rips_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + rips_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology rips_pcoh(rips_stree); + // initializes the coefficient field for homology + rips_pcoh.init_coefficients(p); + rips_pcoh.compute_persistent_cohomology(min_persistence); + + // rips_pcoh.output_diagram(); + + // -------------------------------------------- + // Alpha persistence + // -------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex(off_reader.get_point_cloud()); + + Simplex_tree alpha_stree; + alpha_complex.create_complex(alpha_stree, threshold * threshold); + std::cout << "The Alpha complex contains " << alpha_stree.num_simplices() << " simplices and has dimension " + << alpha_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + alpha_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology alpha_pcoh(alpha_stree); + // initializes the coefficient field for homology + alpha_pcoh.init_coefficients(p); + alpha_pcoh.compute_persistent_cohomology(min_persistence * min_persistence); + + // alpha_pcoh.output_diagram(); + + // -------------------------------------------- + // Bottleneck distance between both persistence + // -------------------------------------------- + double max_b_distance {}; + for (int dim = 0; dim < dim_max; dim ++) { + std::vector< std::pair< Filtration_value , Filtration_value > > rips_intervals; + std::vector< std::pair< Filtration_value , Filtration_value > > alpha_intervals; + rips_intervals = rips_pcoh.intervals_in_dimension(dim); + alpha_intervals = alpha_pcoh.intervals_in_dimension(dim); + std::transform(alpha_intervals.begin(), alpha_intervals.end(), alpha_intervals.begin(), compute_root_square); + + double bottleneck_distance = Gudhi::persistence_diagram::bottleneck_distance(rips_intervals, alpha_intervals); + std::cout << "In dimension " << dim << ", bottleneck distance = " << bottleneck_distance << std::endl; + if (bottleneck_distance > max_b_distance) + max_b_distance = bottleneck_distance; + } + std::cout << "================================================================================" << std::endl; + std::cout << "Bottleneck distance is " << max_b_distance << std::endl; + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , 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") + ("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; + std::abort(); + } +} diff --git a/example/Bottleneck_distance/bottleneck_basic_example.cpp b/example/Bottleneck_distance/bottleneck_basic_example.cpp new file mode 100644 index 00000000..d0ca4e20 --- /dev/null +++ b/example/Bottleneck_distance/bottleneck_basic_example.cpp @@ -0,0 +1,50 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 INRIA + * + * 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/Bottleneck.h> + +#include <iostream> +#include <vector> +#include <utility> // for pair +#include <limits> // for numeric_limits + +int main() { + std::vector< std::pair<double, double> > v1, v2; + + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits<double>::infinity()); + + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits<double>::infinity()); + + + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + + std::cout << "Approx bottleneck distance = " << b << std::endl; +} diff --git a/example/Bottleneck_distance/bottleneck_read_file_example.cpp b/example/Bottleneck_distance/bottleneck_read_file_example.cpp new file mode 100644 index 00000000..bde05825 --- /dev/null +++ b/example/Bottleneck_distance/bottleneck_read_file_example.cpp @@ -0,0 +1,72 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 INRIA + * + * 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/>. + */ + +#define CGAL_HAS_THREADS + +#include <gudhi/Bottleneck.h> +#include <iostream> +#include <vector> +#include <utility> // for pair +#include <fstream> +#include <sstream> +#include <string> + +std::vector< std::pair<double, double> > read_diagram_from_file(const char* filename) { + std::ifstream in; + in.open(filename); + std::vector< std::pair<double, double> > result; + if (!in.is_open()) { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } + + std::string line; + while (!in.eof()) { + getline(in, line); + if (line.length() != 0) { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back(std::make_pair(beginn, endd)); + } + } + in.close(); + return result; +} // read_diagram_from_file + +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << + "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << + " distance (set by default to zero). The program will now terminate \n"; + } + std::vector< std::pair< double, double > > diag1 = read_diagram_from_file(argv[1]); + std::vector< std::pair< double, double > > diag2 = read_diagram_from_file(argv[2]); + double tolerance = 0.; + if (argc == 4) { + tolerance = atof(argv[3]); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; +} diff --git a/example/Contraction/CMakeLists.txt b/example/Contraction/CMakeLists.txt index 4c09a0a7..51a6832d 100644 --- a/example/Contraction/CMakeLists.txt +++ b/example/Contraction/CMakeLists.txt @@ -9,7 +9,10 @@ target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBR target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) -add_test(RipsContraction.tore3D.0.2 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 0.2) +add_test(NAME Contraction_example_tore3D_0.2 COMMAND $<TARGET_FILE:RipsContraction> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "0.2") # TODO(DS) : These tests are too long under Windows -#add_test(RipsContraction.sphere.0.2 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/sphere3D_2646.off 0.2) -#add_test(RipsContraction.S0310000 ${CMAKE_CURRENT_BINARY_DIR}/RipsContraction ${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off 0.3) +#add_test(NAME Contraction_example_sphere_0.2 COMMAND $<TARGET_FILE:RipsContraction> +# "${CMAKE_SOURCE_DIR}/data/points/sphere3D_2646.off" "0.2") +#add_test(NAME Contraction_example_SO3_0.3 COMMAND $<TARGET_FILE:RipsContraction> +# "${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off" "0.3") diff --git a/example/Contraction/Garland_heckbert.cpp b/example/Contraction/Garland_heckbert.cpp index cbc46e91..8b5a6a6c 100644 --- a/example/Contraction/Garland_heckbert.cpp +++ b/example/Contraction/Garland_heckbert.cpp @@ -35,11 +35,6 @@ #include "Garland_heckbert/Error_quadric.h" -using namespace std; -using namespace Gudhi; -using namespace skeleton_blocker; -using namespace contraction; - struct Geometry_trait { typedef Point_d Point; }; @@ -47,7 +42,8 @@ struct Geometry_trait { /** * The vertex stored in the complex contains a quadric. */ -struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits<Geometry_trait> { +struct Garland_heckbert_traits + : public Gudhi::skeleton_blocker::Skeleton_blocker_simple_geometric_traits<Geometry_trait> { public: struct Garland_heckbert_vertex : public Simple_geometric_vertex { Error_quadric<Geometry_trait::Point> quadric; @@ -55,9 +51,9 @@ struct Garland_heckbert_traits : public Skeleton_blocker_simple_geometric_traits typedef Garland_heckbert_vertex Graph_vertex; }; -typedef Skeleton_blocker_geometric_complex< Garland_heckbert_traits > Complex; -typedef Edge_profile<Complex> EdgeProfile; -typedef Skeleton_blocker_contractor<Complex> Complex_contractor; +using Complex = Gudhi::skeleton_blocker::Skeleton_blocker_geometric_complex< Garland_heckbert_traits >; +using EdgeProfile = Gudhi::contraction::Edge_profile<Complex>; +using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor<Complex>; /** * How the new vertex is placed after an edge collapse : here it is placed at @@ -69,7 +65,7 @@ class GH_placement : public Gudhi::contraction::Placement_policy<EdgeProfile> { public: typedef Gudhi::contraction::Placement_policy<EdgeProfile>::Placement_type Placement_type; - GH_placement(Complex& complex) : complex_(complex) { } + GH_placement(Complex& complex) : complex_(complex) { (void)complex_; } Placement_type operator()(const EdgeProfile& profile) const override { auto sum_quad(profile.v0().quadric); @@ -93,7 +89,7 @@ class GH_cost : public Gudhi::contraction::Cost_policy<EdgeProfile> { public: typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; - GH_cost(Complex& complex) : complex_(complex) { } + GH_cost(Complex& complex) : complex_(complex) { (void)complex_; } Cost_type operator()(EdgeProfile const& profile, boost::optional<Point> const& new_point) const override { Cost_type res; @@ -115,7 +111,7 @@ class GH_visitor : public Gudhi::contraction::Contraction_visitor<EdgeProfile> { Complex& complex_; public: - GH_visitor(Complex& complex) : complex_(complex) { } + GH_visitor(Complex& complex) : complex_(complex) { (void)complex_; } // Compute quadrics for every vertex v // The quadric of v consists in the sum of quadric @@ -154,7 +150,7 @@ int main(int argc, char *argv[]) { typedef Complex::Vertex_handle Vertex_handle; // load the points - Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex); + Gudhi::skeleton_blocker::Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex); if (!off_reader.is_valid()) { std::cerr << "Unable to read file:" << argv[1] << std::endl; return EXIT_FAILURE; @@ -175,7 +171,7 @@ int main(int argc, char *argv[]) { Complex_contractor contractor(complex, new GH_cost(complex), new GH_placement(complex), - contraction::make_link_valid_contraction<EdgeProfile>(), + Gudhi::contraction::make_link_valid_contraction<EdgeProfile>(), new GH_visitor(complex)); std::cout << "Contract " << num_contractions << " edges" << std::endl; @@ -187,7 +183,7 @@ int main(int argc, char *argv[]) { complex.num_triangles() << " triangles." << std::endl; // write simplified complex - Skeleton_blocker_off_writer<Complex> off_writer(argv[2], complex); + Gudhi::skeleton_blocker::Skeleton_blocker_off_writer<Complex> off_writer(argv[2], complex); return EXIT_SUCCESS; } diff --git a/example/Contraction/Garland_heckbert/Error_quadric.h b/example/Contraction/Garland_heckbert/Error_quadric.h index 076f1be0..e7dafaa0 100644 --- a/example/Contraction/Garland_heckbert/Error_quadric.h +++ b/example/Contraction/Garland_heckbert/Error_quadric.h @@ -21,8 +21,8 @@ * */ -#ifndef ERROR_QUADRIC_H_ -#define ERROR_QUADRIC_H_ +#ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_ +#define GARLAND_HECKBERT_ERROR_QUADRIC_H_ #include <boost/optional/optional.hpp> @@ -179,4 +179,4 @@ template <typename Point> class Error_quadric { } }; -#endif // ERROR_QUADRIC_H_ +#endif // GARLAND_HECKBERT_ERROR_QUADRIC_H_ diff --git a/example/Contraction/Rips_contraction.cpp b/example/Contraction/Rips_contraction.cpp index 978dd1cb..8289b1d3 100644 --- a/example/Contraction/Rips_contraction.cpp +++ b/example/Contraction/Rips_contraction.cpp @@ -27,20 +27,15 @@ #include <boost/timer/timer.hpp> #include <iostream> -using namespace std; -using namespace Gudhi; -using namespace skeleton_blocker; -using namespace contraction; - struct Geometry_trait { typedef Point_d Point; }; -typedef Geometry_trait::Point Point; -typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits; -typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex; -typedef Edge_profile<Complex> Profile; -typedef Skeleton_blocker_contractor<Complex> Complex_contractor; +using Complex_geometric_traits = Gudhi::skeleton_blocker::Skeleton_blocker_simple_geometric_traits<Geometry_trait>; +using Complex = Gudhi::skeleton_blocker::Skeleton_blocker_geometric_complex< Complex_geometric_traits >; +using Profile = Gudhi::contraction::Edge_profile<Complex>; +using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor<Complex>; + template<typename ComplexType> void build_rips(ComplexType& complex, double offset) { @@ -63,7 +58,7 @@ int main(int argc, char *argv[]) { Complex complex; // load only the points - Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex, true); + Gudhi::skeleton_blocker::Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex, true); if (!off_reader.is_valid()) { std::cerr << "Unable to read file:" << argv[1] << std::endl; return EXIT_FAILURE; @@ -80,10 +75,10 @@ int main(int argc, char *argv[]) { complex.num_edges() << " edges" << std::endl; Complex_contractor contractor(complex, - new Edge_length_cost<Profile>, - contraction::make_first_vertex_placement<Profile>(), - contraction::make_link_valid_contraction<Profile>(), - contraction::make_remove_popable_blockers_visitor<Profile>()); + new Gudhi::contraction::Edge_length_cost<Profile>, + Gudhi::contraction::make_first_vertex_placement<Profile>(), + Gudhi::contraction::make_link_valid_contraction<Profile>(), + Gudhi::contraction::make_remove_popable_blockers_visitor<Profile>()); contractor.contract_edges(); std::cout << "Counting final number of simplices \n"; diff --git a/example/Persistent_cohomology/CMakeLists.txt b/example/Persistent_cohomology/CMakeLists.txt index d97d1b63..3c45e79b 100644 --- a/example/Persistent_cohomology/CMakeLists.txt +++ b/example/Persistent_cohomology/CMakeLists.txt @@ -1,19 +1,21 @@ cmake_minimum_required(VERSION 2.6) project(Persistent_cohomology_examples) -# problem with Visual Studio link on Boost program_options -add_definitions( -DBOOST_ALL_NO_LIB ) -add_definitions( -DBOOST_ALL_DYN_LINK ) - add_executable(plain_homology plain_homology.cpp) target_link_libraries(plain_homology ${Boost_SYSTEM_LIBRARY}) add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp) target_link_libraries(persistence_from_simple_simplex_tree ${Boost_SYSTEM_LIBRARY}) +add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp) +target_link_libraries(rips_distance_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) +add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp) +target_link_libraries(rips_persistence_step_by_step ${Boost_SYSTEM_LIBRARY} ${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_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) @@ -23,62 +25,82 @@ target_link_libraries(persistence_from_file ${Boost_SYSTEM_LIBRARY} ${Boost_PROG if (TBB_FOUND) target_link_libraries(plain_homology ${TBB_LIBRARIES}) target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence ${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(plain_homology ${CMAKE_CURRENT_BINARY_DIR}/plain_homology) -add_test(persistence_from_simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_simple_simplex_tree 1 0) -add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.2 -d 3 -p 3 -m 100) -add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_via_boundary_matrix ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.txt -r 0.3 -d 3 -p 3 -m 100) -add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 2 -m 0) -add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100) +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_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence> + "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") +add_test(NAME Persistent_cohomology_example_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +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") if(GMP_FOUND) if(GMPXX_FOUND) add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) - target_link_libraries(rips_multifield_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) - add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) - target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) + target_link_libraries(rips_multifield_persistence + ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if (TBB_FOUND) target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) - target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) - - add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.2 -d 3 -p 2 -q 71 -m 100) + 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") endif(GMPXX_FOUND) endif(GMP_FOUND) if(CGAL_FOUND) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) target_link_libraries(alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) + target_link_libraries(exact_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) - add_test(alpha_complex_3d_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off 2 0.45) - - - if (NOT CGAL_VERSION VERSION_LESS 4.7.0) - if (EIGEN3_FOUND) - add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) - target_link_libraries(alpha_complex_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) - - add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp) - target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - - add_executable(custom_persistence_sort custom_persistence_sort.cpp) - target_link_libraries(custom_persistence_sort ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - - if (TBB_FOUND) - target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) - target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) - target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES}) - endif(TBB_FOUND) - add_test(alpha_complex_persistence_2_0_45 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -m 0.45 -p 2) - add_test(periodic_alpha_complex_3d_persistence_2_0 ${CMAKE_CURRENT_BINARY_DIR}/periodic_alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off ${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt 2 0) - add_test(custom_persistence_sort ${CMAKE_CURRENT_BINARY_DIR}/custom_persistence_sort) - endif(EIGEN3_FOUND) - endif (NOT CGAL_VERSION VERSION_LESS 4.7.0) + add_test(NAME Persistent_cohomology_example_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + add_test(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $<TARGET_FILE:exact_alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") + + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) + target_link_libraries(alpha_complex_persistence + ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp) + target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + + add_executable(custom_persistence_sort custom_persistence_sort.cpp) + target_link_libraries(custom_persistence_sort ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistent_cohomology_example_alpha_complex COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") + add_test(NAME Persistent_cohomology_example_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:periodic_alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "2" "0") + add_test(NAME Persistent_cohomology_example_custom_persistence_sort COMMAND $<TARGET_FILE:custom_persistence_sort>) + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) endif(CGAL_FOUND) diff --git a/example/Persistent_cohomology/README b/example/Persistent_cohomology/README index 7803e5ab..2ac79398 100644 --- a/example/Persistent_cohomology/README +++ b/example/Persistent_cohomology/README @@ -10,13 +10,13 @@ Example of use of RIPS: 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 +./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2 output: -210 0 0 inf -210 1 0.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +2 0 0 inf +2 1 0.0983494 inf +2 1 0.104347 inf +2 2 0.138335 inf Every line is of this format: p1*...*pr dim b d @@ -29,31 +29,45 @@ where with Z/3Z coefficients: -./rips_persistence ../../data/points/Kl.txt -r 0.25 -d 3 -p 3 -m 100 +./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3 output: -3 0 0 inf -3 1 0.0702103 inf +3 0 0 inf +3 1 0.0983494 inf +3 1 0.104347 inf +3 2 0.138335 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 +./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.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +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 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 + ./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.0702103 inf -2 1 0.0702103 inf -2 2 0.159992 inf +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 ALPHA: diff --git a/example/Persistent_cohomology/alpha_complex_3d_helper.h b/example/Persistent_cohomology/alpha_complex_3d_helper.h new file mode 100644 index 00000000..7865e4ec --- /dev/null +++ b/example/Persistent_cohomology/alpha_complex_3d_helper.h @@ -0,0 +1,76 @@ +/* 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/>. + */ + +#ifndef ALPHA_COMPLEX_3D_HELPER_H_ +#define ALPHA_COMPLEX_3D_HELPER_H_ + +template<class Vertex_list, class Cell_handle> +Vertex_list from_cell(const Cell_handle& ch) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { +#ifdef DEBUG_TRACES + std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(ch->vertex(i)); + } + return the_list; +} + +template<class Vertex_list, class Facet> +Vertex_list from_facet(const Facet& fct) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { + if (fct.second != i) { +#ifdef DEBUG_TRACES + std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(fct.first->vertex(i)); + } + } + return the_list; +} + +template<class Vertex_list, class Edge_3> +Vertex_list from_edge(const Edge_3& edg) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { + if ((edg.second == i) || (edg.third == i)) { +#ifdef DEBUG_TRACES + std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(edg.first->vertex(i)); + } + } + return the_list; +} + +template<class Vertex_list, class Vertex_handle> +Vertex_list from_vertex(const Vertex_handle& vh) { + Vertex_list the_list; +#ifdef DEBUG_TRACES + std::cout << "from vertex=" << vh->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(vh); + return the_list; +} + +#endif // ALPHA_COMPLEX_3D_HELPER_H_ diff --git a/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp index 48fbb91a..fd227b82 100644 --- a/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp +++ b/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Saclay (France) + * Copyright (C) 2014 INRIA * * 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 @@ -20,10 +20,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#include <boost/variant.hpp> + #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> #include <gudhi/Points_3D_off_io.h> -#include <boost/variant.hpp> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Delaunay_triangulation_3.h> @@ -39,84 +40,42 @@ #include <list> #include <vector> +#include "alpha_complex_3d_helper.h" + // Alpha_shape_3 templates type definitions -typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef CGAL::Alpha_shape_vertex_base_3<Kernel> Vb; -typedef CGAL::Alpha_shape_cell_base_3<Kernel> Fb; -typedef CGAL::Triangulation_data_structure_3<Vb, Fb> Tds; -typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> Triangulation_3; -typedef CGAL::Alpha_shape_3<Triangulation_3> Alpha_shape_3; +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>; +using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; // From file type definition -typedef Kernel::Point_3 Point_3; +using Point_3 = Kernel::Point_3; // filtration with alpha values needed type definition -typedef Alpha_shape_3::FT Alpha_value_type; -typedef CGAL::Object Object; -typedef CGAL::Dispatch_output_iterator< -CGAL::cpp11::tuple<Object, Alpha_value_type>, -CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > > Dispatch; -typedef Alpha_shape_3::Cell_handle Cell_handle; -typedef Alpha_shape_3::Facet Facet; -typedef Alpha_shape_3::Edge Edge_3; -typedef std::list<Alpha_shape_3::Vertex_handle> Vertex_list; +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; // gudhi type definition -typedef Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence> ST; -typedef ST::Vertex_handle Simplex_tree_vertex; -typedef std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex > Alpha_shape_simplex_tree_map; -typedef std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex> Alpha_shape_simplex_tree_pair; -typedef std::vector< Simplex_tree_vertex > Simplex_tree_vector_vertex; -typedef Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp > PCOH; - -Vertex_list from(const Cell_handle& ch) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { -#ifdef DEBUG_TRACES - std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(ch->vertex(i)); - } - return the_list; -} - -Vertex_list from(const Facet& fct) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if (fct.second != i) { -#ifdef DEBUG_TRACES - std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(fct.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Edge_3& edg) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if ((edg.second == i) || (edg.third == i)) { -#ifdef DEBUG_TRACES - std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(edg.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) { - Vertex_list the_list; -#ifdef DEBUG_TRACES - std::cout << "from vertex=" << vh->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(vh); - return the_list; -} - -void usage(char * const progName) { +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(const std::string& progName) { std::cerr << "Usage: " << progName << " path_to_file_graph coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); @@ -132,7 +91,7 @@ int main(int argc, char * const argv[]) { int coeff_field_characteristic = atoi(argv[2]); Filtration_value min_persistence = 0.0; - int returnedScanValue = sscanf(argv[3], "%lf", &min_persistence); + int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { std::cerr << "Error: " << argv[3] << " is not correct\n"; usage(argv[0]); @@ -184,30 +143,29 @@ int main(int argc, char * const argv[]) { for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { - vertex_list = from(*cell); + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { - vertex_list = from(*facet); + vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { - vertex_list = from(*edge); + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Alpha_shape_3::Vertex_handle * vertex = - CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; - vertex_list = from(*vertex); + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex Simplex_tree_vector_vertex the_simplex_tree; diff --git a/example/Persistent_cohomology/alpha_complex_persistence.cpp b/example/Persistent_cohomology/alpha_complex_persistence.cpp index cb181936..9e84e91f 100644 --- a/example/Persistent_cohomology/alpha_complex_persistence.cpp +++ b/example/Persistent_cohomology/alpha_complex_persistence.cpp @@ -4,11 +4,16 @@ #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 <string> #include <limits> // for numeric_limits +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; + void program_options(int argc, char * argv[] , std::string & off_file_points , std::string & output_file_diag @@ -30,35 +35,38 @@ int main(int argc, char **argv) { // Init of an alpha complex from an OFF file // ---------------------------------------------------------------------------- using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; - Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points, alpha_square_max_value); - - // ---------------------------------------------------------------------------- - // Display information about the alpha complex - // ---------------------------------------------------------------------------- - std::cout << "Alpha complex is of dimension " << alpha_complex_from_file.dimension() << - " - " << alpha_complex_from_file.num_simplices() << " simplices - " << - alpha_complex_from_file.num_vertices() << " vertices." << std::endl; - - // Sort the simplices in the order of the filtration - alpha_complex_from_file.initialize_filtration(); - - std::cout << "Simplex_tree dim: " << alpha_complex_from_file.dimension() << std::endl; - // Compute the persistence diagram of the complex - Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::alpha_complex::Alpha_complex<Kernel>, - Gudhi::persistent_cohomology::Field_Zp > pcoh(alpha_complex_from_file); - // initializes the coefficient field for homology - pcoh.init_coefficients(coeff_field_characteristic); - - pcoh.compute_persistent_cohomology(min_persistence); - - // Output the diagram in filediag - if (output_file_diag.empty()) { - pcoh.output_diagram(); - } else { - std::cout << "Result in file: " << output_file_diag << std::endl; - std::ofstream out(output_file_diag); - pcoh.output_diagram(out); - out.close(); + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points); + + Simplex_tree simplex; + if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + // ---------------------------------------------------------------------------- + // 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; + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, + Gudhi::persistent_cohomology::Field_Zp > pcoh(simplex); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (output_file_diag.empty()) { + pcoh.output_diagram(); + } else { + std::cout << "Result in file: " << output_file_diag << std::endl; + std::ofstream out(output_file_diag); + pcoh.output_diagram(out); + out.close(); + } } return 0; diff --git a/example/Persistent_cohomology/custom_persistence_sort.cpp b/example/Persistent_cohomology/custom_persistence_sort.cpp index 9af38611..64f2a4dc 100644 --- a/example/Persistent_cohomology/custom_persistence_sort.cpp +++ b/example/Persistent_cohomology/custom_persistence_sort.cpp @@ -27,6 +27,8 @@ #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> @@ -38,6 +40,9 @@ 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 @@ -60,7 +65,7 @@ std::vector<Point> random_points() { * Compare two intervals by dimension, then by length. */ struct cmp_intervals_by_dim_then_length { - explicit cmp_intervals_by_dim_then_length(Alpha_complex * sc) + explicit cmp_intervals_by_dim_then_length(Simplex_tree * sc) : sc_(sc) { } template<typename Persistent_interval> @@ -71,46 +76,62 @@ struct cmp_intervals_by_dim_then_length { else return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); } - Alpha_complex* sc_; + 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, 0.6); - - using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Alpha_complex, - Gudhi::persistent_cohomology::Field_Zp >; - Persistent_cohomology pcoh(alpha_complex_from_points); - - // initializes the coefficient field for homology - Z/3Z - pcoh.init_coefficients(3); - pcoh.compute_persistent_cohomology(0.2); - - // Custom sort and output persistence - cmp_intervals_by_dim_then_length cmp(&alpha_complex_from_points); - auto persistent_pairs = pcoh.get_persistent_pairs(); - std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); - for (auto pair : persistent_pairs) { - std::cout << alpha_complex_from_points.dimension(get<0>(pair)) << " " - << alpha_complex_from_points.filtration(get<0>(pair)) << " " - << alpha_complex_from_points.filtration(get<1>(pair)) << std::endl; + 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; } - - // Persistent Betti numbers - std::cout << "The persistent Betti numbers in interval [0.40, 0.41] are : "; - for (int dim = 0; dim < alpha_complex_from_points.dimension(); dim++) - std::cout << "b" << dim << " = " << pcoh.persistent_betti_number(dim, 0.40, 0.41) << " ; "; - std::cout << std::endl; - - // Betti numbers - std::vector<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/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..8a335075 --- /dev/null +++ b/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp @@ -0,0 +1,245 @@ +/* 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 <boost/variant.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Delaunay_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/iterator.h> + +#include <fstream> +#include <cmath> +#include <string> +#include <tuple> +#include <map> +#include <utility> +#include <list> +#include <vector> + +#include "alpha_complex_3d_helper.h" + +// Alpha_shape_3 templates type definitions +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Exact_tag = CGAL::Tag_true; +using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, CGAL::Default, Exact_tag>; +using Fb = CGAL::Alpha_shape_cell_base_3<Kernel, CGAL::Default, Exact_tag>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3, Exact_tag>; + +// From file type definition +using Point_3 = Kernel::Point_3; + +// filtration with alpha values needed type definition +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Vertex_handle>; + +// gudhi type definition +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << + " path_to_file_graph 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 != 4) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[2]); + + Filtration_value min_persistence = 0.0; + int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); + if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { + std::cerr << "Error: " << argv[3] << " is not correct\n"; + usage(argv[0]); + } + + // Read points from file + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_3> lp = off_reader.get_point_cloud(); + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. + Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL); +#ifdef DEBUG_TRACES + std::cout << "Alpha shape computed in GENERAL mode" << std::endl; +#endif // DEBUG_TRACES + + // filtration with alpha values from alpha shape + std::vector<Object> the_objects; + std::vector<Alpha_value_type> the_alpha_values; + + Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; + + // Loop on objects vector + Vertex_list vertex_list; + ST simplex_tree; + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin(); + int dim_max = 0; + Filtration_value filtration_max = 0.0; + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + vertex_list = from_facet<Vertex_list, Facet>(*facet); + count_facets++; + if (dim_max < 2) { + // Facet is of dim 2 + dim_max = 2; + } + } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); + count_edges++; + if (dim_max < 1) { + // Edge_3 is of dim 1 + dim_max = 1; + } + } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex_tree; + for (auto the_alpha_shape_vertex : vertex_list) { + Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Simplex_tree_vertex vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex)); + } else { + // alpha shape found + Simplex_tree_vertex vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + } + } + // Construction of the simplex_tree + // you can also use the_alpha_value_iterator->exact() + Filtration_value filtr = /*std::sqrt*/CGAL::to_double(the_alpha_value_iterator->exact()); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + if (filtr > filtration_max) { + filtration_max = filtr; + } + simplex_tree.insert_simplex(the_simplex_tree, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + simplex_tree.set_filtration(filtration_max); + simplex_tree.set_dimension(dim_max); + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; + std::cout << " filtration = " << simplex_tree.filtration() << std::endl << std::endl; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + // Compute the persistence diagram of the complex + PCOH pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + pcoh.output_diagram(); + + return 0; +} diff --git a/example/Persistent_cohomology/performance_rips_persistence.cpp b/example/Persistent_cohomology/performance_rips_persistence.cpp deleted file mode 100644 index b4d282ac..00000000 --- a/example/Persistent_cohomology/performance_rips_persistence.cpp +++ /dev/null @@ -1,214 +0,0 @@ -/* 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> -#include <string> -#include <vector> - -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -/* Compute the persistent homology of the complex cpx with coefficients in Z/pZ. */ -template< typename FilteredComplex> -void timing_persistence(FilteredComplex & cpx - , int p); - -/* Compute multi-field persistent homology of the complex cpx with coefficients in - * Z/rZ for all prime number r in [p;q].*/ -template< typename FilteredComplex> -void timing_persistence(FilteredComplex & cpx - , int p - , int q); - -/* Timings for the computation of persistent homology with different - * representations of a Rips complex and different coefficient fields. The - * Rips complex is built on a set of 10000 points sampling a Klein bottle embedded - * in dimension 5. - * We represent complexes with a simplex tree and - * with a Hasse diagram. The Hasse diagram represents explicitly all - * codimension 1 incidence relations in the complex, and hence leads to - * a faster computation of persistence because boundaries are precomputed. - * Hovewer, the simplex tree may be constructed directly from a point cloud and - * is more compact. - * We compute persistent homology with coefficient fields Z/2Z and Z/1223Z. - * We present also timings for the computation of multi-field persistent - * homology in all fields Z/rZ for r prime between 2 and 1223. - */ -int main(int argc, char * argv[]) { - std::chrono::time_point<std::chrono::system_clock> start, end; - int elapsed_sec; - { - - std::string filepoints = "../../../data/points/Kl.txt"; - Filtration_value threshold = 0.27; - 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(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; - - // Construct the Rips complex in a Simplex Tree - Simplex_tree<Simplex_tree_options_fast_persistence> st; - start = std::chrono::system_clock::now(); - - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Compute Rips complex in " << elapsed_sec << " ms.\n"; - std::cout << " - dimension = " << st.dimension() << std::endl; - std::cout << " - number of simplices = " << st.num_simplices() << std::endl; - - // Sort the simplices in the order of the filtration - start = std::chrono::system_clock::now(); - st.initialize_filtration(); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Order the simplices of the filtration in " << elapsed_sec << " ms.\n"; - - // Copy the keys inside the simplices - start = std::chrono::system_clock::now(); - { - int count = 0; - for (auto sh : st.filtration_simplex_range()) - st.assign_key(sh, count++); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Copied the keys inside the simplices in " << elapsed_sec << " ms.\n"; - - // Convert the simplex tree into a hasse diagram - start = std::chrono::system_clock::now(); - Hasse_complex<> hcpx(st); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Convert the simplex tree into a Hasse diagram in " << elapsed_sec << " ms.\n"; - - - std::cout << "Timings when using a simplex tree: \n"; - timing_persistence(st, p); - timing_persistence(st, q); - timing_persistence(st, p, q); - - std::cout << "Timings when using a Hasse complex: \n"; - timing_persistence(hcpx, p); - timing_persistence(hcpx, q); - timing_persistence(hcpx, p, q); - - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << "Running the complex destructors in " << elapsed_sec << " ms.\n"; - return 0; -} - -template< typename FilteredComplex> -void -timing_persistence(FilteredComplex & cpx - , int p) { - std::chrono::time_point<std::chrono::system_clock> start, end; - int elapsed_sec; - { - start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Field_Zp > pcoh(cpx); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; - // initializes the coefficient field for homology - start = std::chrono::system_clock::now(); - pcoh.init_coefficients(p); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; - - start = std::chrono::system_clock::now(); - - pcoh.compute_persistent_cohomology(INFINITY); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Compute persistent homology in Z/" << p << "Z in " << elapsed_sec << " ms.\n"; - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; -} - -template< typename FilteredComplex> -void -timing_persistence(FilteredComplex & cpx - , int p - , int q) { - std::chrono::time_point<std::chrono::system_clock> start, end; - int elapsed_sec; - { - start = std::chrono::system_clock::now(); - Persistent_cohomology< FilteredComplex, Multi_field > pcoh(cpx); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Initialize pcoh in " << elapsed_sec << " ms.\n"; - // initializes the coefficient field for homology - start = std::chrono::system_clock::now(); - pcoh.init_coefficients(p, q); - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Initialize the coefficient field in " << elapsed_sec << " ms.\n"; - // compute persistent homology, disgarding persistent features of life shorter than min_persistence - - start = std::chrono::system_clock::now(); - - pcoh.compute_persistent_cohomology(INFINITY); - - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Compute multi-field persistent homology in all coefficient fields Z/pZ " - << "with p in [" << p << ";" << q << "] in " << elapsed_sec << " ms.\n"; - start = std::chrono::system_clock::now(); - } - end = std::chrono::system_clock::now(); - elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); - std::cout << " Run the persistence destructors in " << elapsed_sec << " ms.\n"; -} diff --git a/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp index a199fea1..8928cfc2 100644 --- a/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp +++ b/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Saclay (France) + * Copyright (C) 2014 INRIA * * 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 @@ -20,10 +20,11 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ +#include <boost/variant.hpp> + #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> #include <gudhi/Points_3D_off_io.h> -#include <boost/variant.hpp> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h> @@ -39,6 +40,9 @@ #include <utility> #include <list> #include <vector> +#include <cstdlib> + +#include "alpha_complex_3d_helper.h" // Traits using K = CGAL::Exact_predicates_inexact_constructions_kernel; @@ -66,10 +70,12 @@ using Dispatch = CGAL::Dispatch_output_iterator< using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; // gudhi type definition using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; @@ -77,52 +83,6 @@ using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; -Vertex_list from(const Cell_handle& ch) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { -#ifdef DEBUG_TRACES - std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(ch->vertex(i)); - } - return the_list; -} - -Vertex_list from(const Facet& fct) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if (fct.second != i) { -#ifdef DEBUG_TRACES - std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(fct.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Edge_3& edg) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if ((edg.second == i) || (edg.third == i)) { -#ifdef DEBUG_TRACES - std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(edg.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) { - Vertex_list the_list; -#ifdef DEBUG_TRACES - std::cout << "from vertex=" << vh->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(vh); - return the_list; -} - void usage(char * const progName) { std::cerr << "Usage: " << progName << " path_to_file_graph path_to_iso_cuboid_3_file coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; @@ -136,19 +96,8 @@ int main(int argc, char * const argv[]) { usage(argv[0]); } - int coeff_field_characteristic = 0; - int returnedScanValue = sscanf(argv[3], "%d", &coeff_field_characteristic); - if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) { - std::cerr << "Error: " << argv[3] << " is not correct\n"; - usage(argv[0]); - } - - Filtration_value min_persistence = 0.0; - returnedScanValue = sscanf(argv[4], "%lf", &min_persistence); - if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { - std::cerr << "Error: " << argv[4] << " is not correct\n"; - usage(argv[0]); - } + int coeff_field_characteristic = atoi(argv[3]); + Filtration_value min_persistence = strtof(argv[4], nullptr); // Read points from file std::string offInputFile(argv[1]); @@ -212,21 +161,21 @@ int main(int argc, char * const argv[]) { for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { - vertex_list = from(*cell); + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { - vertex_list = from(*facet); + vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { - vertex_list = from(*edge); + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 @@ -235,7 +184,7 @@ int main(int argc, char * const argv[]) { } else if (const Alpha_shape_3::Vertex_handle * vertex = CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { count_vertices++; - vertex_list = from(*vertex); + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex Simplex_tree_vector_vertex the_simplex_tree; diff --git a/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp b/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp index ba772f04..7ca9410a 100644 --- a/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp +++ b/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * 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 @@ -29,13 +29,12 @@ #include <utility> #include <vector> -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; -typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; -typedef Simplex_tree<> typeST; +// 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"; @@ -66,7 +65,7 @@ int main(int argc, char * const argv[]) { // TEST OF INSERTION std::cout << "********************************************************************" << std::endl; std::cout << "TEST OF INSERTION" << std::endl; - typeST st; + Simplex_tree st; // ++ FIRST std::cout << " - INSERT (0,1,2)" << std::endl; @@ -166,7 +165,7 @@ int main(int argc, char * const argv[]) { std::cout << "**************************************************************" << std::endl; // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(st); + Persistent_cohomology pcoh(st); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); diff --git a/example/Persistent_cohomology/plain_homology.cpp b/example/Persistent_cohomology/plain_homology.cpp index ae82c817..50f692f2 100644 --- a/example/Persistent_cohomology/plain_homology.cpp +++ b/example/Persistent_cohomology/plain_homology.cpp @@ -27,13 +27,11 @@ #include <vector> #include <cstdint> // for std::uint8_t -using namespace Gudhi; - /* 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 : Simplex_tree_options_full_featured { +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 @@ -43,7 +41,10 @@ struct MyOptions : Simplex_tree_options_full_featured { // Maximum number of simplices to compute persistence is 2^8 - 1 = 255. One is reserved for null_key typedef std::uint8_t Simplex_key; }; -typedef Simplex_tree<MyOptions> ST; + +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; @@ -70,7 +71,7 @@ int main() { st.initialize_filtration(); // Class for homology computation - persistent_cohomology::Persistent_cohomology<ST, persistent_cohomology::Field_Zp> pcoh(st); + Persistent_cohomology pcoh(st); // Initialize the coefficient field Z/2Z for homology pcoh.init_coefficients(2); diff --git a/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp b/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp new file mode 100644 index 00000000..8517e7f6 --- /dev/null +++ b/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp @@ -0,0 +1,144 @@ +/* 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): Pawel Dlotko, Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * 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/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/reader_utils.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Distance_matrix = std::vector<std::vector<Filtration_value>>; + +void program_options(int argc, char * argv[] + , std::string & csv_matrix_file + , std::string & filediag + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); + + Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file); + Rips_complex rips_complex_from_file(distances, threshold); + + // 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(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 & csv_matrix_file + , 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>(&csv_matrix_file), + "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'."); + + 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 distance matrix.\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/example/Persistent_cohomology/rips_multifield_persistence.cpp b/example/Persistent_cohomology/rips_multifield_persistence.cpp index c5cd775d..dae36ed2 100644 --- a/example/Persistent_cohomology/rips_multifield_persistence.cpp +++ b/example/Persistent_cohomology/rips_multifield_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * 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 @@ -20,26 +20,29 @@ * 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/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> -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// 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 & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -48,7 +51,7 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; std::string filediag; Filtration_value threshold; int dim_max; @@ -56,33 +59,26 @@ int main(int argc, char * argv[]) { 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); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance<Point_t>); + 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 - typedef Simplex_tree<Simplex_tree_options_fast_persistence> ST; - ST 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); + 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 - st.initialize_filtration(); + simplex_tree.initialize_filtration(); // Compute the persistence diagram of the complex - Persistent_cohomology<ST, Multi_field > pcoh(st); + Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(min_p, max_p); - // compute persistent homology, disgarding persistent features of life shorter than min_persistence + pcoh.compute_persistent_cohomology(min_persistence); // Output the diagram in filediag @@ -98,7 +94,7 @@ int main(int argc, char * argv[]) { } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -108,8 +104,8 @@ void program_options(int argc, char * argv[] 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"); + ("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() diff --git a/example/Persistent_cohomology/rips_persistence.cpp b/example/Persistent_cohomology/rips_persistence.cpp index cab49395..d504798b 100644 --- a/example/Persistent_cohomology/rips_persistence.cpp +++ b/example/Persistent_cohomology/rips_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * 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 @@ -20,11 +20,11 @@ * 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/Rips_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> @@ -32,14 +32,17 @@ #include <vector> #include <limits> // infinity -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using 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 & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -47,40 +50,30 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; 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); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance<Point_t>); + 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 - typedef Simplex_tree<Simplex_tree_options_fast_persistence> ST; - ST 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); + Simplex_tree simplex_tree; - std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; - std::cout << " and has dimension " << st.dimension() << " \n"; + 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 - st.initialize_filtration(); + simplex_tree.initialize_filtration(); // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology<ST, Field_Zp > pcoh(st); + Persistent_cohomology pcoh(simplex_tree); // initializes the coefficient field for homology pcoh.init_coefficients(p); @@ -99,7 +92,7 @@ int main(int argc, char * argv[]) { } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -108,15 +101,16 @@ void program_options(int argc, char * argv[] 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 "); + ("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()), + ("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.") diff --git a/example/Persistent_cohomology/rips_persistence_step_by_step.cpp b/example/Persistent_cohomology/rips_persistence_step_by_step.cpp new file mode 100644 index 00000000..554eeba6 --- /dev/null +++ b/example/Persistent_cohomology/rips_persistence_step_by_step.cpp @@ -0,0 +1,217 @@ +/* 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/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 Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS +, boost::property < vertex_filtration_t, Filtration_value > +, boost::property < edge_filtration_t, Filtration_value > +>; +using Edge_t = std::pair< Vertex_handle, Vertex_handle >; + +template< typename InputPointRange, typename Distance > +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance); + +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 + Graph_t prox_graph = compute_proximity_graph(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; + std::abort(); + } +} + +/** Output the proximity graph of the points. + * + * If points contains n elements, the proximity graph is the graph + * with n vertices, and an edge [u,v] iff the distance function between + * points u and v is smaller than threshold. + * + * The type PointCloud furnishes .begin() and .end() methods, that return + * iterators with value_type Point. + */ +template< typename InputPointRange, typename Distance > +Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance) { + std::vector< Edge_t > edges; + std::vector< Filtration_value > edges_fil; + + Vertex_handle idx_u, idx_v; + Filtration_value fil; + idx_u = 0; + for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { + idx_v = idx_u + 1; + for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { + fil = distance(*it_u, *it_v); + if (fil <= threshold) { + edges.emplace_back(idx_u, idx_v); + edges_fil.push_back(fil); + } + } + ++idx_u; + } + + Graph_t skel_graph(edges.begin() + , edges.end() + , edges_fil.begin() + , idx_u); // number of points labeled from 0 to idx_u-1 + + auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); + + boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; + for (std::tie(vi, vi_end) = boost::vertices(skel_graph); + vi != vi_end; ++vi) { + boost::put(vertex_prop, *vi, 0.); + } + + return skel_graph; +} diff --git a/example/Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp b/example/Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp index 4c6656f5..9618f278 100644 --- a/example/Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp +++ b/example/Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp @@ -4,8 +4,7 @@ * * Author(s): Clément Maria, Marc Glisse * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France), - * 2015 INRIA Saclay Île de France) + * Copyright (C) 2014 INRIA * * 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 @@ -21,12 +20,12 @@ * 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/Rips_complex.h> #include <gudhi/Hasse_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> #include <boost/program_options.hpp> @@ -44,14 +43,16 @@ // // //////////////////////////////////////////////////////////////// -using namespace Gudhi; -using namespace Gudhi::persistent_cohomology; - -typedef int Vertex_handle; -typedef double Filtration_value; +// 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 & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -59,30 +60,21 @@ void program_options(int argc, char * argv[] , Filtration_value & min_persistence); int main(int argc, char * argv[]) { - std::string filepoints; + std::string off_file_points; 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); + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); - // Compute the proximity graph of the points - Graph_t prox_graph = compute_proximity_graph(points, threshold - , euclidean_distance<Point_t>); + 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<>; - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion(dim_max); + 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"; @@ -99,7 +91,7 @@ int main(int argc, char * argv[]) { st.assign_key(sh, count++); // Convert to a more convenient representation. - Hasse_complex<> hcpx(st); + Gudhi::Hasse_complex<> hcpx(st); #ifdef GUDHI_USE_TBB ts.terminate(); @@ -109,7 +101,7 @@ int main(int argc, char * argv[]) { delete &st; // Compute the persistence diagram of the complex - persistent_cohomology::Persistent_cohomology< Hasse_complex<>, Field_Zp > pcoh(hcpx); + Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::Hasse_complex<>, Field_Zp > pcoh(hcpx); // initializes the coefficient field for homology pcoh.init_coefficients(p); @@ -126,7 +118,7 @@ int main(int argc, char * argv[]) { } void program_options(int argc, char * argv[] - , std::string & filepoints + , std::string & off_file_points , std::string & filediag , Filtration_value & threshold , int & dim_max @@ -135,7 +127,7 @@ void program_options(int argc, char * argv[] namespace po = boost::program_options; po::options_description hidden("Hidden options"); hidden.add_options() - ("input-file", po::value<std::string>(&filepoints), + ("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); diff --git a/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..34b90933 --- /dev/null +++ b/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp @@ -0,0 +1,263 @@ +/* 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 + * + * 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 <boost/variant.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Regular_triangulation_euclidean_traits_3.h> +#include <CGAL/Regular_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/iterator.h> + +#include <fstream> +#include <cmath> +#include <string> +#include <tuple> +#include <map> +#include <utility> +#include <list> +#include <vector> +#include <cstdlib> + +#include "alpha_complex_3d_helper.h" + +// Traits +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>; +using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>; +using Fb = CGAL::Alpha_shape_cell_base_3<Gt>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Regular_triangulation_3<Gt, Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; + +// From file type definition +using Point_3 = Gt::Bare_point; +using Weighted_point_3 = Gt::Weighted_point; + +// filtration with alpha values needed type definition +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; + +// gudhi type definition +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< + ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << + " path_to_file_graph path_to_weight_file 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 != 5) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[3]); + Filtration_value min_persistence = strtof(argv[4], nullptr); + + // Read points from file + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_3> lp = off_reader.get_point_cloud(); + + // Read weights information from file + std::ifstream weights_ifstr(argv[2]); + std::vector<Weighted_point_3> wp; + if (weights_ifstr.good()) { + double weight = 0.0; + std::size_t index = 0; + // Attempt read the weight in a double format, return false if it fails + while ((weights_ifstr >> weight) && (index < lp.size())) { + wp.push_back(Weighted_point_3(lp[index], weight)); + index++; + } + if (index != lp.size()) { + std::cerr << "Bad number of weights in file " << argv[2] << std::endl; + usage(argv[0]); + } + } else { + std::cerr << "Unable to read file " << argv[2] << std::endl; + usage(argv[0]); + } + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. + Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL); +#ifdef DEBUG_TRACES + std::cout << "Alpha shape computed in GENERAL mode" << std::endl; +#endif // DEBUG_TRACES + + // filtration with alpha values from alpha shape + std::vector<Object> the_objects; + std::vector<Alpha_value_type> the_alpha_values; + + Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; + + // Loop on objects vector + Vertex_list vertex_list; + ST simplex_tree; + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin(); + int dim_max = 0; + Filtration_value filtration_max = 0.0; + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + vertex_list = from_facet<Vertex_list, Facet>(*facet); + count_facets++; + if (dim_max < 2) { + // Facet is of dim 2 + dim_max = 2; + } + } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); + count_edges++; + if (dim_max < 1) { + // Edge_3 is of dim 1 + dim_max = 1; + } + } else if (const Alpha_shape_3::Vertex_handle * vertex = + CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex_tree; + for (auto the_alpha_shape_vertex : vertex_list) { + Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Simplex_tree_vertex vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex)); + } else { + // alpha shape found + Simplex_tree_vertex vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + if (filtr > filtration_max) { + filtration_max = filtr; + } + simplex_tree.insert_simplex(the_simplex_tree, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + simplex_tree.set_filtration(filtration_max); + simplex_tree.set_dimension(dim_max); + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; + std::cout << " filtration = " << simplex_tree.filtration() << std::endl << std::endl; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree, true); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + pcoh.output_diagram(); + + return 0; +} diff --git a/example/Rips_complex/CMakeLists.txt b/example/Rips_complex/CMakeLists.txt new file mode 100644 index 00000000..8aee79e2 --- /dev/null +++ b/example/Rips_complex/CMakeLists.txt @@ -0,0 +1,58 @@ +cmake_minimum_required(VERSION 2.6) +project(Rips_complex_examples) + +# Point cloud +add_executable ( Rips_complex_example_from_off example_rips_complex_from_off_file.cpp ) +target_link_libraries(Rips_complex_example_from_off ${Boost_SYSTEM_LIBRARY}) + +add_executable ( Rips_complex_example_one_skeleton_from_points example_one_skeleton_rips_from_points.cpp ) +target_link_libraries(Rips_complex_example_one_skeleton_from_points ${Boost_SYSTEM_LIBRARY}) + +# Distance matrix +add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_one_skeleton_rips_from_distance_matrix.cpp ) +target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${Boost_SYSTEM_LIBRARY}) + +add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp ) +target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${Boost_SYSTEM_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(Rips_complex_example_from_off ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES}) +endif() + +add_test(NAME Rips_complex_example_one_skeleton_from_points + COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_points>) +add_test(NAME Rips_complex_example_one_skeleton_from_distance_matrix + COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_distance_matrix>) + +add_test(NAME Rips_complex_example_from_off_doc_12_1 COMMAND $<TARGET_FILE:Rips_complex_example_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "12.0" "1" "${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt") +add_test(NAME Rips_complex_example_from_off_doc_12_3 COMMAND $<TARGET_FILE:Rips_complex_example_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "12.0" "3" "${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_3.txt") + +add_test(NAME Rips_complex_example_from_csv_distance_matrix_doc_12_1 COMMAND $<TARGET_FILE:Rips_complex_example_from_csv_distance_matrix> + "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "12.0" "1" "${CMAKE_CURRENT_BINARY_DIR}/ripscsvreader_result_12_1.txt") +add_test(NAME Rips_complex_example_from_csv_distance_matrix_doc_12_3 COMMAND $<TARGET_FILE:Rips_complex_example_from_csv_distance_matrix> + "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "12.0" "3" "${CMAKE_CURRENT_BINARY_DIR}/ripscsvreader_result_12_3.txt") + + +if (DIFF_PATH) + # Do not forget to copy test results files in current binary dir + file(COPY "one_skeleton_rips_for_doc.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY "full_skeleton_rips_for_doc.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(Rips_complex_example_from_off_doc_12_1_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt + ${CMAKE_CURRENT_BINARY_DIR}/one_skeleton_rips_for_doc.txt) + add_test(Rips_complex_example_from_off_doc_12_3_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_3.txt + ${CMAKE_CURRENT_BINARY_DIR}/full_skeleton_rips_for_doc.txt) + add_test(Rips_complex_example_from_csv_distance_matrix_doc_12_1_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/ripscsvreader_result_12_1.txt + ${CMAKE_CURRENT_BINARY_DIR}/one_skeleton_rips_for_doc.txt) + add_test(Rips_complex_example_from_csv_distance_matrix_doc_12_3_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/ripscsvreader_result_12_3.txt + ${CMAKE_CURRENT_BINARY_DIR}/full_skeleton_rips_for_doc.txt) +endif() diff --git a/example/Rips_complex/example_one_skeleton_rips_from_distance_matrix.cpp b/example/Rips_complex/example_one_skeleton_rips_from_distance_matrix.cpp new file mode 100644 index 00000000..bbc3c755 --- /dev/null +++ b/example/Rips_complex/example_one_skeleton_rips_from_distance_matrix.cpp @@ -0,0 +1,58 @@ +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for std::numeric_limits + +int main() { + // Type definitions + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + using Distance_matrix = std::vector<std::vector<Filtration_value>>; + + // User defined distance matrix is: + // | 0 0.94 0.77 0.99 0.11 | + // | 0.94 0 0.26 0.99 0.39 | + // | 0.77 0.26 0 0.28 0.97 | + // | 0.99 0.99 0.28 0 0.30 | + // | 0.11 0.39 0.97 0.30 0 | + + Distance_matrix distances; + distances.push_back({}); + distances.push_back({0.94}); + distances.push_back({0.77, 0.26}); + distances.push_back({0.99, 0.99, 0.28}); + distances.push_back({0.11, 0.39, 0.97, 0.30}); + + // ---------------------------------------------------------------------------- + // Init of a Rips complex from points + // ---------------------------------------------------------------------------- + double threshold = 1.0; + Rips_complex rips_complex_from_points(distances, threshold); + + Simplex_tree stree; + rips_complex_from_points.create_complex(stree, 1); + // ---------------------------------------------------------------------------- + // Display information about the one skeleton Rips complex + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + + return 0; +} diff --git a/example/Rips_complex/example_one_skeleton_rips_from_points.cpp b/example/Rips_complex/example_one_skeleton_rips_from_points.cpp new file mode 100644 index 00000000..a1db8910 --- /dev/null +++ b/example/Rips_complex/example_one_skeleton_rips_from_points.cpp @@ -0,0 +1,52 @@ +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for std::numeric_limits + +int main() { + // Type definitions + using Point = std::vector<double>; + 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>; + + std::vector<Point> points; + points.push_back({1.0, 1.0}); + points.push_back({7.0, 0.0}); + points.push_back({4.0, 6.0}); + points.push_back({9.0, 6.0}); + points.push_back({0.0, 14.0}); + points.push_back({2.0, 19.0}); + points.push_back({9.0, 17.0}); + + // ---------------------------------------------------------------------------- + // Init of a Rips complex from points + // ---------------------------------------------------------------------------- + double threshold = 12.0; + Rips_complex rips_complex_from_points(points, threshold, Gudhi::Euclidean_distance()); + + Simplex_tree stree; + rips_complex_from_points.create_complex(stree, 1); + // ---------------------------------------------------------------------------- + // Display information about the one skeleton Rips complex + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + return 0; +} diff --git a/example/Rips_complex/example_rips_complex_from_csv_distance_matrix_file.cpp b/example/Rips_complex/example_rips_complex_from_csv_distance_matrix_file.cpp new file mode 100644 index 00000000..7ae8126f --- /dev/null +++ b/example/Rips_complex/example_rips_complex_from_csv_distance_matrix_file.cpp @@ -0,0 +1,72 @@ +#include <gudhi/Rips_complex.h> +// to construct Rips_complex from a csv file representing a distance matrix +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.csv threshold dim_max [ouput_file.txt]\n"; + std::cerr << " i.e.: " << progName << " ../../data/distance_matrix/full_square_distance_matrix.csv 1.0 3\n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 4) && (argc != 5)) usage(argc, (argv[0] - 1)); + + std::string csv_file_name(argv[1]); + double threshold = atof(argv[2]); + int dim_max = atoi(argv[3]); + + // Type definitions + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + using Distance_matrix = std::vector<std::vector<Filtration_value>>; + + // ---------------------------------------------------------------------------- + // Init of a Rips complex from a distance matrix in a csv file + // Default separator is ';' + // ---------------------------------------------------------------------------- + Distance_matrix distances = read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name); + Rips_complex rips_complex_from_file(distances, threshold); + + std::streambuf* streambufffer; + std::ofstream ouput_file_stream; + + if (argc == 5) { + ouput_file_stream.open(std::string(argv[4])); + streambufffer = ouput_file_stream.rdbuf(); + } else { + streambufffer = std::cout.rdbuf(); + } + + Simplex_tree stree; + rips_complex_from_file.create_complex(stree, dim_max); + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the Rips complex + // ---------------------------------------------------------------------------- + output_stream << "Rips complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + output_stream << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + output_stream << std::endl; + } + + ouput_file_stream.close(); + return 0; +} diff --git a/example/Rips_complex/example_rips_complex_from_off_file.cpp b/example/Rips_complex/example_rips_complex_from_off_file.cpp new file mode 100644 index 00000000..de2e4ea4 --- /dev/null +++ b/example/Rips_complex/example_rips_complex_from_off_file.cpp @@ -0,0 +1,71 @@ +#include <gudhi/Rips_complex.h> +// to construct Rips_complex from a OFF file of points +#include <gudhi/Points_off_io.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off threshold dim_max [ouput_file.txt]\n"; + std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0\n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 4) && (argc != 5)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + double threshold = atof(argv[2]); + int dim_max = atoi(argv[3]); + + // Type definitions + using Point = std::vector<float>; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + + // ---------------------------------------------------------------------------- + // Init of a Rips complex from an OFF file + // ---------------------------------------------------------------------------- + Gudhi::Points_off_reader<Point> off_reader(off_file_name); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + std::streambuf* streambufffer; + std::ofstream ouput_file_stream; + + if (argc == 5) { + ouput_file_stream.open(std::string(argv[4])); + streambufffer = ouput_file_stream.rdbuf(); + } else { + streambufffer = std::cout.rdbuf(); + } + + Simplex_tree stree; + rips_complex_from_file.create_complex(stree, dim_max); + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the Rips complex + // ---------------------------------------------------------------------------- + output_stream << "Rips complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + output_stream << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + output_stream << std::endl; + } + + ouput_file_stream.close(); + return 0; +} diff --git a/example/Rips_complex/full_skeleton_rips_for_doc.txt b/example/Rips_complex/full_skeleton_rips_for_doc.txt new file mode 100644 index 00000000..55de4ab8 --- /dev/null +++ b/example/Rips_complex/full_skeleton_rips_for_doc.txt @@ -0,0 +1,26 @@ +Rips complex is of dimension 3 - 24 simplices - 7 vertices. +Iterator on Rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [5] + ( 5 4 ) -> [5.38516] + ( 2 0 ) -> [5.83095] + ( 1 0 ) -> [6.08276] + ( 3 1 ) -> [6.32456] + ( 2 1 ) -> [6.7082] + ( 2 1 0 ) -> [6.7082] + ( 3 2 1 ) -> [6.7082] + ( 6 5 ) -> [7.28011] + ( 4 2 ) -> [8.94427] + ( 3 0 ) -> [9.43398] + ( 3 1 0 ) -> [9.43398] + ( 3 2 0 ) -> [9.43398] + ( 3 2 1 0 ) -> [9.43398] + ( 6 4 ) -> [9.48683] + ( 6 5 4 ) -> [9.48683] + ( 6 3 ) -> [11] diff --git a/example/Rips_complex/one_skeleton_rips_for_doc.txt b/example/Rips_complex/one_skeleton_rips_for_doc.txt new file mode 100644 index 00000000..706512a5 --- /dev/null +++ b/example/Rips_complex/one_skeleton_rips_for_doc.txt @@ -0,0 +1,20 @@ +Rips complex is of dimension 1 - 18 simplices - 7 vertices. +Iterator on Rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [5] + ( 5 4 ) -> [5.38516] + ( 2 0 ) -> [5.83095] + ( 1 0 ) -> [6.08276] + ( 3 1 ) -> [6.32456] + ( 2 1 ) -> [6.7082] + ( 6 5 ) -> [7.28011] + ( 4 2 ) -> [8.94427] + ( 3 0 ) -> [9.43398] + ( 6 4 ) -> [9.48683] + ( 6 3 ) -> [11] diff --git a/example/Simplex_tree/CMakeLists.txt b/example/Simplex_tree/CMakeLists.txt index 9314a805..b1ea98d4 100644 --- a/example/Simplex_tree/CMakeLists.txt +++ b/example/Simplex_tree/CMakeLists.txt @@ -1,29 +1,32 @@ cmake_minimum_required(VERSION 2.6) project(Simplex_tree_examples) -add_executable ( simplex_tree_from_cliques_of_graph simplex_tree_from_cliques_of_graph.cpp ) +add_executable ( Simplex_tree_example_from_cliques_of_graph simplex_tree_from_cliques_of_graph.cpp ) if (TBB_FOUND) - target_link_libraries(simplex_tree_from_cliques_of_graph ${TBB_LIBRARIES}) + target_link_libraries(Simplex_tree_example_from_cliques_of_graph ${TBB_LIBRARIES}) endif() -add_test(simplex_tree_from_cliques_of_graph_2 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_cliques_of_graph ${CMAKE_SOURCE_DIR}/data/points/Klein_bottle_complex.txt 2) -add_test(simplex_tree_from_cliques_of_graph_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_cliques_of_graph ${CMAKE_SOURCE_DIR}/data/points/Klein_bottle_complex.txt 3) +add_test(NAME Simplex_tree_example_from_cliques_of_graph_2 COMMAND $<TARGET_FILE:Simplex_tree_example_from_cliques_of_graph> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/Klein_bottle_complex.fsc" "2") +add_test(NAME Simplex_tree_example_from_cliques_of_graph_3 COMMAND $<TARGET_FILE:Simplex_tree_example_from_cliques_of_graph> + "${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/Klein_bottle_complex.fsc" "3") -add_executable ( simple_simplex_tree simple_simplex_tree.cpp ) +add_executable ( Simplex_tree_example_simple_simplex_tree simple_simplex_tree.cpp ) if (TBB_FOUND) - target_link_libraries(simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(Simplex_tree_example_simple_simplex_tree ${TBB_LIBRARIES}) endif() -add_test(simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/simple_simplex_tree) +add_test(NAME Simplex_tree_example_simple_simplex_tree COMMAND $<TARGET_FILE:Simplex_tree_example_simple_simplex_tree>) -add_executable ( mini_simplex_tree mini_simplex_tree.cpp ) -add_test(mini_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/mini_simplex_tree) +add_executable ( Simplex_tree_example_mini_simplex_tree mini_simplex_tree.cpp ) +add_test(NAME Simplex_tree_example_mini_simplex_tree COMMAND $<TARGET_FILE:Simplex_tree_example_mini_simplex_tree>) # An example with Simplex-tree using CGAL alpha_shapes_3 if(GMP_FOUND AND CGAL_FOUND) - add_executable ( simplex_tree_from_alpha_shapes_3 simplex_tree_from_alpha_shapes_3.cpp ) - target_link_libraries(simplex_tree_from_alpha_shapes_3 ${GMP_LIBRARIES} ${CGAL_LIBRARY} ${Boost_SYSTEM_LIBRARY}) + add_executable ( Simplex_tree_example_alpha_shapes_3_from_off example_alpha_shapes_3_simplex_tree_from_off_file.cpp ) + target_link_libraries(Simplex_tree_example_alpha_shapes_3_from_off ${GMP_LIBRARIES} ${CGAL_LIBRARY} ${Boost_SYSTEM_LIBRARY}) if (TBB_FOUND) - target_link_libraries(simplex_tree_from_alpha_shapes_3 ${TBB_LIBRARIES}) + target_link_libraries(Simplex_tree_example_alpha_shapes_3_from_off ${TBB_LIBRARIES}) endif() - add_test(simplex_tree_from_alpha_shapes_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_alpha_shapes_3 ${CMAKE_SOURCE_DIR}/data/points/bunny_5000) + add_test(NAME Simplex_tree_example_alpha_shapes_3_from_off COMMAND $<TARGET_FILE:Simplex_tree_example_alpha_shapes_3_from_off> + "${CMAKE_SOURCE_DIR}/data/points/bunny_5000.off") endif() diff --git a/example/Simplex_tree/README b/example/Simplex_tree/README index e37af790..a9498173 100644 --- a/example/Simplex_tree/README +++ b/example/Simplex_tree/README @@ -10,7 +10,7 @@ Example of use : *** Simple simplex tree construction -./simple_simplex_tree +./Simplex_tree_example_simple_simplex_tree ******************************************************************** EXAMPLE OF SIMPLE INSERTION @@ -52,7 +52,7 @@ EXAMPLE OF SIMPLE INSERTION *** Simplex tree construction with Z/2Z coefficients on weighted graph Klein bottle file: -./simplex_tree_from_cliques_of_graph ../../../data/points/Klein_bottle_complex.txt 2 +./Simplex_tree_example_from_cliques_of_graph ../../../data/points/Klein_bottle_complex.txt 2 Insert the 1-skeleton in the simplex tree in 0 s. Expand the simplex tree in 0 s. Information of the Simplex Tree: @@ -60,7 +60,7 @@ Information of the Simplex Tree: with Z/3Z coefficients: -./simplex_tree_from_cliques_of_graph ../../../data/points/Klein_bottle_complex.txt 3 +./Simplex_tree_example_from_cliques_of_graph ../../../data/points/Klein_bottle_complex.txt 3 Insert the 1-skeleton in the simplex tree in 0 s. Expand the simplex tree in 0 s. @@ -70,4 +70,4 @@ Information of the Simplex Tree: *** Simplex_tree computed and displayed from a 3D alpha complex: [ Requires CGAL, GMP and GMPXX to be installed] -./simplex_tree_from_alpha_shapes_3 ../../../data/points/bunny_5000 +./Simplex_tree_example_alpha_shapes_3_from_off ../../../data/points/bunny_5000 diff --git a/example/Simplex_tree/simplex_tree_from_alpha_shapes_3.cpp b/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp index 49d358ab..ff2eebcb 100644 --- a/example/Simplex_tree/simplex_tree_from_alpha_shapes_3.cpp +++ b/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp @@ -20,8 +20,9 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#include <gudhi/graph_simplicial_complex.h> #include <gudhi/Simplex_tree.h> +#include <gudhi/Points_3D_off_io.h> + #include <boost/variant.hpp> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> @@ -118,24 +119,21 @@ int main(int argc, char * const argv[]) { // program args management if (argc != 2) { std::cerr << "Usage: " << argv[0] - << " path_to_file_graph \n"; + << " path_to_off_file \n"; return 0; } // Read points from file - std::string filegraph = argv[1]; - std::list<Point> lp; - std::ifstream is(filegraph.c_str()); - int n; - is >> n; -#ifdef DEBUG_TRACES - std::cout << "Reading " << n << " points " << std::endl; -#endif // DEBUG_TRACES - Point p; - for (; n > 0; n--) { - is >> p; - lp.push_back(p); + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << argv[1] << std::endl; + return 0; } + // Retrieve the triangulation + std::vector<Point> lp = off_reader.get_point_cloud(); // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL); diff --git a/example/Simplex_tree/mini_simplex_tree.cpp b/example/Simplex_tree/mini_simplex_tree.cpp index 7e48aaaf..ad99df23 100644 --- a/example/Simplex_tree/mini_simplex_tree.cpp +++ b/example/Simplex_tree/mini_simplex_tree.cpp @@ -24,19 +24,18 @@ #include <iostream> #include <initializer_list> -using namespace Gudhi; - -struct MyOptions : Simplex_tree_options_full_featured { +struct MyOptions : Gudhi::Simplex_tree_options_full_featured { // Not doing persistence, so we don't need those static const bool store_key = false; static const bool store_filtration = false; // I have few vertices typedef short Vertex_handle; }; -typedef Simplex_tree<MyOptions> ST; + +using ST = Gudhi::Simplex_tree<MyOptions>; // Dictionary should be private, but for now this is the easiest way. -static_assert(sizeof(ST::Dictionary::value_type) < sizeof(Simplex_tree<>::Dictionary::value_type), +static_assert(sizeof(ST::Dictionary::value_type) < sizeof(Gudhi::Simplex_tree<>::Dictionary::value_type), "Not storing the filtration and key should save some space"); int main() { diff --git a/example/Simplex_tree/simple_simplex_tree.cpp b/example/Simplex_tree/simple_simplex_tree.cpp index 5146b906..60f9a35e 100644 --- a/example/Simplex_tree/simple_simplex_tree.cpp +++ b/example/Simplex_tree/simple_simplex_tree.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 * * 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 @@ -27,10 +27,11 @@ #include <utility> // for pair #include <vector> -using namespace Gudhi; - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using typeVectorVertex = std::vector< Vertex_handle >; +using typePairSimplexBool = std::pair< Simplex_tree::Simplex_handle, bool >; int main(int argc, char * const argv[]) { const Filtration_value FIRST_FILTRATION_VALUE = 0.1; @@ -42,7 +43,7 @@ int main(int argc, char * const argv[]) { std::cout << "********************************************************************" << std::endl; std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl; // Construct the Simplex Tree - Simplex_tree<> simplexTree; + Simplex_tree simplexTree; /* Simplex to be inserted: */ /* 1 */ @@ -212,7 +213,7 @@ int main(int argc, char * const argv[]) { // ------------------------------------------------------------------------------------------------------------------ // Find in the simplex_tree // ------------------------------------------------------------------------------------------------------------------ - Simplex_tree<>::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector); + Simplex_tree::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector); std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n"; if (simplexFound != simplexTree.null_simplex()) std::cout << "***+ YES IT IS!\n"; diff --git a/example/Simplex_tree/simplex_tree_from_cliques_of_graph.cpp b/example/Simplex_tree/simplex_tree_from_cliques_of_graph.cpp index 58085014..d1b8b2de 100644 --- a/example/Simplex_tree/simplex_tree_from_cliques_of_graph.cpp +++ b/example/Simplex_tree/simplex_tree_from_cliques_of_graph.cpp @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * 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 @@ -26,9 +26,16 @@ #include <iostream> #include <ctime> #include <string> +#include <utility> // for std::pair using namespace Gudhi; +typedef int Vertex_handle; +typedef double Filtration_value; +typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS, + boost::property < vertex_filtration_t, Filtration_value >, + boost::property < edge_filtration_t, Filtration_value > > Graph_t; + int main(int argc, char * const argv[]) { if (argc != 3) { std::cerr << "Usage: " << argv[0] @@ -43,7 +50,7 @@ int main(int argc, char * const argv[]) { Simplex_tree<> st; start = clock(); - auto g = read_graph(filegraph); + auto g = read_graph<Graph_t, Filtration_value, Vertex_handle>(filegraph); // insert the graph in the simplex tree as 1-skeleton st.insert_graph(g); end = clock(); diff --git a/example/Skeleton_blocker/CMakeLists.txt b/example/Skeleton_blocker/CMakeLists.txt index cc7f37f3..c887e408 100644 --- a/example/Skeleton_blocker/CMakeLists.txt +++ b/example/Skeleton_blocker/CMakeLists.txt @@ -1,12 +1,12 @@ cmake_minimum_required(VERSION 2.6) project(Skeleton_blocker_examples) -add_executable(SkeletonBlockerFromSimplices Skeleton_blocker_from_simplices.cpp) -add_executable(SkeletonBlockerIteration Skeleton_blocker_iteration.cpp) -add_executable(SkeletonBlockerLink Skeleton_blocker_link.cpp) +add_executable(Skeleton_blocker_example_from_simplices Skeleton_blocker_from_simplices.cpp) +add_executable(Skeleton_blocker_example_iteration Skeleton_blocker_iteration.cpp) +add_executable(Skeleton_blocker_example_link Skeleton_blocker_link.cpp) -target_link_libraries(SkeletonBlockerIteration ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) +target_link_libraries(Skeleton_blocker_example_iteration ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) -add_test(SkeletonBlockerFromSimplices ${CMAKE_CURRENT_BINARY_DIR}/SkeletonBlockerFromSimplices) -add_test(SkeletonBlockerIteration ${CMAKE_CURRENT_BINARY_DIR}/SkeletonBlockerIteration) -add_test(SkeletonBlockerLink ${CMAKE_CURRENT_BINARY_DIR}/SkeletonBlockerLink) +add_test(NAME Skeleton_blocker_example_from_simplices COMMAND $<TARGET_FILE:Skeleton_blocker_example_from_simplices>) +add_test(NAME Skeleton_blocker_example_iteration COMMAND $<TARGET_FILE:Skeleton_blocker_example_iteration>) +add_test(NAME Skeleton_blocker_example_link COMMAND $<TARGET_FILE:Skeleton_blocker_example_link>) diff --git a/example/Skeleton_blocker/Skeleton_blocker_from_simplices.cpp b/example/Skeleton_blocker/Skeleton_blocker_from_simplices.cpp index 171f35f2..076c6a53 100644 --- a/example/Skeleton_blocker/Skeleton_blocker_from_simplices.cpp +++ b/example/Skeleton_blocker/Skeleton_blocker_from_simplices.cpp @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * Copyright (C) 2014 INRIA * * 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 @@ -29,11 +29,8 @@ #include <sstream> #include <vector> -using namespace std; -using namespace Gudhi; -using namespace skeleton_blocker; - -typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_simple_traits Traits; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_complex<Traits> Complex; typedef Complex::Vertex_handle Vertex_handle; typedef Complex::Simplex Simplex; @@ -47,7 +44,7 @@ int main(int argc, char *argv[]) { simplices.push_back(Simplex(Vertex_handle(3), Vertex_handle(0), Vertex_handle(1))); // get complex from top faces - Complex complex(make_complex_from_top_faces<Complex>(simplices.begin(), simplices.end())); + Complex complex(Gudhi::skeleton_blocker::make_complex_from_top_faces<Complex>(simplices.begin(), simplices.end())); std::cout << "Simplices:" << std::endl; diff --git a/example/Skeleton_blocker/Skeleton_blocker_iteration.cpp b/example/Skeleton_blocker/Skeleton_blocker_iteration.cpp index 8d9d1a67..6a1bc480 100644 --- a/example/Skeleton_blocker/Skeleton_blocker_iteration.cpp +++ b/example/Skeleton_blocker/Skeleton_blocker_iteration.cpp @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * Copyright (C) 2014 INRIA * * 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 @@ -30,12 +30,8 @@ #include <fstream> #include <sstream> - -using namespace std; -using namespace Gudhi; -using namespace skeleton_blocker; - -typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_simple_traits Traits; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_complex<Traits> Complex; typedef Complex::Vertex_handle Vertex_handle; typedef Complex::Simplex Simplex; diff --git a/example/Skeleton_blocker/Skeleton_blocker_link.cpp b/example/Skeleton_blocker/Skeleton_blocker_link.cpp index 1f937170..58322038 100644 --- a/example/Skeleton_blocker/Skeleton_blocker_link.cpp +++ b/example/Skeleton_blocker/Skeleton_blocker_link.cpp @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (France) + * Copyright (C) 2014 INRIA * * 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 @@ -28,11 +28,8 @@ #include <fstream> #include <sstream> -using namespace std; -using namespace Gudhi; -using namespace skeleton_blocker; - -typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_simple_traits Traits; +typedef Gudhi::skeleton_blocker::Skeleton_blocker_complex<Traits> Complex; typedef Complex::Vertex_handle Vertex_handle; typedef Complex::Root_vertex_handle Root_vertex_handle; typedef Complex::Simplex Simplex; @@ -47,25 +44,26 @@ int main(int argc, char *argv[]) { Simplex tetrahedron(Vertex_handle(0), Vertex_handle(1), Vertex_handle(2), Vertex_handle(3)); complex.add_simplex(tetrahedron); - cout << "complex:" << complex.to_string() << endl; + std::cout << "complex:" << complex.to_string() << std::endl; // build the link of vertex 1, eg a triangle {0,2,3} auto link = complex.link(Vertex_handle(1)); - cout << "link:" << link.to_string() << endl; + std::cout << "link:" << link.to_string() << std::endl; // Internally link is a subcomplex of 'complex' and its vertices are stored in a vector. // They can be accessed via Vertex_handle(x) where x is an index of the vector. // In that example, link has three vertices and thus it contains only // Vertex_handle(0),Vertex_handle(1) and Vertex_handle(2) are). for (int i = 0; i < 5; ++i) - cout << "link.contains_vertex(Vertex_handle(" << i << ")):" << link.contains_vertex(Vertex_handle(i)) << endl; - cout << endl; + std::cout << "link.contains_vertex(Vertex_handle(" << i << ")):" << link.contains_vertex(Vertex_handle(i)) << + std::endl; + std::cout << std::endl; // To access to the initial vertices eg (0,1,2,3,4), Root_vertex_handle must be used. // For instance, to test if the link contains the vertex that was labeled i: for (int i = 0; i < 5; ++i) - cout << "link.contains_vertex(Root_vertex_handle(" << i << ")):" << - link.contains_vertex(Root_vertex_handle(i)) << endl; + std::cout << "link.contains_vertex(Root_vertex_handle(" << i << ")):" << + link.contains_vertex(Root_vertex_handle(i)) << std::endl; return EXIT_SUCCESS; } diff --git a/example/Spatial_searching/CMakeLists.txt b/example/Spatial_searching/CMakeLists.txt new file mode 100644 index 00000000..f4b9f3cb --- /dev/null +++ b/example/Spatial_searching/CMakeLists.txt @@ -0,0 +1,9 @@ +cmake_minimum_required(VERSION 2.6) +project(Spatial_searching_examples) + +if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable( Spatial_searching_example_spatial_searching example_spatial_searching.cpp ) + target_link_libraries(Spatial_searching_example_spatial_searching ${CGAL_LIBRARY}) + add_test(NAME Spatial_searching_example_spatial_searching + COMMAND $<TARGET_FILE:Spatial_searching_example_spatial_searching>) +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/example/Spatial_searching/example_spatial_searching.cpp b/example/Spatial_searching/example_spatial_searching.cpp new file mode 100644 index 00000000..14b324ae --- /dev/null +++ b/example/Spatial_searching/example_spatial_searching.cpp @@ -0,0 +1,52 @@ +#include <gudhi/Kd_tree_search.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> + +namespace gss = Gudhi::spatial_searching; + +int main(void) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef typename K::Point_d Point; + typedef std::vector<Point> Points; + + typedef gss::Kd_tree_search<K, Points> Points_ds; + + CGAL::Random rd; + + Points points; + for (int i = 0; i < 500; ++i) + points.push_back(Point(rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1))); + + Points_ds points_ds(points); + + // 10-nearest neighbor query + std::cout << "10 nearest neighbors from points[20]:\n"; + auto knn_range = points_ds.query_k_nearest_neighbors(points[20], 10, true); + for (auto const& nghb : knn_range) + std::cout << nghb.first << " (sq. dist. = " << nghb.second << ")\n"; + + // Incremental nearest neighbor query + std::cout << "Incremental nearest neighbors:\n"; + auto inn_range = points_ds.query_incremental_nearest_neighbors(points[45]); + // Get the neighbors in distance order until we hit the first point + for (auto ins_iterator = inn_range.begin(); ins_iterator->first != 0; ++ins_iterator) + std::cout << ins_iterator->first << " (sq. dist. = " << ins_iterator->second << ")\n"; + + // 10-farthest neighbor query + std::cout << "10 farthest neighbors from points[20]:\n"; + auto kfn_range = points_ds.query_k_farthest_neighbors(points[20], 10, true); + for (auto const& nghb : kfn_range) + std::cout << nghb.first << " (sq. dist. = " << nghb.second << ")\n"; + + // Incremental farthest neighbor query + std::cout << "Incremental farthest neighbors:\n"; + auto ifn_range = points_ds.query_incremental_farthest_neighbors(points[45]); + // Get the neighbors in distance reverse order until we hit the first point + for (auto ifs_iterator = ifn_range.begin(); ifs_iterator->first != 0; ++ifs_iterator) + std::cout << ifs_iterator->first << " (sq. dist. = " << ifs_iterator->second << ")\n"; + + return 0; +} diff --git a/example/Subsampling/CMakeLists.txt b/example/Subsampling/CMakeLists.txt new file mode 100644 index 00000000..71b8d2e8 --- /dev/null +++ b/example/Subsampling/CMakeLists.txt @@ -0,0 +1,17 @@ +cmake_minimum_required(VERSION 2.6) +project(Subsampling_examples) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable(Subsampling_example_pick_n_random_points example_pick_n_random_points.cpp) + add_executable(Subsampling_example_choose_n_farthest_points example_choose_n_farthest_points.cpp) + add_executable(Subsampling_example_custom_kernel example_custom_kernel.cpp) + add_executable(Subsampling_example_sparsify_point_set example_sparsify_point_set.cpp) + target_link_libraries(Subsampling_example_sparsify_point_set ${CGAL_LIBRARY}) + + add_test(NAME Subsampling_example_pick_n_random_points + COMMAND $<TARGET_FILE:Subsampling_example_pick_n_random_points>) + add_test(NAME Subsampling_example_choose_n_farthest_points + COMMAND $<TARGET_FILE:Subsampling_example_choose_n_farthest_points>) + add_test(NAME Subsampling_example_sparsify_point_set + COMMAND $<TARGET_FILE:Subsampling_example_sparsify_point_set>) +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/example/Subsampling/example_choose_n_farthest_points.cpp b/example/Subsampling/example_choose_n_farthest_points.cpp new file mode 100644 index 00000000..ebf631fc --- /dev/null +++ b/example/Subsampling/example_choose_n_farthest_points.cpp @@ -0,0 +1,29 @@ +#include <gudhi/choose_n_farthest_points.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + +int main(void) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef typename K::Point_d Point_d; + + CGAL::Random rd; + + std::vector<Point_d> points; + for (int i = 0; i < 500; ++i) + points.push_back(Point_d(rd.get_double(-1., 1), rd.get_double(-1., 1), + rd.get_double(-1., 1), rd.get_double(-1., 1))); + + K k; + std::vector<Point_d> results; + Gudhi::subsampling::choose_n_farthest_points(k, points, 100, + Gudhi::subsampling::random_starting_point, + std::back_inserter(results)); + std::cout << "Before sparsification: " << points.size() << " points.\n"; + std::cout << "After sparsification: " << results.size() << " points.\n"; + + return 0; +} diff --git a/example/Subsampling/example_custom_kernel.cpp b/example/Subsampling/example_custom_kernel.cpp new file mode 100644 index 00000000..2d42bdde --- /dev/null +++ b/example/Subsampling/example_custom_kernel.cpp @@ -0,0 +1,65 @@ +#include <gudhi/choose_n_farthest_points.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + + +/* The class Kernel contains a distance function defined on the set of points {0, 1, 2, 3} + * and computes a distance according to the matrix: + * 0 1 2 4 + * 1 0 4 2 + * 2 4 0 1 + * 4 2 1 0 + */ +class Kernel { + public: + typedef double FT; + typedef unsigned Point_d; + + // Class Squared_distance_d + class Squared_distance_d { + private: + std::vector<std::vector<FT>> matrix_; + + public: + Squared_distance_d() { + matrix_.push_back(std::vector<FT>({0, 1, 2, 4})); + matrix_.push_back(std::vector<FT>({1, 0, 4, 2})); + matrix_.push_back(std::vector<FT>({2, 4, 0, 1})); + matrix_.push_back(std::vector<FT>({4, 2, 1, 0})); + } + + FT operator()(Point_d p1, Point_d p2) { + return matrix_[p1][p2]; + } + }; + + // Constructor + Kernel() {} + + // Object of type Squared_distance_d + Squared_distance_d squared_distance_d_object() const { + return Squared_distance_d(); + } +}; + +int main(void) { + typedef Kernel K; + typedef typename K::Point_d Point_d; + + K k; + std::vector<Point_d> points = {0, 1, 2, 3}; + std::vector<Point_d> results; + + Gudhi::subsampling::choose_n_farthest_points(k, points, 2, + Gudhi::subsampling::random_starting_point, + std::back_inserter(results)); + std::cout << "Before sparsification: " << points.size() << " points.\n"; + std::cout << "After sparsification: " << results.size() << " points.\n"; + std::cout << "Result table: {" << results[0] << "," << results[1] << "}\n"; + + return 0; +} diff --git a/example/Subsampling/example_pick_n_random_points.cpp b/example/Subsampling/example_pick_n_random_points.cpp new file mode 100644 index 00000000..1e38e405 --- /dev/null +++ b/example/Subsampling/example_pick_n_random_points.cpp @@ -0,0 +1,27 @@ +#include <gudhi/pick_n_random_points.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + +int main(void) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef typename K::Point_d Point_d; + + CGAL::Random rd; + + std::vector<Point_d> points; + for (int i = 0; i < 500; ++i) + points.push_back(Point_d(rd.get_double(-1., 1), rd.get_double(-1., 1), + rd.get_double(-1., 1), rd.get_double(-1., 1))); + + K k; + std::vector<Point_d> results; + Gudhi::subsampling::pick_n_random_points(points, 100, std::back_inserter(results)); + std::cout << "Before sparsification: " << points.size() << " points.\n"; + std::cout << "After sparsification: " << results.size() << " points.\n"; + + return 0; +} diff --git a/example/Subsampling/example_sparsify_point_set.cpp b/example/Subsampling/example_sparsify_point_set.cpp new file mode 100644 index 00000000..b35a18d9 --- /dev/null +++ b/example/Subsampling/example_sparsify_point_set.cpp @@ -0,0 +1,27 @@ +#include <gudhi/sparsify_point_set.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> +#include <iterator> + +int main(void) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef typename K::Point_d Point_d; + + CGAL::Random rd; + + std::vector<Point_d> points; + for (int i = 0; i < 500; ++i) + points.push_back(Point_d(rd.get_double(-1., 1), rd.get_double(-1., 1), + rd.get_double(-1., 1), rd.get_double(-1., 1))); + + K k; + std::vector<Point_d> results; + Gudhi::subsampling::sparsify_point_set(k, points, 0.4, std::back_inserter(results)); + std::cout << "Before sparsification: " << points.size() << " points.\n"; + std::cout << "After sparsification: " << results.size() << " points.\n"; + + return 0; +} diff --git a/example/Tangential_complex/CMakeLists.txt b/example/Tangential_complex/CMakeLists.txt new file mode 100644 index 00000000..339d0581 --- /dev/null +++ b/example/Tangential_complex/CMakeLists.txt @@ -0,0 +1,18 @@ +cmake_minimum_required(VERSION 2.6) +project(Tangential_complex_examples) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable( Tangential_complex_example_basic example_basic.cpp ) + target_link_libraries(Tangential_complex_example_basic ${CGAL_LIBRARY} ${Boost_DATE_TIME_LIBRARY}) + add_executable( Tangential_complex_example_with_perturb example_with_perturb.cpp ) + target_link_libraries(Tangential_complex_example_with_perturb ${CGAL_LIBRARY} ${Boost_DATE_TIME_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Tangential_complex_example_basic ${TBB_LIBRARIES}) + target_link_libraries(Tangential_complex_example_with_perturb ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Tangential_complex_example_basic + COMMAND $<TARGET_FILE:Tangential_complex_example_basic>) + add_test(NAME Tangential_complex_example_with_perturb + COMMAND $<TARGET_FILE:Tangential_complex_example_with_perturb>) +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/example/Tangential_complex/example_basic.cpp b/example/Tangential_complex/example_basic.cpp new file mode 100644 index 00000000..4f2b859e --- /dev/null +++ b/example/Tangential_complex/example_basic.cpp @@ -0,0 +1,46 @@ +#include <gudhi/Tangential_complex.h> +#include <gudhi/sparsify_point_set.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <array> +#include <vector> + +namespace tc = Gudhi::tangential_complex; + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_d Point; +typedef Kernel::Vector_d Vector; +typedef tc::Tangential_complex< +Kernel, CGAL::Dynamic_dimension_tag, +CGAL::Parallel_tag> TC; + +int main(void) { + const int INTRINSIC_DIM = 2; + const int AMBIENT_DIM = 3; + const int NUM_POINTS = 1000; + + Kernel k; + + // Generate points on a 2-sphere + CGAL::Random_points_on_sphere_d<Point> generator(AMBIENT_DIM, 3.); + std::vector<Point> points; + points.reserve(NUM_POINTS); + for (int i = 0; i < NUM_POINTS; ++i) + points.push_back(*generator++); + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.compute_tangential_complex(); + + // Export the TC into a Simplex_tree + Gudhi::Simplex_tree<> stree; + tc.create_complex(stree); + + // Display stats about inconsistencies + tc.number_of_inconsistent_simplices(true); // verbose + + return 0; +} diff --git a/example/Tangential_complex/example_with_perturb.cpp b/example/Tangential_complex/example_with_perturb.cpp new file mode 100644 index 00000000..d0d877ea --- /dev/null +++ b/example/Tangential_complex/example_with_perturb.cpp @@ -0,0 +1,53 @@ +#include <gudhi/Tangential_complex.h> +#include <gudhi/sparsify_point_set.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <array> +#include <vector> + +namespace subsampl = Gudhi::subsampling; +namespace tc = Gudhi::tangential_complex; + +typedef CGAL::Epick_d<CGAL::Dimension_tag < 3 >> Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_d Point; +typedef Kernel::Vector_d Vector; +typedef tc::Tangential_complex< +Kernel, CGAL::Dimension_tag<2>, +CGAL::Parallel_tag> TC; + +int main(void) { + const int INTRINSIC_DIM = 2; + const int AMBIENT_DIM = 3; + const int NUM_POINTS = 50; + + Kernel k; + + // Generate points on a 2-sphere + CGAL::Random_points_on_sphere_d<Point> generator(AMBIENT_DIM, 3.); + std::vector<Point> points; + points.reserve(NUM_POINTS); + for (int i = 0; i < NUM_POINTS; ++i) + points.push_back(*generator++); + + // Sparsify the point set + std::vector<Point> sparsified_points; + subsampl::sparsify_point_set(k, points, 0.1 * 0.1, + std::back_inserter(sparsified_points)); + sparsified_points.swap(points); + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.compute_tangential_complex(); + + // Try to fix inconsistencies. Give it 10 seconds to succeed + tc.fix_inconsistencies_using_perturbation(0.05, 10); + + // Export the TC into a Simplex_tree + Gudhi::Simplex_tree<> stree; + tc.create_complex(stree); + + return 0; +} diff --git a/example/Witness_complex/CMakeLists.txt b/example/Witness_complex/CMakeLists.txt index 4d67e0d0..670651ce 100644 --- a/example/Witness_complex/CMakeLists.txt +++ b/example/Witness_complex/CMakeLists.txt @@ -1,16 +1,46 @@ cmake_minimum_required(VERSION 2.6) project(Witness_complex_examples) -# A simple example - add_executable( witness_complex_from_file witness_complex_from_file.cpp ) - add_test( witness_complex_from_bunny ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) - -if(CGAL_FOUND) - if (NOT CGAL_VERSION VERSION_LESS 4.6.0) - if (EIGEN3_FOUND) - add_executable ( witness_complex_sphere witness_complex_sphere.cpp ) - target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - add_test( witness_complex_sphere_10 ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere 10) - endif(EIGEN3_FOUND) - endif (NOT CGAL_VERSION VERSION_LESS 4.6.0) +add_executable ( Witness_complex_example_nearest_landmark_table example_nearest_landmark_table.cpp ) +target_link_libraries(Witness_complex_example_nearest_landmark_table ${Boost_SYSTEM_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(Witness_complex_example_nearest_landmark_table ${TBB_LIBRARIES}) endif() +add_test(NAME Witness_complex_example_nearest_landmark_table + COMMAND $<TARGET_FILE:Witness_complex_example_nearest_landmark_table>) + +# CGAL and Eigen3 are required for Euclidean version of Witness +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0) + add_executable( Witness_complex_example_off example_witness_complex_off.cpp ) + target_link_libraries(Witness_complex_example_off ${Boost_SYSTEM_LIBRARY}) + add_executable( Witness_complex_example_strong_off example_strong_witness_complex_off.cpp ) + target_link_libraries(Witness_complex_example_strong_off ${Boost_SYSTEM_LIBRARY}) + add_executable ( Witness_complex_example_sphere example_witness_complex_sphere.cpp ) + target_link_libraries(Witness_complex_example_sphere ${Boost_SYSTEM_LIBRARY}) + + add_executable ( Witness_complex_example_witness_persistence example_witness_complex_persistence.cpp ) + target_link_libraries(Witness_complex_example_witness_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + add_executable ( Witness_complex_example_strong_witness_persistence example_strong_witness_persistence.cpp ) + target_link_libraries(Witness_complex_example_strong_witness_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(Witness_complex_example_witness_persistence ${TBB_LIBRARIES}) + target_link_libraries(Witness_complex_example_strong_witness_persistence ${TBB_LIBRARIES}) + endif() + + add_test(NAME Witness_complex_example_off_test_torus + COMMAND $<TARGET_FILE:Witness_complex_example_off> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3") + add_test(NAME Witness_complex_example_strong_off_test_torus + COMMAND $<TARGET_FILE:Witness_complex_example_strong_off> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3") + add_test(NAME Witness_complex_example_test_sphere_10 + COMMAND $<TARGET_FILE:Witness_complex_example_sphere> "10") + add_test(NAME Witness_complex_example_test_torus_persistence + COMMAND $<TARGET_FILE:Witness_complex_example_witness_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5") + add_test(NAME Witness_complex_example_strong_test_torus_persistence + COMMAND $<TARGET_FILE:Witness_complex_example_strong_witness_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5") +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0) diff --git a/example/Witness_complex/example_nearest_landmark_table.cpp b/example/Witness_complex/example_nearest_landmark_table.cpp new file mode 100644 index 00000000..b8594212 --- /dev/null +++ b/example/Witness_complex/example_nearest_landmark_table.cpp @@ -0,0 +1,69 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA (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/>. + */ + +#define BOOST_PARAMETER_MAX_ARITY 12 + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <fstream> +#include <utility> +#include <string> +#include <vector> + +int main(int argc, char * const argv[]) { + using Nearest_landmark_range = std::vector<std::pair<std::size_t, double>>; + using Nearest_landmark_table = std::vector<Nearest_landmark_range>; + using Witness_complex = Gudhi::witness_complex::Witness_complex<Nearest_landmark_table>; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; + + Simplex_tree simplex_tree; + Nearest_landmark_table nlt; + + // Example contains 5 witnesses and 5 landmarks + Nearest_landmark_range w0 = {std::make_pair(0, 0), std::make_pair(1, 1), std::make_pair(2, 2), + std::make_pair(3, 3), std::make_pair(4, 4)}; nlt.push_back(w0); + Nearest_landmark_range w1 = {std::make_pair(1, 0), std::make_pair(2, 1), std::make_pair(3, 2), + std::make_pair(4, 3), std::make_pair(0, 4)}; nlt.push_back(w1); + Nearest_landmark_range w2 = {std::make_pair(2, 0), std::make_pair(3, 1), std::make_pair(4, 2), + std::make_pair(0, 3), std::make_pair(1, 4)}; nlt.push_back(w2); + Nearest_landmark_range w3 = {std::make_pair(3, 0), std::make_pair(4, 1), std::make_pair(0, 2), + std::make_pair(1, 3), std::make_pair(2, 4)}; nlt.push_back(w3); + Nearest_landmark_range w4 = {std::make_pair(4, 0), std::make_pair(0, 1), std::make_pair(1, 2), + std::make_pair(2, 3), std::make_pair(3, 4)}; nlt.push_back(w4); + + Witness_complex witness_complex(nlt); + witness_complex.create_complex(simplex_tree, 4.1); + + std::cout << "Number of simplices: " << simplex_tree.num_simplices() << std::endl; + + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(11); + + pcoh.compute_persistent_cohomology(-0.1); + pcoh.output_diagram(); +} diff --git a/example/Witness_complex/example_strong_witness_complex_off.cpp b/example/Witness_complex/example_strong_witness_complex_off.cpp new file mode 100644 index 00000000..0ee9ee90 --- /dev/null +++ b/example/Witness_complex/example_strong_witness_complex_off.cpp @@ -0,0 +1,79 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA (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/Simplex_tree.h> +#include <gudhi/Euclidean_strong_witness_complex.h> +#include <gudhi/pick_n_random_points.h> +#include <gudhi/Points_off_io.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = typename K::Point_d; +using Witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>; +using Point_vector = std::vector<Point_d>; + +int main(int argc, char * const argv[]) { + if (argc != 5) { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]); + double alpha2 = atof(argv[3]); + clock_t start, end; + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_vector point_vector, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Strong witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + point_vector = Point_vector(off_reader.get_point_cloud()); + + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n"; + + // Choose landmarks + Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + + // Compute witness complex + start = clock(); + Witness_complex witness_complex(landmarks, + point_vector); + + witness_complex.create_complex(simplex_tree, alpha2, lim_dim); + end = clock(); + std::cout << "Strong witness complex took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; +} diff --git a/example/Witness_complex/example_strong_witness_persistence.cpp b/example/Witness_complex/example_strong_witness_persistence.cpp new file mode 100644 index 00000000..f786fe7b --- /dev/null +++ b/example/Witness_complex/example_strong_witness_persistence.cpp @@ -0,0 +1,171 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA (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/Simplex_tree.h> +#include <gudhi/Euclidean_strong_witness_complex.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/pick_n_random_points.h> + +#include <boost/program_options.hpp> + +#include <CGAL/Epick_d.h> + +#include <string> +#include <vector> +#include <limits> // infinity + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = K::Point_d; + +using Point_vector = std::vector<Point_d>; +using Strong_witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>; +using SimplexTree = Gudhi::Simplex_tree<>; + +using Filtration_value = SimplexTree::Filtration_value; + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>; + +void program_options(int argc, char * argv[] + , int & nbL + , std::string & file_name + , std::string & filediag + , Filtration_value & max_squared_alpha + , int & p + , int & dim_max + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string file_name; + std::string filediag; + Filtration_value max_squared_alpha; + int p, nbL, lim_d; + Filtration_value min_persistence; + SimplexTree simplex_tree; + + program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence); + + // Extract the points from the file file_name + Point_vector witnesses, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + witnesses = Point_vector(off_reader.get_point_cloud()); + std::cout << "Successfully read " << witnesses.size() << " points.\n"; + std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n"; + + // Choose landmarks from witnesses + Gudhi::subsampling::pick_n_random_points(witnesses, nbL, std::back_inserter(landmarks)); + + // Compute witness complex + Strong_witness_complex strong_witness_complex(landmarks, + witnesses); + + strong_witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d); + + 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(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[] + , int & nbL + , std::string & file_name + , std::string & filediag + , Filtration_value & max_squared_alpha + , int & p + , int & dim_max + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&file_name), + "Name of file containing a point set in off format."); + + po::options_description visible("Allowed options", 100); + Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity(); + visible.add_options() + ("help,h", "produce help message") + ("landmarks,l", po::value<int>(&nbL), + "Number of landmarks to choose from the point cloud.") + ("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-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha), + "Maximal squared relaxation parameter.") + ("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)->default_value(0), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()), + "Maximal dimension of the strong witness complex we want to compute."); + + 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 Strong witness 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/example/Witness_complex/example_witness_complex_off.cpp b/example/Witness_complex/example_witness_complex_off.cpp new file mode 100644 index 00000000..b36dac0d --- /dev/null +++ b/example/Witness_complex/example_witness_complex_off.cpp @@ -0,0 +1,60 @@ +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Euclidean_witness_complex.h> +#include <gudhi/pick_n_random_points.h> +#include <gudhi/Points_off_io.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = K::Point_d; +using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<K>; +using Point_vector = std::vector< Point_d >; + +int main(int argc, char * const argv[]) { + if (argc != 5) { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]); + double alpha2 = atof(argv[3]); + clock_t start, end; + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_vector point_vector, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + point_vector = Point_vector(off_reader.get_point_cloud()); + + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n"; + + // Choose landmarks + Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + + // Compute witness complex + start = clock(); + Witness_complex witness_complex(landmarks, + point_vector); + + witness_complex.create_complex(simplex_tree, alpha2, lim_dim); + end = clock(); + std::cout << "Witness complex took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; +} diff --git a/example/Witness_complex/example_witness_complex_persistence.cpp b/example/Witness_complex/example_witness_complex_persistence.cpp new file mode 100644 index 00000000..a1146922 --- /dev/null +++ b/example/Witness_complex/example_witness_complex_persistence.cpp @@ -0,0 +1,171 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA (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/Simplex_tree.h> +#include <gudhi/Euclidean_witness_complex.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/pick_n_random_points.h> + +#include <boost/program_options.hpp> + +#include <CGAL/Epick_d.h> + +#include <string> +#include <vector> +#include <limits> // infinity + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = K::Point_d; + +using Point_vector = std::vector<Point_d>; +using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<K>; +using SimplexTree = Gudhi::Simplex_tree<>; + +using Filtration_value = SimplexTree::Filtration_value; + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>; + +void program_options(int argc, char * argv[] + , int & nbL + , std::string & file_name + , std::string & filediag + , Filtration_value & max_squared_alpha + , int & p + , int & dim_max + , Filtration_value & min_persistence); + +int main(int argc, char * argv[]) { + std::string file_name; + std::string filediag; + Filtration_value max_squared_alpha; + int p, nbL, lim_d; + Filtration_value min_persistence; + SimplexTree simplex_tree; + + program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence); + + // Extract the points from the file file_name + Point_vector witnesses, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + witnesses = Point_vector(off_reader.get_point_cloud()); + std::cout << "Successfully read " << witnesses.size() << " points.\n"; + std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n"; + + // Choose landmarks from witnesses + Gudhi::subsampling::pick_n_random_points(witnesses, nbL, std::back_inserter(landmarks)); + + // Compute witness complex + Witness_complex witness_complex(landmarks, + witnesses); + + witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d); + + 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(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[] + , int & nbL + , std::string & file_name + , std::string & filediag + , Filtration_value & max_squared_alpha + , int & p + , int & dim_max + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&file_name), + "Name of file containing a point set in off format."); + + Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity(); + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("landmarks,l", po::value<int>(&nbL), + "Number of landmarks to choose from the point cloud.") + ("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-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha), + "Maximal squared relaxation parameter.") + ("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)->default_value(0), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()), + "Maximal dimension of the weak witness complex we want to compute."); + + 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 Weak witness 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/example/Witness_complex/witness_complex_sphere.cpp b/example/Witness_complex/example_witness_complex_sphere.cpp index b26c9f36..124fd99b 100644 --- a/example/Witness_complex/witness_complex_sphere.cpp +++ b/example/Witness_complex/example_witness_complex_sphere.cpp @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2016 INRIA (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 @@ -19,17 +19,16 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#define BOOST_PARAMETER_MAX_ARITY 12 - -#include <sys/types.h> -#include <sys/stat.h> +#define BOOST_PARAMETER_MAX_ARITY 12 #include <gudhi/Simplex_tree.h> -#include <gudhi/Witness_complex.h> -#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/Euclidean_witness_complex.h> +#include <gudhi/pick_n_random_points.h> #include <gudhi/reader_utils.h> +#include <CGAL/Epick_d.h> + #include <iostream> #include <fstream> #include <ctime> @@ -51,6 +50,9 @@ void write_data(Data_range & data, std::string filename) { } int main(int argc, char * const argv[]) { + using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<Kernel>; + if (argc != 2) { std::cerr << "Usage: " << argv[0] << " number_of_landmarks \n"; @@ -58,27 +60,27 @@ int main(int argc, char * const argv[]) { } int number_of_landmarks = atoi(argv[1]); - clock_t start, end; - - // Construct the Simplex Tree - Gudhi::Simplex_tree<> simplex_tree; std::vector< std::pair<int, double> > l_time; - // Read the point file + // Generate points for (int nbP = 500; nbP < 10000; nbP += 500) { - Point_Vector point_vector; + clock_t start, end; + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + Point_Vector point_vector, landmarks; generate_points_sphere(point_vector, nbP, 4); std::cout << "Successfully generated " << point_vector.size() << " points.\n"; std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; // Choose landmarks start = clock(); - std::vector<std::vector< int > > knn; - Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, number_of_landmarks, knn); + Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks)); // Compute witness complex - Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); + Witness_complex witness_complex(landmarks, + point_vector); + witness_complex.create_complex(simplex_tree, 0); end = clock(); double time = static_cast<double>(end - start) / CLOCKS_PER_SEC; std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " diff --git a/example/Witness_complex/generators.h b/example/Witness_complex/generators.h index ac445261..7df43db5 100644 --- a/example/Witness_complex/generators.h +++ b/example/Witness_complex/generators.h @@ -25,17 +25,19 @@ #include <CGAL/Epick_d.h> #include <CGAL/point_generators_d.h> +#include <CGAL/Random.h> #include <fstream> #include <string> #include <vector> +#include <cmath> -typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; -typedef K::FT FT; -typedef K::Point_d Point_d; -typedef std::vector<Point_d> Point_Vector; -typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; -typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using FT = K::FT; +using Point_d = K::Point_d; +using Point_Vector = std::vector<Point_d>; +using Random_cube_iterator = CGAL::Random_points_in_cube_d<Point_d>; +using Random_point_iterator = CGAL::Random_points_in_ball_d<Point_d>; /** * \brief Rock age method of reading off file @@ -144,4 +146,21 @@ void generate_points_sphere(Point_Vector& W, int nbP, int dim) { W.push_back(*rp++); } +/** \brief Generate nbP points on a (flat) d-torus embedded in R^{2d} + * + */ +void generate_points_torus(Point_Vector& W, int nbP, int dim) { + CGAL::Random rand; + const double pi = std::acos(-1); + for (int i = 0; i < nbP; i++) { + std::vector<FT> point; + for (int j = 0; j < dim; j++) { + double alpha = rand.uniform_real(static_cast<double>(0), 2*pi); + point.push_back(sin(alpha)); + point.push_back(cos(alpha)); + } + W.push_back(Point_d(point)); + } +} + #endif // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ diff --git a/example/Witness_complex/witness_complex_from_file.cpp b/example/Witness_complex/witness_complex_from_file.cpp deleted file mode 100644 index 53207ad2..00000000 --- a/example/Witness_complex/witness_complex_from_file.cpp +++ /dev/null @@ -1,100 +0,0 @@ -/* 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): Siargey Kachanovich - * - * Copyright (C) 2015 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 <sys/types.h> -#include <sys/stat.h> - -#include <gudhi/Simplex_tree.h> -#include <gudhi/Witness_complex.h> -#include <gudhi/Landmark_choice_by_random_point.h> -#include <gudhi/reader_utils.h> - -#include <iostream> -#include <fstream> -#include <ctime> -#include <string> -#include <vector> - -typedef std::vector< Vertex_handle > typeVectorVertex; -typedef std::vector< std::vector <double> > Point_Vector; - -/** - * \brief Customized version of read_points - * which takes into account a possible nbP first line - * - */ -inline void -read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { - std::ifstream in_file(file_name.c_str(), std::ios::in); - if (!in_file.is_open()) { - std::cerr << "Unable to open file " << file_name << std::endl; - return; - } - std::string line; - double x; - while (getline(in_file, line)) { - std::vector< double > point; - std::istringstream iss(line); - while (iss >> x) { - point.push_back(x); - } - if (point.size() != 1) - points.push_back(point); - } - in_file.close(); -} - -int main(int argc, char * const argv[]) { - if (argc != 3) { - std::cerr << "Usage: " << argv[0] - << " path_to_point_file nbL \n"; - return 0; - } - - std::string file_name = argv[1]; - int nbL = atoi(argv[2]); - clock_t start, end; - - // Construct the Simplex Tree - Gudhi::Simplex_tree<> simplex_tree; - - // Read the point file - Point_Vector point_vector; - read_points_cust(file_name, point_vector); - std::cout << "Successfully read " << point_vector.size() << " points.\n"; - std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - - // Choose landmarks - start = clock(); - std::vector<std::vector< int > > knn; - Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, nbL, knn); - end = clock(); - std::cout << "Landmark choice for " << nbL << " landmarks took " - << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; - - // Compute witness complex - start = clock(); - Gudhi::witness_complex::witness_complex(knn, nbL, point_vector[0].size(), simplex_tree); - end = clock(); - std::cout << "Witness complex took " - << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; -} diff --git a/example/common/CMakeLists.txt b/example/common/CMakeLists.txt index 0da3dcc0..d5311b18 100644 --- a/example/common/CMakeLists.txt +++ b/example/common/CMakeLists.txt @@ -1,17 +1,22 @@ cmake_minimum_required(VERSION 2.6) project(Common_examples) -# need CGAL 4.7 +add_executable ( vector_double_off_reader example_vector_double_points_off_reader.cpp ) +target_link_libraries(vector_double_off_reader ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) +add_test(NAME Common_example_vector_double_off_reader COMMAND $<TARGET_FILE:vector_double_off_reader> + "${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off") + if(CGAL_FOUND) - add_executable ( cgal3Doffreader CGAL_3D_points_off_reader.cpp ) - target_link_libraries(cgal3Doffreader ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - add_test(cgal3Doffreader ${CMAKE_CURRENT_BINARY_DIR}/cgal3Doffreader ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) + add_executable ( cgal_3D_off_reader example_CGAL_3D_points_off_reader.cpp ) + target_link_libraries(cgal_3D_off_reader ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(NAME Common_example_vector_cgal_3D_off_reader COMMAND $<TARGET_FILE:cgal_3D_off_reader> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off") - if (NOT CGAL_VERSION VERSION_LESS 4.7.0) - if (EIGEN3_FOUND) - add_executable ( cgaloffreader CGAL_points_off_reader.cpp ) - target_link_libraries(cgaloffreader ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) - add_test(cgaloffreader ${CMAKE_CURRENT_BINARY_DIR}/cgaloffreader ${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off) - endif(EIGEN3_FOUND) - endif (NOT CGAL_VERSION VERSION_LESS 4.7.0) + # need CGAL 4.7and Eigen3 + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + add_executable ( cgal_off_reader example_CGAL_points_off_reader.cpp ) + target_link_libraries(cgal_off_reader ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(NAME Common_example_vector_cgal_off_reader COMMAND $<TARGET_FILE:cgal_off_reader> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off") + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) endif() diff --git a/example/common/CGAL_3D_points_off_reader.cpp b/example/common/example_CGAL_3D_points_off_reader.cpp index d48bb17d..665b7a29 100644 --- a/example/common/CGAL_3D_points_off_reader.cpp +++ b/example/common/example_CGAL_3D_points_off_reader.cpp @@ -32,7 +32,7 @@ int main(int argc, char **argv) { // Retrieve the triangulation std::vector<Point_3> point_cloud = off_reader.get_point_cloud(); - int n {0}; + int n {}; for (auto point : point_cloud) { ++n; std::cout << "Point[" << n << "] = (" << point[0] << ", " << point[1] << ", " << point[2] << ")\n"; diff --git a/example/common/CGAL_points_off_reader.cpp b/example/common/example_CGAL_points_off_reader.cpp index d1ca166d..8c6a6b54 100644 --- a/example/common/CGAL_points_off_reader.cpp +++ b/example/common/example_CGAL_points_off_reader.cpp @@ -34,10 +34,10 @@ int main(int argc, char **argv) { // Retrieve the triangulation std::vector<Point_d> point_cloud = off_reader.get_point_cloud(); - int n {0}; + int n {}; for (auto point : point_cloud) { std::cout << "Point[" << n << "] = "; - for (int i {0}; i < point.dimension(); i++) + for (std::size_t i {0}; i < point.size(); i++) std::cout << point[i] << " "; std::cout << "\n"; ++n; diff --git a/example/common/example_vector_double_points_off_reader.cpp b/example/common/example_vector_double_points_off_reader.cpp new file mode 100644 index 00000000..8aecb26e --- /dev/null +++ b/example/common/example_vector_double_points_off_reader.cpp @@ -0,0 +1,41 @@ +#include <gudhi/Points_off_io.h> + +#include <iostream> +#include <string> +#include <vector> + +using Point_d = std::vector<double>; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " inputFile.off" << std::endl; + exit(-1); +} + +int main(int argc, char **argv) { + if (argc != 2) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_off_reader<Point_d> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_d> point_cloud = off_reader.get_point_cloud(); + + int n {0}; + for (auto point : point_cloud) { + std::cout << "Point[" << n << "] = "; + for (std::size_t i {0}; i < point.size(); i++) + std::cout << point[i] << " "; + std::cout << "\n"; + ++n; + } + return 0; +} |