From 9899ae167f281d10b1684dfcd02c6838c5bf28df Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Fri, 2 Feb 2018 13:51:45 +0100 Subject: GUDHI 2.1.0 as released by upstream in a tarball. --- utilities/Alpha_complex/CMakeLists.txt | 65 ++++ utilities/Alpha_complex/alpha_complex_3d_helper.h | 74 ++++ .../Alpha_complex/alpha_complex_3d_persistence.cpp | 269 ++++++++++++++ .../Alpha_complex/alpha_complex_persistence.cpp | 116 ++++++ utilities/Alpha_complex/alphacomplex.md | 158 +++++++++ .../exact_alpha_complex_3d_persistence.cpp | 263 ++++++++++++++ .../periodic_alpha_complex_3d_persistence.cpp | 300 ++++++++++++++++ .../weighted_alpha_complex_3d_persistence.cpp | 314 +++++++++++++++++ ...ghted_periodic_alpha_complex_3d_persistence.cpp | 286 +++++++++++++++ utilities/Bitmap_cubical_complex/CMakeLists.txt | 29 ++ .../cubical_complex_persistence.cpp | 80 +++++ utilities/Bitmap_cubical_complex/cubicalcomplex.md | 29 ++ .../periodic_cubical_complex_persistence.cpp | 82 +++++ utilities/Bottleneck_distance/CMakeLists.txt | 16 + .../Bottleneck_distance/bottleneck_distance.cpp | 50 +++ .../Bottleneck_distance/bottleneckdistance.md | 18 + utilities/Nerve_GIC/CMakeLists.txt | 24 ++ utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py | 89 +++++ utilities/Nerve_GIC/Nerve.cpp | 96 +++++ utilities/Nerve_GIC/Nerve.txt | 63 ++++ utilities/Nerve_GIC/VoronoiGIC.cpp | 90 +++++ utilities/Nerve_GIC/covercomplex.md | 62 ++++ utilities/Nerve_GIC/km.py | 390 +++++++++++++++++++++ utilities/Nerve_GIC/km.py.COPYRIGHT | 26 ++ .../Persistence_representations/CMakeLists.txt | 53 +++ .../persistence_heat_maps/CMakeLists.txt | 15 + .../average_persistence_heat_maps.cpp | 63 ++++ .../compute_distance_of_persistence_heat_maps.cpp | 94 +++++ ...ute_scalar_product_of_persistence_heat_maps.cpp | 85 +++++ ...h_m_weighted_by_arctan_of_their_persistence.cpp | 81 +++++ ...te_p_h_m_weighted_by_distance_from_diagonal.cpp | 81 +++++ ...ate_p_h_m_weighted_by_squared_diag_distance.cpp | 83 +++++ .../create_persistence_heat_maps.cpp | 78 +++++ .../persistence_heat_maps/create_pssk.cpp | 79 +++++ .../plot_persistence_heat_map.cpp | 42 +++ .../persistence_intervals/CMakeLists.txt | 32 ++ ...te_birth_death_range_in_persistence_diagram.cpp | 68 ++++ .../compute_bottleneck_distance.cpp | 95 +++++ .../compute_number_of_dominant_intervals.cpp | 54 +++ .../plot_histogram_of_intervals_lengths.cpp | 77 ++++ .../plot_persistence_Betti_numbers.cpp | 87 +++++ .../plot_persistence_intervals.cpp | 53 +++ .../persistence_landscapes/CMakeLists.txt | 10 + .../persistence_landscapes/average_landscapes.cpp | 63 ++++ .../compute_distance_of_landscapes.cpp | 93 +++++ .../compute_scalar_product_of_landscapes.cpp | 84 +++++ .../persistence_landscapes/create_landscapes.cpp | 65 ++++ .../persistence_landscapes/plot_landscapes.cpp | 43 +++ .../persistence_landscapes_on_grid/CMakeLists.txt | 11 + .../average_landscapes_on_grid.cpp | 63 ++++ .../compute_distance_of_landscapes_on_grid.cpp | 93 +++++ ...ompute_scalar_product_of_landscapes_on_grid.cpp | 85 +++++ .../create_landscapes_on_grid.cpp | 79 +++++ .../plot_landscapes_on_grid.cpp | 44 +++ .../persistence_vectors/CMakeLists.txt | 10 + .../average_persistence_vectors.cpp | 65 ++++ .../compute_distance_of_persistence_vectors.cpp | 94 +++++ ...mpute_scalar_product_of_persistence_vectors.cpp | 86 +++++ .../create_persistence_vectors.cpp | 69 ++++ .../plot_persistence_vectors.cpp | 43 +++ utilities/Rips_complex/CMakeLists.txt | 21 ++ .../rips_distance_matrix_persistence.cpp | 133 +++++++ utilities/Rips_complex/rips_persistence.cpp | 135 +++++++ utilities/Rips_complex/ripscomplex.md | 49 +++ utilities/Witness_complex/CMakeLists.txt | 28 ++ .../Witness_complex/strong_witness_persistence.cpp | 156 +++++++++ .../Witness_complex/weak_witness_persistence.cpp | 156 +++++++++ utilities/Witness_complex/witnesscomplex.md | 66 ++++ utilities/common/README | 19 - utilities/common/pointsetgenerator.md | 33 ++ 70 files changed, 6086 insertions(+), 19 deletions(-) create mode 100644 utilities/Alpha_complex/CMakeLists.txt create mode 100644 utilities/Alpha_complex/alpha_complex_3d_helper.h create mode 100644 utilities/Alpha_complex/alpha_complex_3d_persistence.cpp create mode 100644 utilities/Alpha_complex/alpha_complex_persistence.cpp create mode 100644 utilities/Alpha_complex/alphacomplex.md create mode 100644 utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp create mode 100644 utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp create mode 100644 utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp create mode 100644 utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp create mode 100644 utilities/Bitmap_cubical_complex/CMakeLists.txt create mode 100644 utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp create mode 100644 utilities/Bitmap_cubical_complex/cubicalcomplex.md create mode 100644 utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp create mode 100644 utilities/Bottleneck_distance/CMakeLists.txt create mode 100644 utilities/Bottleneck_distance/bottleneck_distance.cpp create mode 100644 utilities/Bottleneck_distance/bottleneckdistance.md create mode 100644 utilities/Nerve_GIC/CMakeLists.txt create mode 100755 utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py create mode 100644 utilities/Nerve_GIC/Nerve.cpp create mode 100644 utilities/Nerve_GIC/Nerve.txt create mode 100644 utilities/Nerve_GIC/VoronoiGIC.cpp create mode 100644 utilities/Nerve_GIC/covercomplex.md create mode 100755 utilities/Nerve_GIC/km.py create mode 100644 utilities/Nerve_GIC/km.py.COPYRIGHT create mode 100644 utilities/Persistence_representations/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp create mode 100644 utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp create mode 100644 utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp create mode 100644 utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp create mode 100644 utilities/Persistence_representations/persistence_vectors/CMakeLists.txt create mode 100644 utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp create mode 100644 utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp create mode 100644 utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp create mode 100644 utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp create mode 100644 utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp create mode 100644 utilities/Rips_complex/CMakeLists.txt create mode 100644 utilities/Rips_complex/rips_distance_matrix_persistence.cpp create mode 100644 utilities/Rips_complex/rips_persistence.cpp create mode 100644 utilities/Rips_complex/ripscomplex.md create mode 100644 utilities/Witness_complex/CMakeLists.txt create mode 100644 utilities/Witness_complex/strong_witness_persistence.cpp create mode 100644 utilities/Witness_complex/weak_witness_persistence.cpp create mode 100644 utilities/Witness_complex/witnesscomplex.md delete mode 100644 utilities/common/README create mode 100644 utilities/common/pointsetgenerator.md (limited to 'utilities') diff --git a/utilities/Alpha_complex/CMakeLists.txt b/utilities/Alpha_complex/CMakeLists.txt new file mode 100644 index 00000000..a2dfac20 --- /dev/null +++ b/utilities/Alpha_complex/CMakeLists.txt @@ -0,0 +1,65 @@ +cmake_minimum_required(VERSION 2.6) +project(Alpha_complex_utilities) + +if(CGAL_FOUND) + add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) + target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_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(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") + add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") + add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "-p" "2" "-m" "0.45") + + install(TARGETS alpha_complex_3d_persistence DESTINATION bin) + install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin) + install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin) + + 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 ${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 ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") + add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $ + "${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" "-p" "2" "-m" "0") + + install(TARGETS alpha_complex_persistence DESTINATION bin) + install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin) + + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + + if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable(weighted_periodic_alpha_complex_3d_persistence weighted_periodic_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Alpha_complex_utilities_weigted_periodic_alpha_complex_3d COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "3" "1.0") + + install(TARGETS weighted_periodic_alpha_complex_3d_persistence DESTINATION bin) + + endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) + +endif(CGAL_FOUND) diff --git a/utilities/Alpha_complex/alpha_complex_3d_helper.h b/utilities/Alpha_complex/alpha_complex_3d_helper.h new file mode 100644 index 00000000..a59f0654 --- /dev/null +++ b/utilities/Alpha_complex/alpha_complex_3d_helper.h @@ -0,0 +1,74 @@ +/* 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 . + */ + +#ifndef ALPHA_COMPLEX_3D_HELPER_H_ +#define ALPHA_COMPLEX_3D_HELPER_H_ + +template +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 +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 +Vertex_list from_edge(const Edge_3& edg) { + Vertex_list the_list; + for (auto i : {edg.second, edg.third}) { +#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 +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/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..8ef5ffb2 --- /dev/null +++ b/utilities/Alpha_complex/alpha_complex_3d_persistence.cpp @@ -0,0 +1,269 @@ +/* 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 . + */ + +#include +#include +#include + +#if BOOST_VERSION >= 105400 +#include +#endif + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Alpha_shape_3 templates type definitions +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Vb = CGAL::Alpha_shape_vertex_base_3; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; + +// 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 >, + std::back_insert_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; + +#if BOOST_VERSION >= 105400 +using Vertex_list = boost::container::static_vector; +#else +using Vertex_list = std::vector; +#endif + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence); + + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader off_reader(off_file_points); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_points << std::endl; + exit(-1); + } + + // Retrieve the points + std::vector 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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); + count_facets++; + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); + count_edges++; + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + 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.push_back(vertex); + map_cgal_simplex_tree.emplace(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.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 + simplex_tree.insert_simplex(the_simplex, filtr); + GUDHI_CHECK(the_alpha_value_iterator != the_alpha_values.end(), "CGAL provided more simplices than values"); + ++the_alpha_value_iterator; + } + +#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() << " "; +#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); + + // 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; +} + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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 3D Alpha 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/utilities/Alpha_complex/alpha_complex_persistence.cpp b/utilities/Alpha_complex/alpha_complex_persistence.cpp new file mode 100644 index 00000000..2105220a --- /dev/null +++ b/utilities/Alpha_complex/alpha_complex_persistence.cpp @@ -0,0 +1,116 @@ +#include + +#include + +#include +#include +// to construct a simplex_tree from alpha complex +#include + +#include +#include +#include // 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, + Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, + Filtration_value &min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + Filtration_value alpha_square_max_value; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, coeff_field_characteristic, + min_persistence); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from an OFF file + // ---------------------------------------------------------------------------- + using Kernel = CGAL::Epick_d; + Gudhi::alpha_complex::Alpha_complex 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 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; +} + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, + Filtration_value &min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-alpha-square-value,r", po::value(&alpha_square_max_value) + ->default_value(std::numeric_limits::infinity()), + "Maximal alpha square value for the Alpha complex construction.")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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 an Alpha 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/utilities/Alpha_complex/alphacomplex.md b/utilities/Alpha_complex/alphacomplex.md new file mode 100644 index 00000000..aace85d3 --- /dev/null +++ b/utilities/Alpha_complex/alphacomplex.md @@ -0,0 +1,158 @@ + + +# Alpha complex # + + +## alpha_complex_persistence ## + +This program computes the persistent homology with coefficient field Z/pZ of the dD alpha complex built from a dD point cloud. +The output diagram contains one bar per line, written with the convention: + +``` + p dim birth death +``` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, +and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). + +**Usage** + +``` + alpha_complex_persistence [options] +``` + +where +`` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html). + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-r [ --max-alpha-square-value ]` (default = inf) Maximal alpha square value for the Alpha complex construction. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +**Example** + +``` + alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off +``` + +N.B.: + +* Filtration values are alpha square values. + + +## alpha_complex_3d_persistence ## +This program computes the persistent homology with coefficient field Z/pZ of the 3D alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention: + +``` +p dim birth death +``` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). + +**Usage** + +``` + alpha_complex_3d_persistence [options] +``` + +where `` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html). + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +**Example** + +``` +alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45 +``` + +N.B.: + +* `alpha_complex_3d_persistence` only accepts OFF files in dimension 3. +* Filtration values are alpha square values. + + +## exact_alpha_complex_3d_persistence ## + +Same as `alpha_complex_3d_persistence`, but using exact computation. +It is slower, but it is necessary when points are on a grid for instance. + + + +## weighted_alpha_complex_3d_persistence ## + +Same as `alpha_complex_3d_persistence`, but using weighted points. + +**Usage** + +``` + weighted_alpha_complex_3d_persistence [options] +``` + +where + +* `` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html). +* `` is the path to the file containing the weights of the points (one value per line). + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +**Example** + +``` + weighted_alpha_complex_3d_persistence ../../data/points/tore3D_300.off ../../data/points/tore3D_300.weights -p 2 -m 0.45 +``` + + +N.B.: + +* Weights values are explained on CGAL [Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0) +and [Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation. +* Filtration values are alpha square values. + + +## periodic_alpha_complex_3d_persistence ## +Same as `alpha_complex_3d_persistence`, but using periodic alpha shape 3d. +Refer to the [CGAL's 3D Periodic Triangulations User Manual](https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html) for more details. + +**Usage** + +``` + periodic_alpha_complex_3d_persistence [options] +``` + +where + +* `` is the path to the input point cloud in [nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html). +* `` is the path to the file describing the periodic domain. It must be in the format described +[here](/doc/latest/fileformats.html#FileFormatsIsoCuboid). + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals + + +**Example** + +``` +periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt -p 3 -m 1.0 +``` + +N.B.: + +* Cuboid file must be in the format described [here](/doc/latest/fileformats.html#FileFormatsIsoCuboid). +* Filtration values are alpha square values. diff --git a/utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/exact_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..cceac46e --- /dev/null +++ b/utilities/Alpha_complex/exact_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 . + */ + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#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; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; + +// 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 >, + std::back_insert_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::vector; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence); + + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader off_reader(off_file_points); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_points << std::endl; + exit(-1); + } + + // Retrieve the points + std::vector 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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); + count_facets++; + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); + count_edges++; + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + 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.push_back(vertex); + map_cgal_simplex_tree.emplace(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.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 + simplex_tree.insert_simplex(the_simplex, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + +#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() << " "; +#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); + + // 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; +} + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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 3D Alpha 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/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..188cf604 --- /dev/null +++ b/utilities/Alpha_complex/periodic_alpha_complex_3d_persistence.cpp @@ -0,0 +1,300 @@ +/* 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 + * Pawel Dlotko - 2017 - Swansea University, UK + * + * 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 . + */ + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Traits +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3; +// Vertex type +using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; +using Vb = CGAL::Triangulation_vertex_base_3; +using AsVb = CGAL::Alpha_shape_vertex_base_3; +// Cell type +using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; +using Cb = CGAL::Triangulation_cell_base_3; +using AsCb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; +using Point_3 = PK::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 >, + std::back_insert_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::vector; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file, + std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string cuboid_file; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, cuboid_file, output_file_diag, coeff_field_characteristic, + min_persistence); + + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader off_reader(off_file_points); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read OFF file " << off_file_points << std::endl; + exit(-1); + } + + // Read iso_cuboid_3 information from file + std::ifstream iso_cuboid_str(cuboid_file); + double x_min, y_min, z_min, x_max, y_max, z_max; + if (iso_cuboid_str.good()) { + iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max; + } else { + std::cerr << "Unable to read file " << cuboid_file << std::endl; + exit(-1); + } + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + if ((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min)) { + std::cerr << "The size of the cuboid in every directions is not the same." << std::endl; + exit(-1); + } + + // Retrieve the points + std::vector lp = off_reader.get_point_cloud(); + + // Define the periodic cube + P3DT3 pdt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + pdt.insert(lp.begin(), lp.end(), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + if (pdt.is_triangulation_in_1_sheet()) { + pdt.convert_to_1_sheeted_covering(); + } else { + std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl; + exit(-1); + } + std::cout << "Periodic Delaunay computed." << std::endl; + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL); + + // filtration with alpha values from alpha shape + std::vector the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); + count_facets++; + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); + count_edges++; + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + 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.push_back(vertex); + map_cgal_simplex_tree.emplace(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.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 + simplex_tree.insert_simplex(the_simplex, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + +#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() << " "; +#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); + + // 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; +} + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file, + std::string &output_file_diag, int &coeff_field_characteristic, + Filtration_value &min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd ")( + "cuboid-file", po::value(&cuboid_file), + "Name of file describing the periodic domain. Format is: min_hx min_hy min_hz\nmax_hx max_hy max_hz"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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); + pos.add("cuboid-file", 2); + + 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") || !vm.count("cuboid-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a periodic 3D Alpha 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 cuboid-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..93be8a05 --- /dev/null +++ b/utilities/Alpha_complex/weighted_alpha_complex_3d_persistence.cpp @@ -0,0 +1,314 @@ +/* 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 . + */ + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +// For CGAL < 4.11 +#if CGAL_VERSION_NR < 1041100000 +#include +#endif // CGAL_VERSION_NR < 1041100000 + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Alpha_shape_3 templates type definitions +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; + +// For CGAL < 4.11 +#if CGAL_VERSION_NR < 1041100000 +using Gt = CGAL::Regular_triangulation_euclidean_traits_3; +using Vb = CGAL::Alpha_shape_vertex_base_3; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Regular_triangulation_3; + +// From file type definition +using Point_3 = Gt::Bare_point; +using Weighted_point_3 = Gt::Weighted_point; + +// For CGAL >= 4.11 +#else // CGAL_VERSION_NR < 1041100000 +using Rvb = CGAL::Regular_triangulation_vertex_base_3; +using Vb = CGAL::Alpha_shape_vertex_base_3; +using Rcb = CGAL::Regular_triangulation_cell_base_3; +using Cb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Regular_triangulation_3; + +// From file type definition +using Point_3 = Triangulation_3::Bare_point; +using Weighted_point_3 = Triangulation_3::Weighted_point; +#endif // CGAL_VERSION_NR < 1041100000 + +using Alpha_shape_3 = CGAL::Alpha_shape_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 >, + std::back_insert_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::vector; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file, + std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string weight_file; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, weight_file, output_file_diag, coeff_field_characteristic, + min_persistence); + + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader off_reader(off_file_points); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read OFF file " << off_file_points << std::endl; + exit(-1); + } + + // Retrieve the points + std::vector lp = off_reader.get_point_cloud(); + + // Read weights information from file + std::ifstream weights_ifstr(weight_file); + std::vector wp; + if (weights_ifstr.good()) { + double weight = 0.0; + std::size_t index = 0; + wp.reserve(lp.size()); + // 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 " << weight_file << std::endl; + exit(-1); + } + } else { + std::cerr << "Unable to read weights file " << weight_file << std::endl; + exit(-1); + } + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. + Alpha_shape_3 as(wp.begin(), wp.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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); + count_facets++; + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); + count_edges++; + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + 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.push_back(vertex); + map_cgal_simplex_tree.emplace(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.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 + simplex_tree.insert_simplex(the_simplex, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + +#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() << " "; +#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); + + // 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; +} + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file, + std::string &output_file_diag, int &coeff_field_characteristic, + Filtration_value &min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd ")( + "weight-file", po::value(&weight_file), + "Name of file containing a point weights. Format is one weigt per line: W1\n...\nWn "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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); + pos.add("weight-file", 2); + + 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") || !vm.count("weight-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a weighted 3D Alpha 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 weight-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp b/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..5321bb0a --- /dev/null +++ b/utilities/Alpha_complex/weighted_periodic_alpha_complex_3d_persistence.cpp @@ -0,0 +1,286 @@ +/* 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 + * Pawel Dlotko - 2017 - Swansea University, UK + * + * 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 . + */ + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Traits +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using PK = CGAL::Periodic_3_regular_triangulation_traits_3; + +// Vertex type +using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; +using Vb = CGAL::Regular_triangulation_vertex_base_3; +using AsVb = CGAL::Alpha_shape_vertex_base_3; +// Cell type +using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; +using Cb = CGAL::Regular_triangulation_cell_base_3; +using AsCb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using P3RT3 = CGAL::Periodic_3_regular_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; + +using Point_3 = P3RT3::Bare_point; +using Weighted_point_3 = P3RT3::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 >, + std::back_insert_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::vector; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; + +void usage(const std::string& progName) { + std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_weight_file path_to_the_cuboid_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 != 6) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[4]); + Filtration_value min_persistence = strtof(argv[5], 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 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 points + std::vector lp = off_reader.get_point_cloud(); + + // Read iso_cuboid_3 information from file + std::ifstream iso_cuboid_str(argv[3]); + double x_min, y_min, z_min, x_max, y_max, z_max; + if (iso_cuboid_str.is_open()) { + if (!(iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max)) { + std::cerr << argv[3] << " - Bad file format." << std::endl; + usage(argv[0]); + } + + } else { + std::cerr << "Unable to read file " << argv[3] << std::endl; + usage(argv[0]); + } + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + if ((x_max - x_min != y_max - y_min) || (x_max - x_min != z_max - z_min) || (z_max - z_min != y_max - y_min)) { + std::cerr << "The size of the cuboid in every directions is not the same." << std::endl; + exit(-1); + } + + double maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min); + + // Read weights information from file + std::ifstream weights_ifstr(argv[2]); + std::vector wp; + if (weights_ifstr.is_open()) { + double weight = 0.0; + std::size_t index = 0; + wp.reserve(lp.size()); + // Attempt read the weight in a double format, return false if it fails + while ((weights_ifstr >> weight) && (index < lp.size())) { + if ((weight >= maximal_possible_weight) || (weight < 0)) { + std::cerr << "At line " << (index + 1) << ", the weight (" << weight + << ") is negative or more than or equal to maximal possible weight (" << maximal_possible_weight + << ") = 1/64*cuboid length squared, which is not an acceptable input." << std::endl; + exit(-1); + } + + 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]); + } + + // Define the periodic cube + P3RT3 prt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + prt.insert(wp.begin(), wp.end(), true); + // As prt won't be modified anymore switch to 1-sheeted cover if possible + if (prt.is_triangulation_in_1_sheet()) { + prt.convert_to_1_sheeted_covering(); + } else { + std::cerr << "ERROR: we were not able to construct a triangulation within a single periodic domain." << std::endl; + exit(-1); + } + std::cout << "Weighted Periodic Delaunay computed." << std::endl; + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + Alpha_shape_3 as(prt, 0, Alpha_shape_3::GENERAL); + + // filtration with alpha values from alpha shape + std::vector the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::iterator the_alpha_value_iterator = the_alpha_values.begin(); + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle* cell = CGAL::object_cast(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + } else if (const Facet* facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_facet(*facet); + count_facets++; + } else if (const Edge_3* edge = CGAL::object_cast(&object_iterator)) { + vertex_list = from_edge(*edge); + count_edges++; + } else if (const Vertex_handle* vertex = CGAL::object_cast(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + 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.push_back(vertex); + map_cgal_simplex_tree.emplace(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.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 + simplex_tree.insert_simplex(the_simplex, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + +#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() << " "; +#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/utilities/Bitmap_cubical_complex/CMakeLists.txt b/utilities/Bitmap_cubical_complex/CMakeLists.txt new file mode 100644 index 00000000..676a730a --- /dev/null +++ b/utilities/Bitmap_cubical_complex/CMakeLists.txt @@ -0,0 +1,29 @@ +cmake_minimum_required(VERSION 2.6) +project(Bitmap_cubical_complex_utilities) + +add_executable ( cubical_complex_persistence cubical_complex_persistence.cpp ) +if (TBB_FOUND) + target_link_libraries(cubical_complex_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Bitmap_cubical_complex_utility_persistence_one_sphere COMMAND $ + "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt") + +add_test(NAME Bitmap_cubical_complex_utility_persistence_two_sphere COMMAND $ + "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt") + +add_executable ( periodic_cubical_complex_persistence periodic_cubical_complex_persistence.cpp ) +if (TBB_FOUND) + target_link_libraries(periodic_cubical_complex_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Bitmap_cubical_complex_utility_periodic_boundary_conditions_2d_torus + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt") + +add_test(NAME Bitmap_cubical_complex_utility_periodic_boundary_conditions_3d_torus + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt") + +install(TARGETS cubical_complex_persistence DESTINATION bin) +install(TARGETS periodic_cubical_complex_persistence DESTINATION bin) diff --git a/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp b/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp new file mode 100644 index 00000000..9d1bc08c --- /dev/null +++ b/utilities/Bitmap_cubical_complex/cubical_complex_persistence.cpp @@ -0,0 +1,80 @@ +/* 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 + * + * Copyright (C) 2015 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 . + */ + +#include +#include +#include + +// standard stuff +#include +#include +#include +#include + +int main(int argc, char** argv) { + std::cout + << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " + << "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a" + << "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let " + << "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are " + << "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See " + << "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n" + << std::endl; + + if (argc != 2) { + std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " + << "the input. The program will now terminate.\n"; + return 1; + } + + typedef Gudhi::cubical_complex::Bitmap_cubical_complex_base Bitmap_cubical_complex_base; + typedef Gudhi::cubical_complex::Bitmap_cubical_complex Bitmap_cubical_complex; + typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp; + typedef Gudhi::persistent_cohomology::Persistent_cohomology Persistent_cohomology; + + Bitmap_cubical_complex b(argv[1]); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(b); + int p = 11; + double min_persistence = 0; + + pcoh.init_coefficients(p); // initializes the coefficient field for homology + pcoh.compute_persistent_cohomology(min_persistence); + + std::string output_file_name(argv[1]); + output_file_name += "_persistence"; + + std::size_t last_in_path = output_file_name.find_last_of("/\\"); + + if (last_in_path != std::string::npos) { + output_file_name = output_file_name.substr(last_in_path + 1); + } + + std::ofstream out(output_file_name.c_str()); + pcoh.output_diagram(out); + out.close(); + + std::cout << "Result in file: " << output_file_name << "\n"; + + return 0; +} diff --git a/utilities/Bitmap_cubical_complex/cubicalcomplex.md b/utilities/Bitmap_cubical_complex/cubicalcomplex.md new file mode 100644 index 00000000..6e1b2578 --- /dev/null +++ b/utilities/Bitmap_cubical_complex/cubicalcomplex.md @@ -0,0 +1,29 @@ + + +# Cubical complex# + +## cubical_complex_persistence ## +This program computes persistent homology, by using the Bitmap_cubical_complex class, of cubical complexes provided in text files in Perseus style. +See [here](/doc/latest/fileformats.html#FileFormatsPerseus) for a description of the file format. + +**Example** + +``` + cubical_complex_persistence data/bitmap/CubicalTwoSphere.txt +``` + +* Creates a Cubical Complex from the Perseus style file `CubicalTwoSphere.txt`, +computes Persistence cohomology from it and writes the results in a persistence file `CubicalTwoSphere.txt_persistence`. + +## periodic_cubical_complex_persistence ## + +Same as above, but with periodic boundary conditions. + +**Example** + +``` + periodic_cubical_complex_persistence data/bitmap/3d_torus.txt +``` + +* Creates a Periodical Cubical Complex from the Perseus style file `3d_torus.txt`, +computes Persistence cohomology from it and writes the results in a persistence file `3d_torus.txt_persistence`. diff --git a/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp b/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp new file mode 100644 index 00000000..c812cb3a --- /dev/null +++ b/utilities/Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp @@ -0,0 +1,82 @@ +/* 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 + * + * Copyright (C) 2015 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 . + */ + +#include +#include +#include +#include + +// standard stuff +#include +#include +#include +#include + +int main(int argc, char** argv) { + std::cout + << "This program computes persistent homology, by using " + << "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in " + << "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 " + << "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers " + << "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We " + << "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for" + << " example.\n" + << std::endl; + + if (argc != 2) { + std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " + << "the input. The program will now terminate.\n"; + return 1; + } + + typedef Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base Bitmap_base; + typedef Gudhi::cubical_complex::Bitmap_cubical_complex Bitmap_cubical_complex; + + Bitmap_cubical_complex b(argv[1]); + + typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp; + typedef Gudhi::persistent_cohomology::Persistent_cohomology Persistent_cohomology; + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(b, true); + + int p = 11; + double min_persistence = 0; + pcoh.init_coefficients(p); // initializes the coefficient field for homology + pcoh.compute_persistent_cohomology(min_persistence); + + std::string output_file_name(argv[1]); + output_file_name += "_persistence"; + + std::size_t last_in_path = output_file_name.find_last_of("/\\"); + + if (last_in_path != std::string::npos) { + output_file_name = output_file_name.substr(last_in_path + 1); + } + + std::ofstream out(output_file_name.c_str()); + pcoh.output_diagram(out); + out.close(); + + std::cout << "Result in file: " << output_file_name << "\n"; + + return 0; +} diff --git a/utilities/Bottleneck_distance/CMakeLists.txt b/utilities/Bottleneck_distance/CMakeLists.txt new file mode 100644 index 00000000..d19e3b1c --- /dev/null +++ b/utilities/Bottleneck_distance/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_utilities) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable (bottleneck_distance bottleneck_distance.cpp) + if (TBB_FOUND) + target_link_libraries(bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers") + + install(TARGETS bottleneck_distance DESTINATION bin) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/utilities/Bottleneck_distance/bottleneck_distance.cpp b/utilities/Bottleneck_distance/bottleneck_distance.cpp new file mode 100644 index 00000000..9dd52b31 --- /dev/null +++ b/utilities/Bottleneck_distance/bottleneck_distance.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 . + */ + +#include +#include +#include +#include +#include // for pair +#include +#include // for numeric_limits + +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 the bottleneck" << + " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << + " aware this version is exact but expensive. The program will now terminate \n"; + return -1; + } + std::vector> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); + std::vector> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); + + double tolerance = std::numeric_limits::min(); + 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; + + return 0; +} diff --git a/utilities/Bottleneck_distance/bottleneckdistance.md b/utilities/Bottleneck_distance/bottleneckdistance.md new file mode 100644 index 00000000..526f5822 --- /dev/null +++ b/utilities/Bottleneck_distance/bottleneckdistance.md @@ -0,0 +1,18 @@ + + +# Bottleneck distance # + +## bottleneck_read_file_example ## + +This program computes the Bottleneck distance between two persistence diagram files. + +**Usage** + +``` + bottleneck_read_file_example [] +``` + +where + +* `` and `` must be in the format described [here](/doc/latest/fileformats.html#FileFormatsPers). +* `` is an error bound on the bottleneck distance (set by default to the smallest positive double value). diff --git a/utilities/Nerve_GIC/CMakeLists.txt b/utilities/Nerve_GIC/CMakeLists.txt new file mode 100644 index 00000000..7762c8a0 --- /dev/null +++ b/utilities/Nerve_GIC/CMakeLists.txt @@ -0,0 +1,24 @@ +cmake_minimum_required(VERSION 2.6) +project(Nerve_GIC_examples) + +if (NOT CGAL_VERSION VERSION_LESS 4.8.1) + + add_executable ( Nerve Nerve.cpp ) + add_executable ( VoronoiGIC VoronoiGIC.cpp ) + + if (TBB_FOUND) + target_link_libraries(Nerve ${TBB_LIBRARIES}) + target_link_libraries(VoronoiGIC ${TBB_LIBRARIES}) + endif() + + file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + # Copy files for not to pollute sources when testing + file(COPY "${CMAKE_SOURCE_DIR}/data/points/human.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(NAME Nerve_GIC_utilities_nerve COMMAND $ + "human.off" "2" "10" "0.3") + + add_test(NAME Nerve_GIC_utilities_VoronoiGIC COMMAND $ + "human.off" "100") + +endif (NOT CGAL_VERSION VERSION_LESS 4.8.1) diff --git a/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py b/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py new file mode 100755 index 00000000..c811f610 --- /dev/null +++ b/utilities/Nerve_GIC/KeplerMapperVisuFromTxtFile.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python + +import km +import numpy as np +from collections import defaultdict +import argparse + +"""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): Mathieu Carriere + + 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 . +""" + +__author__ = "Mathieu Carriere" +__copyright__ = "Copyright (C) 2017 INRIA" +__license__ = "GPL v3" + +parser = argparse.ArgumentParser(description='Creates an html Keppler Mapper ' + 'file to visualize a SC.txt file.', + epilog='Example: ' + './KeplerMapperVisuFromTxtFile.py ' + '-f ../../data/points/human.off_sc.txt' + '- Constructs an human.off_sc.html file.') +parser.add_argument("-f", "--file", type=str, required=True) + +args = parser.parse_args() + +with open(args.file, 'r') as f: + network = {} + mapper = km.KeplerMapper(verbose=0) + data = np.zeros((3,3)) + projected_data = mapper.fit_transform( data, projection="sum", scaler=None ) + + nodes = defaultdict(list) + links = defaultdict(list) + custom = defaultdict(list) + + dat = f.readline() + lens = f.readline() + color = f.readline(); + param = [float(i) for i in f.readline().split(" ")] + + nums = [int(i) for i in f.readline().split(" ")] + num_nodes = nums[0] + num_edges = nums[1] + + for i in range(0,num_nodes): + point = [float(j) for j in f.readline().split(" ")] + nodes[ str(int(point[0])) ] = [ int(point[0]), point[1], int(point[2]) ] + links[ str(int(point[0])) ] = [] + custom[ int(point[0]) ] = point[1] + + m = min([custom[i] for i in range(0,num_nodes)]) + M = max([custom[i] for i in range(0,num_nodes)]) + + for i in range(0,num_edges): + edge = [int(j) for j in f.readline().split(" ")] + links[ str(edge[0]) ].append( str(edge[1]) ) + links[ str(edge[1]) ].append( str(edge[0]) ) + + network["nodes"] = nodes + network["links"] = links + network["meta"] = lens + + html_output_filename = args.file.rsplit('.', 1)[0] + '.html' + mapper.visualize(network, color_function = color, path_html=html_output_filename, title=dat, + graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0, + height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M) + message = repr(html_output_filename) + " is generated. You can now use your favorite web browser to visualize it." + print(message) + + + f.close() diff --git a/utilities/Nerve_GIC/Nerve.cpp b/utilities/Nerve_GIC/Nerve.cpp new file mode 100644 index 00000000..aefc3874 --- /dev/null +++ b/utilities/Nerve_GIC/Nerve.cpp @@ -0,0 +1,96 @@ +/* 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): Mathieu Carrière + * + * 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 . + */ + +#include + +#include +#include + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 5) && (argc != 6)) usage(argc, argv[0]); + + using Point = std::vector; + + std::string off_file_name(argv[1]); + int coord = atoi(argv[2]); + int resolution = atoi(argv[3]); + double gain = atof(argv[4]); + bool verb = 0; + if (argc == 6) verb = 1; + + // -------------------------------- + // Init of a Nerve from an OFF file + // -------------------------------- + + Gudhi::cover_complex::Cover_complex SC; + SC.set_verbose(verb); + + bool check = SC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + SC.set_type("Nerve"); + + SC.set_color_from_coordinate(coord); + SC.set_function_from_coordinate(coord); + + SC.set_graph_from_OFF(); + SC.set_resolution_with_interval_number(resolution); + SC.set_gain(gain); + SC.set_cover_from_function(); + + SC.find_simplices(); + + SC.write_info(); + + Gudhi::Simplex_tree<> stree; + SC.create_complex(stree); + SC.compute_PD(); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if (verb) { + std::cout << "Nerve is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Nerve simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + return 0; +} diff --git a/utilities/Nerve_GIC/Nerve.txt b/utilities/Nerve_GIC/Nerve.txt new file mode 100644 index 00000000..839ff45e --- /dev/null +++ b/utilities/Nerve_GIC/Nerve.txt @@ -0,0 +1,63 @@ +Min function value = -0.979672 and Max function value = 0.816414 +Interval 0 = [-0.979672, -0.761576] +Interval 1 = [-0.838551, -0.581967] +Interval 2 = [-0.658942, -0.402359] +Interval 3 = [-0.479334, -0.22275] +Interval 4 = [-0.299725, -0.0431415] +Interval 5 = [-0.120117, 0.136467] +Interval 6 = [0.059492, 0.316076] +Interval 7 = [0.239101, 0.495684] +Interval 8 = [0.418709, 0.675293] +Interval 9 = [0.598318, 0.816414] +Computing preimages... +Computing connected components... +.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox. +5 interval(s) in dimension 0: + [-0.909111, 0.00817529] + [-0.171433, 0.367392] + [-0.171433, 0.367392] + [-0.909111, 0.745853] +0 interval(s) in dimension 1: +Nerve is of dimension 1 - 41 simplices - 21 vertices. +Iterator on Nerve simplices +1 +0 +4 +4 0 +2 +2 1 +8 +8 2 +5 +5 4 +9 +9 8 +13 +13 5 +14 +14 9 +19 +19 13 +25 +32 +20 +32 20 +33 +33 25 +26 +26 14 +26 19 +42 +42 26 +34 +34 33 +27 +27 20 +35 +35 27 +35 34 +42 35 +44 +44 35 +54 +54 44 \ No newline at end of file diff --git a/utilities/Nerve_GIC/VoronoiGIC.cpp b/utilities/Nerve_GIC/VoronoiGIC.cpp new file mode 100644 index 00000000..54bb871e --- /dev/null +++ b/utilities/Nerve_GIC/VoronoiGIC.cpp @@ -0,0 +1,90 @@ +/* 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): Mathieu Carrière + * + * 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 . + */ + +#include + +#include +#include + +void usage(int nbArgs, char *const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off N [-v] \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 100 -v \n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, argv[0]); + + using Point = std::vector; + + std::string off_file_name(argv[1]); + int m = atoi(argv[2]); + bool verb = 0; + if (argc == 4) verb = 1; + + // ---------------------------------------------------------------------------- + // Init of a graph induced complex from an OFF file + // ---------------------------------------------------------------------------- + + Gudhi::cover_complex::Cover_complex GIC; + GIC.set_verbose(verb); + + bool check = GIC.read_point_cloud(off_file_name); + + if (!check) { + std::cout << "Incorrect OFF file." << std::endl; + } else { + GIC.set_type("GIC"); + + GIC.set_color_from_coordinate(); + + GIC.set_graph_from_OFF(); + GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), m); + + GIC.find_simplices(); + + GIC.plot_OFF(); + + Gudhi::Simplex_tree<> stree; + GIC.create_complex(stree); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if (verb) { + std::cout << "Graph induced complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on graph induced complex simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << std::endl; + } + } + } + + return 0; +} diff --git a/utilities/Nerve_GIC/covercomplex.md b/utilities/Nerve_GIC/covercomplex.md new file mode 100644 index 00000000..f33cb2e0 --- /dev/null +++ b/utilities/Nerve_GIC/covercomplex.md @@ -0,0 +1,62 @@ + + +# Cover complex # + + +## Nerve ## +This program builds the Nerve of a point cloud sampled on an OFF file. +The cover C comes from the preimages of intervals covering a coordinate function, +which are then refined into their connected components using the triangulation of the .OFF file. + +The program also writes a file SC.txt. +The first three lines in this file are the location of the input point cloud and the function used to compute the cover. +The fourth line contains the number of vertices nv and edges ne of the Nerve. The next nv lines represent the vertices. +Each line contains the vertex ID, the number of data points it contains, and their average color function value. +Finally, the next ne lines represent the edges, characterized by the ID of their vertices. + +**Usage** + +`Nerve coordinate resolution gain [-v]` + +where + +* `coordinate` is the coordinate function to cover +* `resolution` is the number of the intervals +* `gain` is the gain for each interval +* `-v` is optional, it activates verbose mode. + +**Example** + +`Nerve ../../data/points/human.off 2 10 0.3` + +* Builds the Nerve of a point cloud sampled on a 3D human shape (human.off). +The cover C comes from the preimages of intervals (10 intervals with gain 0.3) covering the height function (coordinate 2). + +`python KeplerMapperVisuFromTxtFile.py -f ../../data/points/human.off_sc.txt` + +* Constructs `human.off_sc.html` file. You can now use your favorite web browser to visualize it. + +## VoronoiGIC ## + +This util builds the Graph Induced Complex (GIC) of a point cloud. +It subsamples *N* points in the point cloud, which act as seeds of a geodesic Voronoï diagram. +Each cell of the diagram is then an element of C. + +The program also writes a file `*_sc.off`, that is an OFF file that can be visualized with GeomView. + +**Usage** + +`VoroniGIC samples_number [-v]` + +where + +* `samples_number` is the number of samples to take from the point cloud +* `-v` is optional, it activates verbose mode. + +**Example** + +`VoroniGIC ../../data/points/human.off 700` + +* Builds the Voronoi Graph Induced Complex with 700 subsamples from `human.off` file. +`../../data/points/human_sc.off` can be visualized with GeomView. + diff --git a/utilities/Nerve_GIC/km.py b/utilities/Nerve_GIC/km.py new file mode 100755 index 00000000..53024aab --- /dev/null +++ b/utilities/Nerve_GIC/km.py @@ -0,0 +1,390 @@ +from __future__ import division +import numpy as np +from collections import defaultdict +import json +import itertools +from sklearn import cluster, preprocessing, manifold +from datetime import datetime +import sys + +class KeplerMapper(object): + # With this class you can build topological networks from (high-dimensional) data. + # + # 1) Fit a projection/lens/function to a dataset and transform it. + # For instance "mean_of_row(x) for x in X" + # 2) Map this projection with overlapping intervals/hypercubes. + # Cluster the points inside the interval + # (Note: we cluster on the inverse image/original data to lessen projection loss). + # If two clusters/nodes have the same members (due to the overlap), then: + # connect these with an edge. + # 3) Visualize the network using HTML and D3.js. + # + # functions + # --------- + # fit_transform: Create a projection (lens) from a dataset + # map: Apply Mapper algorithm on this projection and build a simplicial complex + # visualize: Turns the complex dictionary into a HTML/D3.js visualization + + def __init__(self, verbose=2): + self.verbose = verbose + + self.chunk_dist = [] + self.overlap_dist = [] + self.d = [] + self.nr_cubes = 0 + self.overlap_perc = 0 + self.clusterer = False + + def fit_transform(self, X, projection="sum", scaler=preprocessing.MinMaxScaler()): + # Creates the projection/lens from X. + # + # Input: X. Input features as a numpy array. + # Output: projected_X. original data transformed to a projection (lens). + # + # parameters + # ---------- + # projection: Projection parameter is either a string, + # a scikit class with fit_transform, like manifold.TSNE(), + # or a list of dimension indices. + # scaler: if None, do no scaling, else apply scaling to the projection + # Default: Min-Max scaling + + self.scaler = scaler + self.projection = str(projection) + + # Detect if projection is a class (for scikit-learn) + #if str(type(projection))[1:6] == "class": #TODO: de-ugly-fy + # reducer = projection + # if self.verbose > 0: + # try: + # projection.set_params(**{"verbose":self.verbose}) + # except: + # pass + # print("\n..Projecting data using: \n\t%s\n"%str(projection)) + # X = reducer.fit_transform(X) + + # Detect if projection is a string (for standard functions) + if isinstance(projection, str): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(projection)) + # Stats lenses + if projection == "sum": # sum of row + X = np.sum(X, axis=1).reshape((X.shape[0],1)) + if projection == "mean": # mean of row + X = np.mean(X, axis=1).reshape((X.shape[0],1)) + if projection == "median": # mean of row + X = np.median(X, axis=1).reshape((X.shape[0],1)) + if projection == "max": # max of row + X = np.max(X, axis=1).reshape((X.shape[0],1)) + if projection == "min": # min of row + X = np.min(X, axis=1).reshape((X.shape[0],1)) + if projection == "std": # std of row + X = np.std(X, axis=1).reshape((X.shape[0],1)) + + if projection == "dist_mean": # Distance of x to mean of X + X_mean = np.mean(X, axis=0) + X = np.sum(np.sqrt((X - X_mean)**2), axis=1).reshape((X.shape[0],1)) + + # Detect if projection is a list (with dimension indices) + if isinstance(projection, list): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(str(projection))) + X = X[:,np.array(projection)] + + # Scaling + if scaler is not None: + if self.verbose > 0: + print("\n..Scaling with: %s\n"%str(scaler)) + X = scaler.fit_transform(X) + + return X + + def map(self, projected_X, inverse_X=None, clusterer=cluster.DBSCAN(eps=0.5,min_samples=3), nr_cubes=10, overlap_perc=0.1): + # This maps the data to a simplicial complex. Returns a dictionary with nodes and links. + # + # Input: projected_X. A Numpy array with the projection/lens. + # Output: complex. A dictionary with "nodes", "links" and "meta information" + # + # parameters + # ---------- + # projected_X projected_X. A Numpy array with the projection/lens. Required. + # inverse_X Numpy array or None. If None then the projection itself is used for clustering. + # clusterer Scikit-learn API compatible clustering algorithm. Default: DBSCAN + # nr_cubes Int. The number of intervals/hypercubes to create. + # overlap_perc Float. The percentage of overlap "between" the intervals/hypercubes. + + start = datetime.now() + + # Helper function + def cube_coordinates_all(nr_cubes, nr_dimensions): + # Helper function to get origin coordinates for our intervals/hypercubes + # Useful for looping no matter the number of cubes or dimensions + # Example: if there are 4 cubes per dimension and 3 dimensions + # return the bottom left (origin) coordinates of 64 hypercubes, + # as a sorted list of Numpy arrays + # TODO: elegance-ify... + l = [] + for x in range(nr_cubes): + l += [x] * nr_dimensions + return [np.array(list(f)) for f in sorted(set(itertools.permutations(l,nr_dimensions)))] + + nodes = defaultdict(list) + links = defaultdict(list) + complex = {} + self.nr_cubes = nr_cubes + self.clusterer = clusterer + self.overlap_perc = overlap_perc + + if self.verbose > 0: + print("Mapping on data shaped %s using dimensions\n"%(str(projected_X.shape))) + + # If inverse image is not provided, we use the projection as the inverse image (suffer projection loss) + if inverse_X is None: + inverse_X = projected_X + + # We chop up the min-max column ranges into 'nr_cubes' parts + self.chunk_dist = (np.max(projected_X, axis=0) - np.min(projected_X, axis=0))/nr_cubes + + # We calculate the overlapping windows distance + self.overlap_dist = self.overlap_perc * self.chunk_dist + + # We find our starting point + self.d = np.min(projected_X, axis=0) + + # Use a dimension index array on the projected X + # (For now this uses the entire dimensionality, but we keep for experimentation) + di = np.array([x for x in range(projected_X.shape[1])]) + + # Prefix'ing the data with ID's + ids = np.array([x for x in range(projected_X.shape[0])]) + projected_X = np.c_[ids,projected_X] + inverse_X = np.c_[ids,inverse_X] + + # Subdivide the projected data X in intervals/hypercubes with overlap + if self.verbose > 0: + total_cubes = len(cube_coordinates_all(nr_cubes,projected_X.shape[1])) + print("Creating %s hypercubes."%total_cubes) + + for i, coor in enumerate(cube_coordinates_all(nr_cubes,di.shape[0])): + # Slice the hypercube + hypercube = projected_X[ np.invert(np.any((projected_X[:,di+1] >= self.d[di] + (coor * self.chunk_dist[di])) & + (projected_X[:,di+1] < self.d[di] + (coor * self.chunk_dist[di]) + self.chunk_dist[di] + self.overlap_dist[di]) == False, axis=1 )) ] + + if self.verbose > 1: + print("There are %s points in cube_%s / %s with starting range %s"% + (hypercube.shape[0],i,total_cubes,self.d[di] + (coor * self.chunk_dist[di]))) + + # If at least one sample inside the hypercube + if hypercube.shape[0] > 0: + # Cluster the data point(s) in the cube, skipping the id-column + # Note that we apply clustering on the inverse image (original data samples) that fall inside the cube. + inverse_x = inverse_X[[int(nn) for nn in hypercube[:,0]]] + + clusterer.fit(inverse_x[:,1:]) + + if self.verbose > 1: + print("Found %s clusters in cube_%s\n"%(np.unique(clusterer.labels_[clusterer.labels_ > -1]).shape[0],i)) + + #Now for every (sample id in cube, predicted cluster label) + for a in np.c_[hypercube[:,0],clusterer.labels_]: + if a[1] != -1: #if not predicted as noise + cluster_id = str(coor[0])+"_"+str(i)+"_"+str(a[1])+"_"+str(coor)+"_"+str(self.d[di] + (coor * self.chunk_dist[di])) # TODO: de-rudimentary-ify + nodes[cluster_id].append( int(a[0]) ) # Append the member id's as integers + else: + if self.verbose > 1: + print("Cube_%s is empty.\n"%(i)) + + # Create links when clusters from different hypercubes have members with the same sample id. + candidates = itertools.combinations(nodes.keys(),2) + for candidate in candidates: + # if there are non-unique members in the union + if len(nodes[candidate[0]]+nodes[candidate[1]]) != len(set(nodes[candidate[0]]+nodes[candidate[1]])): + links[candidate[0]].append( candidate[1] ) + + # Reporting + if self.verbose > 0: + nr_links = 0 + for k in links: + nr_links += len(links[k]) + print("\ncreated %s edges and %s nodes in %s."%(nr_links,len(nodes),str(datetime.now()-start))) + + complex["nodes"] = nodes + complex["links"] = links + complex["meta"] = self.projection + + return complex + + def visualize(self, complex, color_function="", path_html="mapper_visualization_output.html", title="My Data", + graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=None, width_html=0, + height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=0,gain=0,minimum=0,maximum=0): + # Turns the dictionary 'complex' in a html file with d3.js + # + # Input: complex. Dictionary (output from calling .map()) + # Output: a HTML page saved as a file in 'path_html'. + # + # parameters + # ---------- + # color_function string. Not fully implemented. Default: "" (distance to origin) + # path_html file path as string. Where to save the HTML page. + # title string. HTML page document title and first heading. + # graph_link_distance int. Edge length. + # graph_gravity float. "Gravity" to center of layout. + # graph_charge int. charge between nodes. + # custom_tooltips None or Numpy Array. You could use "y"-label array for this. + # width_html int. Width of canvas. Default: 0 (full width) + # height_html int. Height of canvas. Default: 0 (full height) + # show_tooltips bool. default:True + # show_title bool. default:True + # show_meta bool. default:True + + # Format JSON for D3 graph + json_s = {} + json_s["nodes"] = [] + json_s["links"] = [] + k2e = {} # a key to incremental int dict, used for id's when linking + + for e, k in enumerate(complex["nodes"]): + # Tooltip and node color formatting, TODO: de-mess-ify + if custom_tooltips is not None: + tooltip_s = "

Cluster %s

"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" ")) + if maximum == minimum: + tooltip_i = 0 + else: + tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum)) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i}) + else: + tooltip_s = "

Cluster %s

Contains %s members."%(k,len(complex["nodes"][k])) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])}) + k2e[k] = e + for k in complex["links"]: + for link in complex["links"][k]: + json_s["links"].append({"source": k2e[k], "target":k2e[link],"value":1}) + + # Width and height of graph in HTML output + if width_html == 0: + width_css = "100%" + width_js = 'document.getElementById("holder").offsetWidth-20' + else: + width_css = "%spx" % width_html + width_js = "%s" % width_html + if height_html == 0: + height_css = "100%" + height_js = 'document.getElementById("holder").offsetHeight-20' + else: + height_css = "%spx" % height_html + height_js = "%s" % height_html + + # Whether to show certain UI elements or not + if show_tooltips == False: + tooltips_display = "display: none;" + else: + tooltips_display = "" + + if show_meta == False: + meta_display = "display: none;" + else: + meta_display = "" + + if show_title == False: + title_display = "display: none;" + else: + title_display = "" + + with open(path_html,"wb") as outfile: + html = """ + + + %s | KeplerMapper + + + +
+

%s

+

+ Lens
%s

+ Length of intervals
%s

+ Overlap percentage
%s%%

+ Color Function
%s +

+
+ + """%(title,width_css, height_css, title_display, meta_display, tooltips_display, title,complex["meta"],res,gain*100,color_function,width_js,height_js,graph_charge,graph_link_distance,graph_gravity,json.dumps(json_s)) + outfile.write(html.encode("utf-8")) + if self.verbose > 0: + print("\nWrote d3.js graph to '%s'"%path_html) diff --git a/utilities/Nerve_GIC/km.py.COPYRIGHT b/utilities/Nerve_GIC/km.py.COPYRIGHT new file mode 100644 index 00000000..bef7b121 --- /dev/null +++ b/utilities/Nerve_GIC/km.py.COPYRIGHT @@ -0,0 +1,26 @@ +km.py is a fork of https://github.com/MLWave/kepler-mapper. +Only the visualization part has been kept (Mapper part has been removed). + +This file has te following Copyright : + +The MIT License (MIT) + +Copyright (c) 2015 Triskelion - HJ van Veen - info@mlwave.com + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/utilities/Persistence_representations/CMakeLists.txt b/utilities/Persistence_representations/CMakeLists.txt new file mode 100644 index 00000000..137eb0c1 --- /dev/null +++ b/utilities/Persistence_representations/CMakeLists.txt @@ -0,0 +1,53 @@ +# Copy files, otherwise result files are created in sources +file(COPY "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") +file(COPY "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") + +function(add_persistence_representation_creation_utility creation_utility) + add_executable ( ${creation_utility} ${creation_utility}.cpp ) + + # as the function is called in a subdirectory level, need to '../' to find persistence files + # ARGN will add all the other arguments (except creation_utility) sent to the CMake functions + add_test(NAME Persistence_representation_utilities_${creation_utility} COMMAND $ + ${ARGN} "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" + "${CMAKE_CURRENT_BINARY_DIR}/../second.pers") +endfunction(add_persistence_representation_creation_utility) + +function(add_persistence_representation_plot_utility plot_utility tool_extension) + add_executable ( ${plot_utility} ${plot_utility}.cpp ) + + # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files + add_test(NAME Persistence_representation_utilities_${plot_utility}_first COMMAND $ + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}") + #add_test(NAME Persistence_representation_utilities_${plot_utility}_second COMMAND $ + # "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}") + if(GNUPLOT_PATH) + add_test(NAME Persistence_representation_utilities_${plot_utility}_first_gnuplot COMMAND ${GNUPLOT_PATH} + "-e" "load '${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}_GnuplotScript'") + #add_test(NAME Persistence_representation_utilities_${plot_utility}_second_gnuplot COMMAND ${GNUPLOT_PATH} + # "-e" "load '${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}_GnuplotScript'") + endif() +endfunction(add_persistence_representation_plot_utility) + +function(add_persistence_representation_function_utility function_utility tool_extension) + add_executable ( ${function_utility} ${function_utility}.cpp ) + + # ARGV2 is an optional argument + if (${ARGV2}) + # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files + add_test(NAME Persistence_representation_utilities_${function_utility} COMMAND $ + "${ARGV2}" + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}" + "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}") + else() + # as the function is called in a subdirectory level, need to '../' to find persistence heat maps files + add_test(NAME Persistence_representation_utilities_${function_utility} COMMAND $ + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers${tool_extension}" + "${CMAKE_CURRENT_BINARY_DIR}/../second.pers${tool_extension}") + endif() +endfunction(add_persistence_representation_function_utility) + +add_subdirectory(persistence_heat_maps) +add_subdirectory(persistence_intervals) +add_subdirectory(persistence_landscapes) +add_subdirectory(persistence_landscapes_on_grid) +add_subdirectory(persistence_vectors) diff --git a/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt b/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt new file mode 100644 index 00000000..386e9fa5 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistence_representations_heat_maps_utilities) + +add_persistence_representation_creation_utility(create_pssk "10" "-1" "-1" "4" "-1") +add_persistence_representation_creation_utility(create_p_h_m_weighted_by_arctan_of_their_persistence "10" "-1" "-1" "4" "-1") +add_persistence_representation_creation_utility(create_p_h_m_weighted_by_distance_from_diagonal "10" "-1" "-1" "4" "-1") +add_persistence_representation_creation_utility(create_p_h_m_weighted_by_squared_diag_distance "10" "-1" "-1" "4" "-1") +# Need to set grid min and max for further average, distance and scalar_product +add_persistence_representation_creation_utility(create_persistence_heat_maps "10" "0" "35" "10" "-1") + +add_persistence_representation_plot_utility(plot_persistence_heat_map ".mps") + +add_persistence_representation_function_utility(average_persistence_heat_maps ".mps") +add_persistence_representation_function_utility(compute_distance_of_persistence_heat_maps ".mps" "1") +add_persistence_representation_function_utility(compute_scalar_product_of_persistence_heat_maps ".mps") diff --git a/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp new file mode 100644 index 00000000..6739e0b6 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/average_persistence_heat_maps.cpp @@ -0,0 +1,63 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program computes average of persistence heat maps stored in files (the files needs to be " + << "created beforehand).\n" + << "The parameters of this programs are names of files with persistence heat maps.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector maps; + for (size_t i = 0; i != filenames.size(); ++i) { + Persistence_heat_maps* l = new Persistence_heat_maps; + l->load_from_file(filenames[i]); + maps.push_back(l); + } + + Persistence_heat_maps av; + av.compute_average(maps); + av.print_to_file("average.mps"); + + for (size_t i = 0; i != filenames.size(); ++i) { + delete maps[i]; + } + + std::cout << "Average can be found in 'average.mps' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp new file mode 100644 index 00000000..ed8278a2 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/compute_distance_of_persistence_heat_maps.cpp @@ -0,0 +1,94 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program computes distance of persistence heat maps stored in files (the files needs to be " + << "created beforehand).\n" + << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat " + << "maps. For L^infty distance choose p = -1. \n" + << "The remaining parameters of this program are names of files with persistence heat maps.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + int pp = atoi(argv[1]); + double p = std::numeric_limits::max(); + if (pp != -1) { + p = pp; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector maps; + maps.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_heat_maps l; + l.load_from_file(filenames[file_no]); + maps.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > distance(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + distance[i] = v; + } + + // and now we can compute the distances: + for (size_t i = 0; i != filenames.size(); ++i) { + for (size_t j = i; j != filenames.size(); ++j) { + distance[i][j] = distance[j][i] = maps[i].distance(maps[j], p); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("distance.mps"); + for (size_t i = 0; i != distance.size(); ++i) { + for (size_t j = 0; j != distance.size(); ++j) { + std::cout << distance[i][j] << " "; + out << distance[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'distance.mps' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp new file mode 100644 index 00000000..63626853 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/compute_scalar_product_of_persistence_heat_maps.cpp @@ -0,0 +1,85 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program computes scalar product of persistence heat maps stored in a file (the file needs to be " + << "created beforehand). \n" + << "The parameters of this programs are names of files with persistence heat maps.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector maps; + maps.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_heat_maps l; + l.load_from_file(filenames[file_no]); + maps.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > scalar_product(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + scalar_product[i] = v; + } + + // and now we can compute the scalar product: + for (size_t i = 0; i != maps.size(); ++i) { + for (size_t j = i; j != maps.size(); ++j) { + scalar_product[i][j] = scalar_product[j][i] = maps[i].compute_scalar_product(maps[j]); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("scalar_product.mps"); + for (size_t i = 0; i != scalar_product.size(); ++i) { + for (size_t j = 0; j != scalar_product.size(); ++j) { + std::cout << scalar_product[i][j] << " "; + out << scalar_product[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'scalar_product.mps' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp new file mode 100644 index 00000000..b4a1daa5 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_arctan_of_their_persistence.cpp @@ -0,0 +1,81 @@ +/* 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 +* +* 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 . +*/ + +#include + +#include +#include +#include +#include + +using arc_tan_of_persistence_of_point = Gudhi::Persistence_representations::arc_tan_of_persistence_of_point; +using Persistence_heat_maps = + Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) " + << "provided as an input.The Gaussian kernels are weighted by the arc tangential of their persistence.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number " + << "of pixels.\n" + << "The fifth parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 7) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + size_t stdiv = atof(argv[4]); + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[5]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 6; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1); + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl; + Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension); + + std::stringstream ss; + ss << filenames[i] << ".mps"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp new file mode 100644 index 00000000..c50f9ddb --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_distance_from_diagonal.cpp @@ -0,0 +1,81 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using distance_from_diagonal_scaling = Gudhi::Persistence_representations::distance_from_diagonal_scaling; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) " + << "provided as an input.The Gaussian kernels are weighted by the distance of a center from the " + << "diagonal.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number " + << "of pixels.\n" + << "The fifth parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 7) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + size_t stdiv = atof(argv[4]); + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[5]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 6; i != argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1); + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl; + Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension); + + std::stringstream ss; + ss << filenames[i] << ".mps"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp new file mode 100644 index 00000000..59ff3c24 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/create_p_h_m_weighted_by_squared_diag_distance.cpp @@ -0,0 +1,83 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using squared_distance_from_diagonal_scaling = + Gudhi::Persistence_representations::squared_distance_from_diagonal_scaling; +using Persistence_heat_maps = + Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) " + << "provided as an input.The Gaussian kernels are weighted by the square of distance of a center from the " + << "diagonal.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number " + << "of pixels.\n" + << "The fifth parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 7) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + size_t stdiv = atof(argv[4]); + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[5]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 6; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1); + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a heat map based on a file : " << filenames[i] << std::endl; + Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension); + + std::stringstream ss; + ss << filenames[i] << ".mps"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp new file mode 100644 index 00000000..25cd1067 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/create_persistence_heat_maps.cpp @@ -0,0 +1,78 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence heat map files (*.mps) of persistence diagrams files (*.pers) " + << "provided as an input.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number " + << "of pixels.\n" + << "The fifth parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the fifth parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 7) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + size_t stdiv = atof(argv[4]); + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[5]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + std::vector filenames; + for (int i = 6; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1); + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a heat map based on file : " << filenames[i] << std::endl; + Persistence_heat_maps l(filenames[i], filter, false, size_of_grid, min_, max_, dimension); + + std::stringstream ss; + ss << filenames[i] << ".mps"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp b/utilities/Persistence_representations/persistence_heat_maps/create_pssk.cpp new file mode 100644 index 00000000..97ddb8f0 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/create_pssk.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): Pawel Dlotko + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using PSSK = Gudhi::Persistence_representations::PSSK; + +int main(int argc, char** argv) { + std::cout << "This program creates PSSK files (*.pssk) of persistence diagrams files (*.pers) " + << "provided as an input.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter is an integer, the standard deviation of a Gaussian kernel expressed in a number " + << "of pixels.\n" + << "The fifth parameter of this program is a dimension of persistence that will be used in creation of " + << "the PSSK." + << "If your input files contains persistence pairs of various dimension, as a fifth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 7) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + size_t stdiv = atof(argv[4]); + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[5]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 6; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector > filter = Gudhi::Persistence_representations::create_Gaussian_filter(stdiv, 1); + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a PSSK based on a file : " << filenames[i] << std::endl; + PSSK l(filenames[i], filter, size_of_grid, min_, max_, dimension); + + std::stringstream ss; + ss << filenames[i] << ".pssk"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp b/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp new file mode 100644 index 00000000..63711d83 --- /dev/null +++ b/utilities/Persistence_representations/persistence_heat_maps/plot_persistence_heat_map.cpp @@ -0,0 +1,42 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps; + +int main(int argc, char** argv) { + std::cout << "This program creates a gnuplot script from a persistence heat maps stored in a file (the file needs " + << "to be created beforehand). Please call the code with the name of a single heat maps file \n"; + if (argc != 2) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + Persistence_heat_maps l; + l.load_from_file(argv[1]); + l.plot(argv[1]); + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt b/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt new file mode 100644 index 00000000..897e12a3 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/CMakeLists.txt @@ -0,0 +1,32 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistence_representations_intervals_utilities) + + +add_executable ( plot_histogram_of_intervals_lengths plot_histogram_of_intervals_lengths.cpp ) + +add_test(NAME plot_histogram_of_intervals_lengths COMMAND $ + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" "-1") + +add_persistence_representation_plot_utility(plot_persistence_intervals "") +add_persistence_representation_plot_utility(plot_persistence_Betti_numbers "") + +add_persistence_representation_creation_utility(compute_birth_death_range_in_persistence_diagram "-1") + + +add_executable ( compute_number_of_dominant_intervals compute_number_of_dominant_intervals.cpp ) +add_test(NAME Persistence_representation_utilities_compute_number_of_dominant_intervals + COMMAND $ + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" "-1" "2") + + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) + add_executable ( compute_bottleneck_distance compute_bottleneck_distance.cpp ) + if (TBB_FOUND) + target_link_libraries(compute_bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Persistence_representation_utilities_compute_bottleneck_distance + COMMAND $ + "-1" + "${CMAKE_CURRENT_BINARY_DIR}/../first.pers" + "${CMAKE_CURRENT_BINARY_DIR}/../second.pers") +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp b/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp new file mode 100644 index 00000000..9102da79 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/compute_birth_death_range_in_persistence_diagram.cpp @@ -0,0 +1,68 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals; + +int main(int argc, char** argv) { + std::cout << "This program computes the range of birth and death times of persistence pairs in diagrams provided as " + << "an input.\n" + << "The first parameter is the dimension of persistence to be used to create persistence intervals. " + << "If your file contains the information about dimension of persistence pairs, please provide here the " + << "dimension of persistence pairs you want to use. " + << "If your input files consist only of birth-death pairs, please set this first parameter to -1.\n" + << "The remaining parameters of the program are the names of files with persistence diagrams.\n"; + + if (argc < 3) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[1]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + double min_ = std::numeric_limits::max(); + double max_ = -std::numeric_limits::max(); + + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + std::cout << "Creating diagram based on a file : " << filenames[file_no] << std::endl; + Persistence_intervals p(filenames[file_no], dimension); + std::pair min_max_ = p.get_x_range(); + if (min_max_.first < min_) min_ = min_max_.first; + if (min_max_.second > max_) max_ = min_max_.second; + } + std::cout << "Birth-death range : min: " << min_ << ", max: " << max_ << std::endl; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp b/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp new file mode 100644 index 00000000..c8290845 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/compute_bottleneck_distance.cpp @@ -0,0 +1,95 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals_with_distances = Gudhi::Persistence_representations::Persistence_intervals_with_distances; + +int main(int argc, char** argv) { + std::cout << "This program computes the bottleneck distance of persistence pairs in diagrams provided as " + << "an input.\n" + << "The first parameter is the dimension of persistence to be used to create persistence intervals. " + << "If your file contains the information about dimension of persistence pairs, please provide here the " + << "dimension of persistence pairs you want to use. " + << "If your input files consist only of birth-death pairs, please set this first parameter to -1.\n" + << "The remaining parameters of the program are the names of files with persistence diagrams.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + unsigned dimension = std::numeric_limits::max(); + int dim = atoi(argv[1]); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + // reading the persistence intervals: + std::vector persistence_intervals; + for (size_t i = 0; i != filenames.size(); ++i) { + Persistence_intervals_with_distances pers(filenames[i], dimension); + persistence_intervals.push_back(pers); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > distance(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + distance[i] = v; + } + + // and now we can compute the distances: + for (size_t i = 0; i != persistence_intervals.size(); ++i) { + for (size_t j = i + 1; j != persistence_intervals.size(); ++j) { + distance[i][j] = distance[j][i] = persistence_intervals[i].distance(persistence_intervals[j]); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("distance.itv"); + for (size_t i = 0; i != distance.size(); ++i) { + for (size_t j = 0; j != distance.size(); ++j) { + std::cout << distance[i][j] << " "; + out << distance[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'distance.itv' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp b/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp new file mode 100644 index 00000000..b3d126f0 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/compute_number_of_dominant_intervals.cpp @@ -0,0 +1,54 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals; + +int main(int argc, char** argv) { + std::cout << "This program compute the dominant intervals. A number of intervals to be displayed is a parameter of " + "this program. \n"; + if (argc != 4) { + std::cout << "To run this program, please provide the name of a file with persistence diagram, dimension of " + "intervals that should be taken into account (if your file contains only persistence pairs in a " + "single dimension, set it up to -1) and number of dominant intervals you would like to get \n"; + return 1; + } + int dim = atoi(argv[2]); + unsigned dimension = std::numeric_limits::max(); + if (dim >= 0) { + dimension = (unsigned)dim; + } + Persistence_intervals p(argv[1], dimension); + std::vector > dominant_intervals = p.dominant_intervals(atoi(argv[3])); + std::cout << "Here are the dominant intervals : " << std::endl; + for (size_t i = 0; i != dominant_intervals.size(); ++i) { + std::cout << " " << dominant_intervals[i].first << "," << dominant_intervals[i].second << " " << std::endl; + } + + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp b/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp new file mode 100644 index 00000000..ccb5b645 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/plot_histogram_of_intervals_lengths.cpp @@ -0,0 +1,77 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals; + +int main(int argc, char** argv) { + std::cout << "This program computes a histogram of barcode's length. A number of bins in the histogram is a " + << "parameter of this program. \n"; + if ((argc != 3) && (argc != 4)) { + std::cout << "To run this program, please provide the name of a file with persistence diagram and number of " + << "dominant intervals you would like to get. Set a negative number dominant intervals value " + << "If your file contains only birth-death pairs.\n" + << "The third parameter is the dimension of the persistence that is to be used. If your " + << "file contains only birth-death pairs, you can skip this parameter\n"; + return 1; + } + + unsigned dominant_interval_number = std::numeric_limits::max(); + int nbr = atoi(argv[2]); + if (nbr >= 0) { + dominant_interval_number = static_cast(nbr); + } + + int persistence_dimension = -1; + if (argc == 4) { + persistence_dimension = atoi(argv[3]); + } + + Persistence_intervals p(argv[1], dominant_interval_number); + std::vector > dominant_intervals = p.dominant_intervals(persistence_dimension); + std::vector histogram = p.histogram_of_lengths(10); + + std::stringstream gnuplot_script; + gnuplot_script << argv[1] << "_GnuplotScript"; + std::ofstream out; + out.open(gnuplot_script.str().c_str()); + + out << "set style data histogram" << std::endl; + out << "set style histogram cluster gap 1" << std::endl; + out << "set style fill solid border -1" << std::endl; + out << "plot '-' notitle" << std::endl; + for (size_t i = 0; i != histogram.size(); ++i) { + out << histogram[i] << std::endl; + } + out << std::endl; + out.close(); + + std::cout << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << std::endl; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp b/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp new file mode 100644 index 00000000..b433c2b3 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/plot_persistence_Betti_numbers.cpp @@ -0,0 +1,87 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals; + +int main(int argc, char** argv) { + if ((argc != 3) && (argc != 2)) { + std::cout << "This program creates a gnuplot script of Betti numbers from a single persistence diagram file" + << "(*.pers).\n" + << "To run this program, please provide the name of a file with persistence diagram.\n" + << "The second optional parameter of a program is the dimension of the persistence that is to be used. " + << "If your file contains only birth-death pairs, you can skip this parameter.\n"; + return 1; + } + + unsigned dimension = std::numeric_limits::max(); + int dim = -1; + if (argc == 3) { + dim = atoi(argv[2]); + } + if (dim >= 0) { + dimension = (unsigned)dim; + } + + Persistence_intervals p(argv[1], dimension); + std::vector > pbns = p.compute_persistent_betti_numbers(); + + // set up the ranges so that we see the image well. + double xRangeBegin = pbns[0].first; + double xRangeEnd = pbns[pbns.size() - 1].first; + double yRangeBegin = 0; + double yRangeEnd = 0; + for (size_t i = 0; i != pbns.size(); ++i) { + if (pbns[i].second > yRangeEnd) yRangeEnd = pbns[i].second; + } + xRangeBegin -= (xRangeEnd - xRangeBegin) / 100.0; + xRangeEnd += (xRangeEnd - xRangeBegin) / 100.0; + yRangeEnd += yRangeEnd / 100; + + std::stringstream gnuplot_script; + gnuplot_script << argv[1] << "_GnuplotScript"; + std::ofstream out; + out.open(gnuplot_script.str().c_str()); + + out << "set xrange [" << xRangeBegin << " : " << xRangeEnd << "]" << std::endl; + out << "set yrange [" << yRangeBegin << " : " << yRangeEnd << "]" << std::endl; + out << "plot '-' using 1:2 notitle with lp " << std::endl; + double previous_y = 0; + for (size_t i = 0; i != pbns.size(); ++i) { + out << pbns[i].first << " " << previous_y << std::endl; + out << pbns[i].first << " " << pbns[i].second << std::endl; + previous_y = pbns[i].second; + } + out << std::endl; + out.close(); + + std::cout << "To visualize, install gnuplot and type the command: gnuplot -persist -e \"load \'" + << gnuplot_script.str().c_str() << "\'\"" << std::endl; + + return 0; +} diff --git a/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp b/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp new file mode 100644 index 00000000..33387802 --- /dev/null +++ b/utilities/Persistence_representations/persistence_intervals/plot_persistence_intervals.cpp @@ -0,0 +1,53 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals; + +int main(int argc, char** argv) { + if ((argc != 3) && (argc != 2)) { + std::cout << "This program creates a gnuplot script from a single persistence diagram file (*.pers).\n" + << "To run this program, please provide the name of a file with persistence diagram.\n" + << "The second optional parameter of a program is the dimension of the persistence that is to be used. " + << "If your file contains only birth-death pairs, you can skip this parameter.\n"; + return 1; + } + unsigned dimension = std::numeric_limits::max(); + int dim = -1; + if (argc == 3) { + dim = atoi(argv[2]); + } + if (dim >= 0) { + dimension = (unsigned)dim; + } + std::vector > intervals = + Gudhi::Persistence_representations::read_persistence_intervals_in_one_dimension_from_file(argv[1], dimension); + Persistence_intervals b(intervals); + b.plot(argv[1]); + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt b/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt new file mode 100644 index 00000000..d7087ed8 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistence_representations_landscapes_utilities) + +add_persistence_representation_creation_utility(create_landscapes "-1") + +add_persistence_representation_plot_utility(plot_landscapes ".land") + +add_persistence_representation_function_utility(average_landscapes ".land") +add_persistence_representation_function_utility(compute_distance_of_landscapes ".land" "1") +add_persistence_representation_function_utility(compute_scalar_product_of_landscapes ".land") diff --git a/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp new file mode 100644 index 00000000..1a59be8c --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/average_landscapes.cpp @@ -0,0 +1,63 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape; + +int main(int argc, char** argv) { + std::cout << "This program computes average of persistence landscapes stored in files (the files needs to be " + << "created beforehand).\n" + << "The parameters of this programs are names of files with persistence landscapes.\n"; + std::vector filenames; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector lands; + for (size_t i = 0; i != filenames.size(); ++i) { + Persistence_landscape* l = new Persistence_landscape; + l->load_landscape_from_file(filenames[i]); + lands.push_back(l); + } + + Persistence_landscape av; + av.compute_average(lands); + + av.print_to_file("average.land"); + + for (size_t i = 0; i != filenames.size(); ++i) { + delete lands[i]; + } + + std::cout << "Average can be found in 'average.land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp new file mode 100644 index 00000000..5062f521 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/compute_distance_of_landscapes.cpp @@ -0,0 +1,93 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape; + +int main(int argc, char** argv) { + std::cout << "This program computes distance of persistence landscapes stored in files (the files needs to be " + << "created beforehand).\n" + << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat " + << "maps. For L^infty distance choose p = -1. \n" + << "The remaining parameters of this program are names of files with persistence landscapes.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + int pp = atoi(argv[1]); + double p = std::numeric_limits::max(); + if (pp != -1) { + p = pp; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector landscaspes; + landscaspes.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_landscape l; + l.load_landscape_from_file(filenames[file_no]); + landscaspes.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > distance(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + distance[i] = v; + } + + // and now we can compute the distances: + for (size_t i = 0; i != landscaspes.size(); ++i) { + for (size_t j = i; j != landscaspes.size(); ++j) { + distance[i][j] = distance[j][i] = compute_distance_of_landscapes(landscaspes[i], landscaspes[j], p); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("distance.land"); + for (size_t i = 0; i != distance.size(); ++i) { + for (size_t j = 0; j != distance.size(); ++j) { + std::cout << distance[i][j] << " "; + out << distance[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'distance.land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp new file mode 100644 index 00000000..5b5e9fa3 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/compute_scalar_product_of_landscapes.cpp @@ -0,0 +1,84 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include + +using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape; + +int main(int argc, char** argv) { + std::cout << "This program computes scalar product of persistence landscapes stored in a file (the file needs to be " + << "created beforehand). \n" + << "The parameters of this programs are names of files with persistence landscapes.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector landscaspes; + landscaspes.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_landscape l; + l.load_landscape_from_file(filenames[file_no]); + landscaspes.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > scalar_product(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + scalar_product[i] = v; + } + + // and now we can compute the scalar product: + for (size_t i = 0; i != landscaspes.size(); ++i) { + for (size_t j = i; j != landscaspes.size(); ++j) { + scalar_product[i][j] = scalar_product[j][i] = compute_inner_product(landscaspes[i], landscaspes[j]); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("scalar_product.land"); + for (size_t i = 0; i != scalar_product.size(); ++i) { + for (size_t j = 0; j != scalar_product.size(); ++j) { + std::cout << scalar_product[i][j] << " "; + out << scalar_product[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'scalar_product.land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp new file mode 100644 index 00000000..6030e994 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/create_landscapes.cpp @@ -0,0 +1,65 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence landscapes files (*.land) of persistence diagrams files (*.pers) " + << "provided as an input.\n" + << "The first parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a first parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 3) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + std::vector filenames; + int dim = atoi(argv[1]); + unsigned dimension = std::numeric_limits::max(); + if (dim >= 0) { + dimension = (unsigned)dim; + } + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating a landscape based on file : " << filenames[i] << std::endl; + Persistence_landscape l(filenames[i], dimension); + std::stringstream ss; + ss << filenames[i] << ".land"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp b/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp new file mode 100644 index 00000000..c797a7a8 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes/plot_landscapes.cpp @@ -0,0 +1,43 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape; + +int main(int argc, char** argv) { + std::cout << "This program creates a gnuplot script from a persistence landscape stored in a file (the file needs " + << "to be created beforehand). Please call the code with the name of a single landscape file.\n"; + if (argc != 2) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + Persistence_landscape l; + l.load_landscape_from_file(argv[1]); + l.plot(argv[1]); + + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt b/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt new file mode 100644 index 00000000..c5ea4bbf --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/CMakeLists.txt @@ -0,0 +1,11 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistence_representations_lanscapes_on_grid_utilities) + +# Need to set grid min and max for further average, distance and scalar_product +add_persistence_representation_creation_utility(create_landscapes_on_grid "100" "0" "35" "-1") + +add_persistence_representation_plot_utility(plot_landscapes_on_grid ".g_land") + +add_persistence_representation_function_utility(average_landscapes_on_grid ".g_land") +add_persistence_representation_function_utility(compute_distance_of_landscapes_on_grid ".g_land" "1") +add_persistence_representation_function_utility(compute_scalar_product_of_landscapes_on_grid ".g_land") diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp new file mode 100644 index 00000000..0b098d1a --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/average_landscapes_on_grid.cpp @@ -0,0 +1,63 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid; + +int main(int argc, char** argv) { + std::cout << "This program computes average of persistence landscapes on grid stored in files (the files needs to " + << "be created beforehand).\n" + << "The parameters of this programs are names of files with persistence landscapes on grid.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector lands; + for (size_t i = 0; i != filenames.size(); ++i) { + Persistence_landscape_on_grid* l = new Persistence_landscape_on_grid; + l->load_landscape_from_file(filenames[i]); + lands.push_back(l); + } + + Persistence_landscape_on_grid av; + av.compute_average(lands); + + av.print_to_file("average.g_land"); + + for (size_t i = 0; i != filenames.size(); ++i) { + delete lands[i]; + } + + std::cout << "Average can be found in 'average.g_land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp new file mode 100644 index 00000000..fd0fcd15 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_distance_of_landscapes_on_grid.cpp @@ -0,0 +1,93 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid; + +int main(int argc, char** argv) { + std::cout << "This program computes distance of persistence landscapes on grid stored in files (the files needs to " + << "be created beforehand).\n" + << "The first parameter of a program is an integer p. The program compute L^p distance of the two heat " + << "maps. For L^infty distance choose p = -1. \n" + << "The remaining parameters of this program are names of files with persistence landscapes on grid.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + int pp = atoi(argv[1]); + double p = std::numeric_limits::max(); + if (pp != -1) { + p = pp; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector landscaspes; + landscaspes.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_landscape_on_grid l; + l.load_landscape_from_file(filenames[file_no]); + landscaspes.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > distance(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + distance[i] = v; + } + + // and now we can compute the scalar product: + for (size_t i = 0; i != landscaspes.size(); ++i) { + for (size_t j = i; j != landscaspes.size(); ++j) { + distance[i][j] = distance[j][i] = compute_distance_of_landscapes_on_grid(landscaspes[i], landscaspes[j], p); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("distance.g_land"); + for (size_t i = 0; i != distance.size(); ++i) { + for (size_t j = 0; j != distance.size(); ++j) { + std::cout << distance[i][j] << " "; + out << distance[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'distance.g_land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp new file mode 100644 index 00000000..01de3dee --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/compute_scalar_product_of_landscapes_on_grid.cpp @@ -0,0 +1,85 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include + +using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid; + +int main(int argc, char** argv) { + std::cout + << "This program computes scalar product of persistence landscapes on grid stored in a file (the file needs to " + << "be created beforehand). \n" + << "The parameters of this programs are names of files with persistence landscapes on grid.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector landscaspes; + landscaspes.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Persistence_landscape_on_grid l; + l.load_landscape_from_file(filenames[file_no]); + landscaspes.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > scalar_product(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + scalar_product[i] = v; + } + + // and now we can compute the scalar product: + for (size_t i = 0; i != landscaspes.size(); ++i) { + for (size_t j = i; j != landscaspes.size(); ++j) { + scalar_product[i][j] = scalar_product[j][i] = compute_inner_product(landscaspes[i], landscaspes[j]); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("scalar_product.g_land"); + for (size_t i = 0; i != scalar_product.size(); ++i) { + for (size_t j = 0; j != scalar_product.size(); ++j) { + std::cout << scalar_product[i][j] << " "; + out << scalar_product[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'scalar_product.g_land' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.cpp new file mode 100644 index 00000000..78e8ef57 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/create_landscapes_on_grid.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): Pawel Dlotko + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence landscapes on grid files (*.g_land) of persistence diagrams files " + << "(*.pers) provided as an input.\n" + << "The first parameter of a program is an integer, a size of a grid.\n" + << "The second and third parameters are min and max of the grid. If you want those numbers to be computed " + << "based on the data, set them both to -1 \n" + << "The fourth parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a fourth parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the fourth parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 6) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + size_t size_of_grid = (size_t)atoi(argv[1]); + double min_ = atof(argv[2]); + double max_ = atof(argv[3]); + int dim = atoi(argv[4]); + unsigned dimension = std::numeric_limits::max(); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 5; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + for (size_t i = 0; i != filenames.size(); ++i) { + std::cout << "Creating persistence landscape on a grid based on a file : " << filenames[i] << std::endl; + Persistence_landscape_on_grid l; + if ((min_ != -1) || (max_ != -1)) { + l = Persistence_landscape_on_grid(filenames[i], min_, max_, size_of_grid, dimension); + } else { + // (min_ == -1) && (max_ == -1), in this case the program will find min_ and max_ based on the data. + l = Persistence_landscape_on_grid(filenames[i], size_of_grid, dimension); + } + std::stringstream ss; + ss << filenames[i] << ".g_land"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp b/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp new file mode 100644 index 00000000..dddb3615 --- /dev/null +++ b/utilities/Persistence_representations/persistence_landscapes_on_grid/plot_landscapes_on_grid.cpp @@ -0,0 +1,44 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid; + +int main(int argc, char** argv) { + std::cout << "This program creates a gnuplot script from a persistence landscape on grid stored in a file (the file " + << "needs to be created beforehand). Please call the code with the name of a single landscape on grid file" + << ".\n"; + if (argc != 2) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + Persistence_landscape_on_grid l; + l.load_landscape_from_file(argv[1]); + l.plot(argv[1]); + + return 0; +} diff --git a/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt b/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt new file mode 100644 index 00000000..a401c955 --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required(VERSION 2.6) +project(Persistence_vectors_utilities) + +add_persistence_representation_creation_utility(create_persistence_vectors "-1") + +add_persistence_representation_plot_utility(plot_persistence_vectors ".vect") + +add_persistence_representation_function_utility(average_persistence_vectors ".vect") +add_persistence_representation_function_utility(compute_distance_of_persistence_vectors ".vect" "1") +add_persistence_representation_function_utility(compute_scalar_product_of_persistence_vectors ".vect") diff --git a/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp new file mode 100644 index 00000000..0144e76f --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/average_persistence_vectors.cpp @@ -0,0 +1,65 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Euclidean_distance = Gudhi::Euclidean_distance; +using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram; + +int main(int argc, char** argv) { + std::cout << "This program computes average of persistence vectors stored in files (the files needs to " + << "be created beforehand).\n" + << "The parameters of this programs are names of files with persistence vectors.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + std::vector lands; + for (size_t i = 0; i != filenames.size(); ++i) { + Vector_distances_in_diagram* l = new Vector_distances_in_diagram; + l->load_from_file(filenames[i]); + lands.push_back(l); + } + + Vector_distances_in_diagram av; + av.compute_average(lands); + + av.print_to_file("average.vect"); + + for (size_t i = 0; i != filenames.size(); ++i) { + delete lands[i]; + } + + std::cout << "Done \n"; + + return 0; +} diff --git a/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp new file mode 100644 index 00000000..7e66d25e --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/compute_distance_of_persistence_vectors.cpp @@ -0,0 +1,94 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Euclidean_distance = Gudhi::Euclidean_distance; +using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram; + +int main(int argc, char** argv) { + std::cout << "This program compute distance of persistence vectors stored in a file (the file needs to be created " + "beforehand). \n"; + std::cout << "The first parameter of a program is an integer p. The program compute l^p distance of the vectors. For " + "l^infty distance choose p = -1. \n"; + std::cout << "The remaining parameters of this programs are names of files with persistence vectors.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + int pp = atoi(argv[1]); + double p = std::numeric_limits::max(); + if (pp != -1) { + p = pp; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector vectors; + vectors.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Vector_distances_in_diagram l; + l.load_from_file(filenames[file_no]); + vectors.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > distance(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + distance[i] = v; + } + + // and now we can compute the distances: + for (size_t i = 0; i != vectors.size(); ++i) { + for (size_t j = i + 1; j != vectors.size(); ++j) { + distance[i][j] = distance[j][i] = vectors[i].distance(vectors[j], p); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("distance.vect"); + for (size_t i = 0; i != distance.size(); ++i) { + for (size_t j = 0; j != distance.size(); ++j) { + std::cout << distance[i][j] << " "; + out << distance[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'distance.vect' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp new file mode 100644 index 00000000..303c6e3e --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/compute_scalar_product_of_persistence_vectors.cpp @@ -0,0 +1,86 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Euclidean_distance = Gudhi::Euclidean_distance; +using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram; + +int main(int argc, char** argv) { + std::cout << "This program computes scalar product of persistence vectors stored in a file (the file needs to " + << "be created beforehand). \n" + << "The parameters of this programs are names of files with persistence vectors.\n"; + + if (argc < 3) { + std::cout << "Wrong number of parameters, the program will now terminate \n"; + return 1; + } + + std::vector filenames; + for (int i = 1; i < argc; ++i) { + filenames.push_back(argv[i]); + } + std::vector vectors; + vectors.reserve(filenames.size()); + for (size_t file_no = 0; file_no != filenames.size(); ++file_no) { + Vector_distances_in_diagram l; + l.load_from_file(filenames[file_no]); + vectors.push_back(l); + } + + // and now we will compute the scalar product of landscapes. + + // first we prepare an array: + std::vector > scalar_product(filenames.size()); + for (size_t i = 0; i != filenames.size(); ++i) { + std::vector v(filenames.size(), 0); + scalar_product[i] = v; + } + + // and now we can compute the scalar product: + for (size_t i = 0; i != vectors.size(); ++i) { + for (size_t j = i; j != vectors.size(); ++j) { + scalar_product[i][j] = scalar_product[j][i] = vectors[i].compute_scalar_product(vectors[j]); + } + } + + // and now output the result to the screen and a file: + std::ofstream out; + out.open("scalar_product.vect"); + for (size_t i = 0; i != scalar_product.size(); ++i) { + for (size_t j = 0; j != scalar_product.size(); ++j) { + std::cout << scalar_product[i][j] << " "; + out << scalar_product[i][j] << " "; + } + std::cout << std::endl; + out << std::endl; + } + out.close(); + + std::cout << "Distance can be found in 'scalar_product.vect' file\n"; + return 0; +} diff --git a/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.cpp new file mode 100644 index 00000000..cc5e5393 --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/create_persistence_vectors.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): Pawel Dlotko + * + * 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 . + */ + +#include + +#include +#include +#include +#include + +using Euclidean_distance = Gudhi::Euclidean_distance; +using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram; + +int main(int argc, char** argv) { + std::cout << "This program creates persistence vectors files (*.vect) of persistence diagrams files (*.pers) " + << "provided as an input.\n" + << "The first parameter of this program is a dimension of persistence that will be used in creation of " + << "the persistence heat maps." + << "If your input files contains persistence pairs of various dimension, as a first parameter of the " + << "procedure please provide the dimension of persistence you want to use." + << "If in your files there are only birth-death pairs of the same dimension, set the first parameter to " + << "-1.\n" + << "The remaining parameters are the names of files with persistence diagrams. \n"; + + if (argc < 3) { + std::cout << "Wrong parameter list, the program will now terminate \n"; + return 1; + } + + std::cout << "The remaining parameters are the names of files with persistence diagrams. \n"; + int dim = atoi(argv[1]); + unsigned dimension = std::numeric_limits::max(); + if (dim >= 0) { + dimension = (unsigned)dim; + } + + std::vector filenames; + for (int i = 2; i < argc; ++i) { + filenames.push_back(argv[i]); + } + + for (size_t i = 0; i != filenames.size(); ++i) { + std::cerr << "Creating persistence vectors based on a file : " << filenames[i] << std::endl; + Vector_distances_in_diagram l(filenames[i], dimension); + std::stringstream ss; + ss << filenames[i] << ".vect"; + l.print_to_file(ss.str().c_str()); + } + return 0; +} diff --git a/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp b/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp new file mode 100644 index 00000000..aa33107d --- /dev/null +++ b/utilities/Persistence_representations/persistence_vectors/plot_persistence_vectors.cpp @@ -0,0 +1,43 @@ +/* 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 + * + * 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 . + */ + +#include + +#include +#include + +using Euclidean_distance = Gudhi::Euclidean_distance; +using Vector_distances_in_diagram = Gudhi::Persistence_representations::Vector_distances_in_diagram; + +int main(int argc, char** argv) { + std::cout << "This program create a Gnuplot script to plot persistence vector. Please call this program with the " + "name of file with persistence vector. \n"; + if (argc != 2) { + std::cout << "Wrong number of parameters, the program will now terminate. \n"; + return 1; + } + Vector_distances_in_diagram l; + l.load_from_file(argv[1]); + l.plot(argv[1]); + + return 0; +} diff --git a/utilities/Rips_complex/CMakeLists.txt b/utilities/Rips_complex/CMakeLists.txt new file mode 100644 index 00000000..baa571fa --- /dev/null +++ b/utilities/Rips_complex/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required(VERSION 2.6) +project(Rips_complex_utilities) + +add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp) +target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(rips_persistence rips_persistence.cpp) +target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $ + "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") +add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") + +install(TARGETS rips_distance_matrix_persistence DESTINATION bin) +install(TARGETS rips_persistence DESTINATION bin) diff --git a/utilities/Rips_complex/rips_distance_matrix_persistence.cpp b/utilities/Rips_complex/rips_distance_matrix_persistence.cpp new file mode 100644 index 00000000..ca3c0327 --- /dev/null +++ b/utilities/Rips_complex/rips_distance_matrix_persistence.cpp @@ -0,0 +1,133 @@ +/* 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 . + */ + +#include +#include +#include +#include + +#include + +#include +#include +#include // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using Distance_matrix = std::vector>; + +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 = Gudhi::read_lower_triangular_matrix_from_csv_file(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(&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(&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(&threshold)->default_value(std::numeric_limits::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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/utilities/Rips_complex/rips_persistence.cpp b/utilities/Rips_complex/rips_persistence.cpp new file mode 100644 index 00000000..8405c014 --- /dev/null +++ b/utilities/Rips_complex/rips_persistence.cpp @@ -0,0 +1,135 @@ +/* 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 + * + * 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 . + */ + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using Point = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; + +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree 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& 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(&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(&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(&threshold)->default_value(std::numeric_limits::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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/utilities/Rips_complex/ripscomplex.md b/utilities/Rips_complex/ripscomplex.md new file mode 100644 index 00000000..4291fae7 --- /dev/null +++ b/utilities/Rips_complex/ripscomplex.md @@ -0,0 +1,49 @@ + + +# Rips complex # + +## rips_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* of a Rips complex defined on a set of input points, using Euclidean distance. The output diagram contains one bar per line, written with the convention: + +`p dim birth death` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). + +**Usage** + +`rips_persistence [options] ` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction. +* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +Beware: this program may use a lot of RAM and take a lot of time if `max-edge-length` is set to a large value. + +**Example 1 with Z/2Z coefficients** + +`rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2` + +**Example 2 with Z/3Z coefficients** + +`rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3` + + +## rips_distance_matrix_persistence ## + +Same as `rips_persistence` but taking a distance matrix as input. + +**Usage** + +`rips_persistence [options] ` + +where +`` is the path to the file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'. + +**Example** + +`rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` diff --git a/utilities/Witness_complex/CMakeLists.txt b/utilities/Witness_complex/CMakeLists.txt new file mode 100644 index 00000000..125a41ff --- /dev/null +++ b/utilities/Witness_complex/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 2.6) +project(Witness_complex_utilities) + +# 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_strong_witness_persistence strong_witness_persistence.cpp ) + target_link_libraries(Witness_complex_strong_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + add_executable ( Witness_complex_weak_witness_persistence weak_witness_persistence.cpp ) + target_link_libraries(Witness_complex_weak_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(Witness_complex_strong_witness_persistence ${TBB_LIBRARIES}) + target_link_libraries(Witness_complex_weak_witness_persistence ${TBB_LIBRARIES}) + endif() + + add_test(NAME Witness_complex_strong_test_torus_persistence + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5") + add_test(NAME Witness_complex_weak_test_torus_persistence + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5") + + install(TARGETS Witness_complex_strong_witness_persistence DESTINATION bin) + install(TARGETS Witness_complex_weak_witness_persistence DESTINATION bin) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0) diff --git a/utilities/Witness_complex/strong_witness_persistence.cpp b/utilities/Witness_complex/strong_witness_persistence.cpp new file mode 100644 index 00000000..2fba631b --- /dev/null +++ b/utilities/Witness_complex/strong_witness_persistence.cpp @@ -0,0 +1,156 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include // infinity + +using K = CGAL::Epick_d; +using Point_d = K::Point_d; + +using Point_vector = std::vector; +using Strong_witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex; +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; + +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 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 (decomment one of the following two lines) + // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point, + 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(&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::infinity(); + visible.add_options()("help,h", "produce help message")("landmarks,l", po::value(&nbL), + "Number of landmarks to choose from the point cloud.")( + "output-file,o", po::value(&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(&max_squared_alpha)->default_value(default_alpha), + "Maximal squared relaxation parameter.")( + "field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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(&dim_max)->default_value(std::numeric_limits::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/utilities/Witness_complex/weak_witness_persistence.cpp b/utilities/Witness_complex/weak_witness_persistence.cpp new file mode 100644 index 00000000..23fa93aa --- /dev/null +++ b/utilities/Witness_complex/weak_witness_persistence.cpp @@ -0,0 +1,156 @@ +/* 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 . + */ + +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include // infinity + +using K = CGAL::Epick_d; +using Point_d = K::Point_d; + +using Point_vector = std::vector; +using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex; +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; + +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 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 (decomment one of the following two lines) + // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point, + 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(&file_name), + "Name of file containing a point set in off format."); + + Filtration_value default_alpha = std::numeric_limits::infinity(); + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")("landmarks,l", po::value(&nbL), + "Number of landmarks to choose from the point cloud.")( + "output-file,o", po::value(&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(&max_squared_alpha)->default_value(default_alpha), + "Maximal squared relaxation parameter.")( + "field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::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(&dim_max)->default_value(std::numeric_limits::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/utilities/Witness_complex/witnesscomplex.md b/utilities/Witness_complex/witnesscomplex.md new file mode 100644 index 00000000..2341759b --- /dev/null +++ b/utilities/Witness_complex/witnesscomplex.md @@ -0,0 +1,66 @@ + + +# Witness complex # + + +For more details about the witness complex, please read the [user manual of the package](/doc/latest/group__witness__complex.html). + +## weak_witness_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* of a Weak witness complex defined on a set of input points. +The output diagram contains one bar per line, written with the convention: + +`p dim birth death` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, +and `p` is the characteristic of the field *Z/pZ* used for homology coefficients. + +**Usage** + +`weak_witness_persistence [options] ` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud. +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout. +* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +* `-d [ --cpx-dimension ]` (default = 2147483647) Maximal dimension of the weak witness complex we want to compute. + +**Example** + +`weak_witness_persistence data/points/tore3D_1307.off -l 20 -a 0.5 -m 0.006` + +N.B.: output is random as the 20 landmarks are chosen randomly. + + +## strong_witness_persistence ## + +This program computes the persistent homology with coefficient field *Z/pZ* of a Strong witness complex defined on a set of input points. +The output diagram contains one bar per line, written with the convention: + +`p dim birth death` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, +and `p` is the characteristic of the field *Z/pZ* used for homology coefficients. + +**Usage** + +`strong_witness_persistence [options] ` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud. +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout. +* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +* `-d [ --cpx-dimension ]` (default = 2147483647) Maximal dimension of the weak witness complex we want to compute. + +**Example** + +`strong_witness_persistence data/points/tore3D_1307.off -l 20 -a 0.5 -m 0.06` + +N.B.: output is random as the 20 landmarks are chosen randomly. diff --git a/utilities/common/README b/utilities/common/README deleted file mode 100644 index dc841521..00000000 --- a/utilities/common/README +++ /dev/null @@ -1,19 +0,0 @@ -======================= off_file_from_shape_generator ================================== - -Example of use : - -*** on|in sphere|cube|curve|torus|klein generator - -./off_file_from_shape_generator on sphere onSphere.off 1000 3 15.2 - - => generates a onSphere.off file with 1000 points randomized on a sphere of dimension 3 and radius 15.2 - -./off_file_from_shape_generator in sphere inSphere.off 100 2 - - => generates a inSphere.off file with 100 points randomized in a sphere of dimension 2 (circle) and radius 1.0 (default) - -./off_file_from_shape_generator in cube inCube.off 10000 3 5.8 - - => generates a inCube.off file with 10000 points randomized in a cube of dimension 3 and side 5.8 - -!! Warning: hypegenerator on cube is not available !! diff --git a/utilities/common/pointsetgenerator.md b/utilities/common/pointsetgenerator.md new file mode 100644 index 00000000..284715d4 --- /dev/null +++ b/utilities/common/pointsetgenerator.md @@ -0,0 +1,33 @@ + + +# common # + +## off_file_from_shape_generator ## + +Generates a pointset and save it in an OFF file. Command-line is: + +``` +off_file_from_shape_generator on|in sphere|cube|curve|torus|klein ... +``` + +Warning: "on cube" generator is not available! + +**Examples** + +``` +off_file_from_shape_generator on sphere onSphere.off 1000 3 15.2 +``` + +* Generates an onSphere.off file with 1000 points randomized on a sphere of dimension 3 and radius 15.2. + +``` +off_file_from_shape_generator in sphere inSphere.off 100 2 +``` + +* Generates an inSphere.off file with 100 points randomized in a sphere of dimension 2 (circle) and radius 1.0 (default). + +``` +off_file_from_shape_generator in cube inCube.off 10000 3 5.8 +``` + +* Generates a inCube.off file with 10000 points randomized in a cube of dimension 3 and side 5.8. -- cgit v1.2.3