From 0586a149b5bb3a4b65b63b2ab7d3ecdd9682ee1b Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 20 Feb 2018 16:03:52 +0000 Subject: tests and utils fix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3253 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5786a8a7e4b16750f29fac99ca61926158542cfd --- .../concept/SimplicialComplexForCech.h | 66 +++++ .../example_one_skeleton_cech_from_points.cpp | 5 +- src/Cech_complex/include/Miniball/Miniball.hpp | 5 + src/Cech_complex/include/gudhi/Cech_complex.h | 18 +- .../include/gudhi/Cech_complex_blocker.h | 23 +- src/Cech_complex/test/CMakeLists.txt | 15 ++ src/Cech_complex/test/README | 12 + src/Cech_complex/test/test_cech_complex.cpp | 274 +++++++++++++++++++++ src/Cech_complex/utilities/CMakeLists.txt | 14 ++ src/Cech_complex/utilities/cech_persistence.cpp | 136 ++++++++++ src/Cech_complex/utilities/cechcomplex.md | 33 +++ 11 files changed, 587 insertions(+), 14 deletions(-) create mode 100644 src/Cech_complex/concept/SimplicialComplexForCech.h create mode 100644 src/Cech_complex/test/CMakeLists.txt create mode 100644 src/Cech_complex/test/README create mode 100644 src/Cech_complex/test/test_cech_complex.cpp create mode 100644 src/Cech_complex/utilities/CMakeLists.txt create mode 100644 src/Cech_complex/utilities/cech_persistence.cpp create mode 100644 src/Cech_complex/utilities/cechcomplex.md (limited to 'src/Cech_complex') diff --git a/src/Cech_complex/concept/SimplicialComplexForCech.h b/src/Cech_complex/concept/SimplicialComplexForCech.h new file mode 100644 index 00000000..1954a703 --- /dev/null +++ b/src/Cech_complex/concept/SimplicialComplexForCech.h @@ -0,0 +1,66 @@ +/* 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) 2018 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 . + */ + +#ifndef CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_ +#define CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_ + +namespace Gudhi { + +namespace cech_complex { + +/** \brief The concept SimplicialComplexForCech describes the requirements for a type to implement a simplicial + * complex, that can be created from a `Cech_complex`. + */ +struct SimplicialComplexForCech { + /** Handle to specify a simplex. */ + typedef unspecified Simplex_handle; + /** Handle to specify a vertex. Must be a non-negative integer. */ + typedef unspecified Vertex_handle; + /** Handle to specify the simplex filtration value. */ + typedef unspecified Filtration_value; + + /** Assigns the 'simplex' with the given 'filtration' value. */ + int assign_filtration(Simplex_handle simplex, Filtration_value filtration); + + /** \brief Returns a range over vertices of a given + * simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle const & simplex); + + /** \brief Inserts a given `Gudhi::ProximityGraph` in the simplicial complex. */ + template + void insert_graph(const ProximityGraph& proximity_graph); + + /** \brief Expands the simplicial complex containing only its one skeleton until a given maximal dimension. + * expansion can be blocked by the blocker oracle. */ + template< typename Blocker > + void expansion_with_blockers(int max_dim, Blocker block_simplex); + + /** Returns the number of vertices in the simplicial complex. */ + std::size_t num_vertices(); + +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ diff --git a/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp b/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp index 9b03616c..73679716 100644 --- a/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp +++ b/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp @@ -27,11 +27,12 @@ #include #include #include +#include #include // for std::numeric_limits int main() { // Type definitions - using Point_cloud = std::vector>; + using Point_cloud = std::vector>; using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -52,7 +53,7 @@ int main() { Cech_complex cech_complex_from_points(points, threshold, Gudhi::Euclidean_distance()); Simplex_tree stree; - cech_complex_from_points.create_complex(stree, 2); + cech_complex_from_points.create_complex(stree, 3); // ---------------------------------------------------------------------------- // Display information about the one skeleton Rips complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/include/Miniball/Miniball.hpp b/src/Cech_complex/include/Miniball/Miniball.hpp index cb76c534..a42d62a7 100644 --- a/src/Cech_complex/include/Miniball/Miniball.hpp +++ b/src/Cech_complex/include/Miniball/Miniball.hpp @@ -23,6 +23,9 @@ // CH-8092 Zuerich, Switzerland // http://www.inf.ethz.ch/personal/gaertner +#ifndef MINIBALL_HPP_ +#define MINIBALL_HPP_ + #include #include #include @@ -513,3 +516,5 @@ namespace Miniball { } } // end Namespace Miniball + +#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 94939105..e847c97f 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -28,7 +28,6 @@ #include // for Gudhi::cech_complex::Cech_blocker #include -#include #include // for std::size_t #include // for exception management @@ -40,7 +39,7 @@ namespace cech_complex { * \class Cech_complex * \brief Cech complex data structure. * - * \ingroup Cech_complex + * \ingroup cech_complex * * \details * The data structure is a proximity graph, containing edges when the edge length is less or equal @@ -65,10 +64,9 @@ class Cech_complex { * @param[in] points Range of points. * @param[in] threshold Rips value. * @param[in] distance distance function that returns a `Filtration_value` from 2 given points. - * @exception std::invalid_argument In debug mode, if `points.size()` returns a value ≤ 0. * * \tparam ForwardPointRange must be a range for which `.size()`, `.begin()` and `.end()` methods return input - * iterators on a point. A point must have a `.size()` method available. + * iterators on a point. `.begin()` and `.end()` methods are required for a point. * * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`. @@ -77,12 +75,10 @@ class Cech_complex { Cech_complex(const ForwardPointRange& points, Filtration_value threshold, Distance distance) : threshold_(threshold), point_cloud_(points) { - GUDHI_CHECK(points.size() > 0, - std::invalid_argument("Cech_complex::create_complex - point cloud is empty")); - dimension_ = points.begin()->size(); + dimension_ = points.begin()->end() - points.begin()->begin(); cech_skeleton_graph_ = Gudhi::compute_proximity_graph(point_cloud_, - threshold_, - distance); + threshold_, + distance); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal @@ -116,10 +112,10 @@ class Cech_complex { } /** @param[in] vertex Point position in the range. - * @return Threshold value given at construction. + * @return A const iterator on the point. * @exception std::out_of_range In debug mode, if point position in the range is out. */ - typename ForwardPointRange::const_iterator point(std::size_t vertex) const { + typename ForwardPointRange::const_iterator point_iterator(std::size_t vertex) const { GUDHI_CHECK((point_cloud_.begin() + vertex) < point_cloud_.end(), std::out_of_range("Cech_complex::point - simplicial complex is not empty")); return (point_cloud_.begin() + vertex); diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 25fab909..f8738be0 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -39,6 +39,20 @@ namespace cech_complex { template class Cech_complex; +/** \internal + * \class Cech_blocker + * \brief Cech complex blocker. + * + * \ingroup cech_complex + * + * \details + * Cech blocker is an oracle constructed from a Cech_complex and a simplicial complex. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Simplex_handle` and `Filtration_value` type definition, + * `simplex_vertex_range(Simplex_handle sh)`and `assign_filtration(Simplex_handle sh, Filtration_value filt)` methods. + * + * \tparam ForwardPointRange is required by the pointer on Chech_complex for type definition. + */ template class Cech_blocker { private: @@ -52,10 +66,15 @@ class Cech_blocker { using Cech_complex = Gudhi::cech_complex::Cech_complex; public: + /** \internal \brief Cech complex blocker operator() - the oracle - assigns the filtration value from the simplex + * radius and returns if the simplex expansion must be blocked. + * \param[in] sh The Simplex_handle. + * \return true if the simplex radius is greater than the Cech_complex threshold*/ bool operator()(Simplex_handle sh) { Point_cloud points; for (auto vertex : simplicial_complex_.simplex_vertex_range(sh)) { - points.push_back(*(cc_ptr_->point(vertex))); + points.push_back(Point(cc_ptr_->point_iterator(vertex)->begin(), + cc_ptr_->point_iterator(vertex)->end())); #ifdef DEBUG_TRACES std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES @@ -68,6 +87,8 @@ class Cech_blocker { simplicial_complex_.assign_filtration(sh, radius); return (radius > cc_ptr_->threshold()); } + + /** \internal \brief Cech complex blocker constructor. */ Cech_blocker(SimplicialComplexForCech& simplicial_complex, Cech_complex* cc_ptr) : simplicial_complex_(simplicial_complex), cc_ptr_(cc_ptr) { diff --git a/src/Cech_complex/test/CMakeLists.txt b/src/Cech_complex/test/CMakeLists.txt new file mode 100644 index 00000000..8db51173 --- /dev/null +++ b/src/Cech_complex/test/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_tests) + +include(GUDHI_test_coverage) + +add_executable ( Cech_complex_test_unit test_cech_complex.cpp ) +target_link_libraries(Cech_complex_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES}) +endif() + +# Do not forget to copy test files in current binary dir +file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +gudhi_add_coverage_test(Cech_complex_test_unit) diff --git a/src/Cech_complex/test/README b/src/Cech_complex/test/README new file mode 100644 index 00000000..adf704f7 --- /dev/null +++ b/src/Cech_complex/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./Cech_complex_test_unit --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp new file mode 100644 index 00000000..aa42d322 --- /dev/null +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -0,0 +1,274 @@ +/* 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) 2018 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 . + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "cech_complex" +#include + +#include // float comparison +#include +#include +#include +#include // std::max + +#include +// to construct Cech_complex from a OFF file of points +#include +#include +#include +#include + +#include + +// Type definitions +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector; +using Point_cloud = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; +using Cech_complex = Gudhi::cech_complex::Cech_complex; + +using Point_iterator = Point_cloud::const_iterator; +using Coordinate_iterator = Point::const_iterator; +using Min_sphere = Miniball::Miniball>; + +BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double threshold = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech threshold=" << + threshold << "==========" << std::endl; + + Points_off_reader off_reader(off_file_name); + Point_cloud point_cloud = off_reader.get_point_cloud(); + Cech_complex cech_complex_from_file(point_cloud, threshold, Gudhi::Euclidean_distance()); + + std::size_t i = 0; + for (; i < point_cloud.size(); i++) { + BOOST_CHECK(point_cloud[i] == *(cech_complex_from_file.point_iterator(i))); + } +#ifdef GUDHI_DEBUG + BOOST_CHECK_THROW (cech_complex_from_file.point_iterator(i+1), std::out_of_range); +#endif // GUDHI_DEBUG + + const int DIMENSION_1 = 1; + Simplex_tree st; + cech_complex_from_file.create_complex(st, DIMENSION_1); + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION_1); + + const int NUMBER_OF_VERTICES = 7; + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 18); + + // Check filtration values of vertices is 0.0 + for (auto f_simplex : st.skeleton_simplex_range(0)) { + BOOST_CHECK(st.filtration(f_simplex) == 0.0); + } + + // Check filtration values of edges + for (auto f_simplex : st.skeleton_simplex_range(DIMENSION_1)) { + if (DIMENSION_1 == st.dimension(f_simplex)) { + std::vector vp; + std::cout << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << ","; + vp.push_back(off_reader.get_point_cloud().at(vertex)); + } + std::cout << ") - distance =" << Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)) << + " - filtration =" << st.filtration(f_simplex) << std::endl; + BOOST_CHECK(vp.size() == 2); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), Gudhi::Euclidean_distance()(vp.at(0), vp.at(1))); + } + } + + const int DIMENSION_2 = 2; + Simplex_tree st2; + cech_complex_from_file.create_complex(st2, DIMENSION_2); + std::cout << "st2.dimension()=" << st2.dimension() << std::endl; + BOOST_CHECK(st2.dimension() == DIMENSION_2); + + std::cout << "st2.num_vertices()=" << st2.num_vertices() << std::endl; + BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; + BOOST_CHECK(st2.num_simplices() == 23); + + Point_cloud points012; + for (std::size_t vertex = 0; vertex <= 2; vertex++) { + points012.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), + cech_complex_from_file.point_iterator(vertex)->end())); + } + Min_sphere ms012(cech_complex_from_file.dimension(), points012.begin(),points012.end()); + + Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); + std::cout << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); + + Point_cloud points456; + for (std::size_t vertex = 4; vertex <= 6; vertex++) { + points456.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), + cech_complex_from_file.point_iterator(vertex)->end())); + } + Min_sphere ms456(cech_complex_from_file.dimension(), points456.begin(),points456.end()); + + Simplex_tree::Filtration_value f456 = st2.filtration(st2.find({4, 5, 6})); + std::cout << "f456= " << f456 << " | ms456_radius= " << std::sqrt(ms456.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f456, std::sqrt(ms456.squared_radius())); + + const int DIMENSION_3 = 3; + Simplex_tree st3; + cech_complex_from_file.create_complex(st3, DIMENSION_3); + std::cout << "st3.dimension()=" << st3.dimension() << std::endl; + BOOST_CHECK(st3.dimension() == DIMENSION_3); + + std::cout << "st3.num_vertices()=" << st3.num_vertices() << std::endl; + BOOST_CHECK(st3.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st3.num_simplices()=" << st3.num_simplices() << std::endl; + BOOST_CHECK(st3.num_simplices() == 24); + + Point_cloud points0123; + for (std::size_t vertex = 0; vertex <= 3; vertex++) { + points0123.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), + cech_complex_from_file.point_iterator(vertex)->end())); + } + Min_sphere ms0123(cech_complex_from_file.dimension(), points0123.begin(),points0123.end()); + + Simplex_tree::Filtration_value f0123 = st3.filtration(st3.find({0, 1, 2, 3})); + std::cout << "f0123= " << f0123 << " | ms0123_radius= " << std::sqrt(ms0123.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f0123, std::sqrt(ms0123.squared_radius())); + + + + Point_cloud points01; + for (std::size_t vertex = 0; vertex <= 1; vertex++) { + points01.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), + cech_complex_from_file.point_iterator(vertex)->end())); + } + Min_sphere ms01(cech_complex_from_file.dimension(), points01.begin(),points01.end()); + + Simplex_tree::Filtration_value f01 = st2.filtration(st2.find({0, 1})); + std::cout << "f01= " << f01 << " | ms01_radius= " << std::sqrt(ms01.squared_radius()) << std::endl; + +} + +BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Point_cloud points; + std::vector coords = { 0.0, 0.0, 0.0, 1.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 0.0, 1.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 1.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 1.0, 0.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of a Cech complex from the list of points + // ---------------------------------------------------------------------------- + Cech_complex cech_complex_from_points(points, 2.0, Gudhi::Euclidean_distance()); + + std::cout << "========== cech_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + cech_complex_from_points.create_complex(st, DIMENSION); + + // Another way to check num_simplices + std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : st.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 15); + + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION); + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == 4); + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; + switch (st.dimension(f_simplex)) { + case 0: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0); + break; + case 1: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 1.41421, .00001); + break; + case 2: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.816497, .00001); + break; + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.866025, .00001); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double threshold = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech threshold=" << + threshold << "==========" << std::endl; + + Gudhi::Points_off_reader off_reader(off_file_name); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + Simplex_tree stree; + std::vector simplex = {0, 1, 2}; + stree.insert_simplex_and_subfaces(simplex); + std::cout << "Check exception throw in debug mode" << std::endl; + // throw excpt because stree is not empty + BOOST_CHECK_THROW (cech_complex_from_file.create_complex(stree, 1), std::invalid_argument); +} +#endif diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt new file mode 100644 index 00000000..a4f89d2c --- /dev/null +++ b/src/Cech_complex/utilities/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_utilities) + +add_executable(cech_persistence cech_persistence.cpp) +target_link_libraries(cech_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(cech_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Cech_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 cech_persistence DESTINATION bin) diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp new file mode 100644 index 00000000..e93596d4 --- /dev/null +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -0,0 +1,136 @@ +/* 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 Point = std::vector; +using Point_cloud = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; +using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; + +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); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Cech complex in a Simplex Tree + Simplex_tree simplex_tree; + + cech_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 Cech complex construction.")( + "cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Cech 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 Cech complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Cech_complex/utilities/cechcomplex.md b/src/Cech_complex/utilities/cechcomplex.md new file mode 100644 index 00000000..6330727a --- /dev/null +++ b/src/Cech_complex/utilities/cechcomplex.md @@ -0,0 +1,33 @@ + + +# Cech complex # + +## cech_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* of a Cech 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** + +`cech_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 Cech complex construction. +* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Cech 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** + +`cech_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2` + +**Example 2 with Z/3Z coefficients** + +`cech_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3` -- cgit v1.2.3