diff options
Diffstat (limited to 'src/Rips_complex')
18 files changed, 1217 insertions, 101 deletions
diff --git a/src/Rips_complex/concept/SimplicialComplexForRips.h b/src/Rips_complex/concept/SimplicialComplexForRips.h index 7dab0615..3c5acecf 100644 --- a/src/Rips_complex/concept/SimplicialComplexForRips.h +++ b/src/Rips_complex/concept/SimplicialComplexForRips.h @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2016 INRIA + * 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 diff --git a/src/Rips_complex/doc/COPYRIGHT b/src/Rips_complex/doc/COPYRIGHT index 594b7d03..2c31a0d6 100644 --- a/src/Rips_complex/doc/COPYRIGHT +++ b/src/Rips_complex/doc/COPYRIGHT @@ -4,7 +4,7 @@ computational topology. Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau -Copyright (C) 2015 INRIA +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 diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 124dfec9..712d3b6e 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau * - * Copyright (C) 2016 INRIA + * 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 @@ -29,28 +29,41 @@ namespace rips_complex { /** \defgroup rips_complex Rips complex * - * \author Clément Maria, Pawel Dlotko, Vincent Rouvreau + * \author Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse * * @{ * * \section ripsdefinition Rips complex definition * - * Rips_complex - * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> is a - * one skeleton graph that allows to construct a - * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> - * from it. - * The input can be a point cloud with a given distance function, or a distance matrix. - * - * The filtration value of each edge is computed from a user-given distance function, or directly from the distance - * matrix. + * The Vietoris-Rips complex + * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> + * is an abstract simplicial complex + * defined on a finite metric space, where each simplex corresponds to a subset + * of point whose diameter is smaller that some given threshold. + * Varying the threshold, we can also see the Rips complex as a filtration of + * the \f$(n-1)-\f$dimensional simplex, where the filtration value of each + * simplex is the diameter of the corresponding subset of points. + * + * This filtered complex is most often used as an approximation of the + * Čech complex. After rescaling (Rips using the length of the edges and Čech + * the half-length), they share the same 1-skeleton and are multiplicatively + * 2-interleaved or better. While it is slightly bigger, it is also much + * easier to compute. + * + * The number of simplices in the full Rips complex is exponential in the + * number of vertices, it is thus usually restricted, by excluding all the + * simplices with filtration value larger than some threshold, and keeping only + * the dim_max-skeleton. + * + * In order to build this complex, the algorithm first builds the graph. + * The filtration value of each edge is computed from a user-given distance + * function, or directly read from the distance matrix. + * In a second step, this graph is inserted in a simplicial complex, which then + * gets expanded to a flag complex. * - * All edges that have a filtration value strictly greater than a given threshold value are not inserted into - * the complex. + * The input can be given as a range of points and a distance function, or as a + * distance matrix. * - * When creating a simplicial complex from this one skeleton graph, Rips inserts the one skeleton graph into the data - * structure, and then expands the simplicial complex when required. - * * Vertex name correspond to the index of the point in the given range (aka. the point cloud). * * \image html "rips_complex_representation.png" "Rips-complex one skeleton graph representation" @@ -61,7 +74,36 @@ namespace rips_complex { * * If the Rips_complex interfaces are not detailed enough for your need, please refer to * <a href="_persistent_cohomology_2rips_persistence_step_by_step_8cpp-example.html"> - * rips_persistence_step_by_step.cpp</a> example, where the graph construction over the Simplex_tree is more detailed. + * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and + * the Simplex_tree are more detailed. + * + * \section sparserips Sparse Rips complex + * + * Even truncated in filtration value and dimension, the Rips complex remains + * quite large. However, it is possible to approximate it by a much smaller + * filtered simplicial complex (linear size, with constants that depend on + * ε and the doubling dimension of the space) that is + * \f$(1+O(\epsilon))-\f$interleaved with it (in particular, their persistence + * diagrams are at log-bottleneck distance at most \f$O(\epsilon)\f$). + * + * The sparse Rips filtration was introduced by Don Sheehy + * \cite sheehy13linear. We are using the version described in + * \cite buchet16efficient (except that we multiply all filtration values + * by 2, to match the usual Rips complex), which proves a + * \f$\frac{1+\epsilon}{1-\epsilon}\f$-interleaving, although in practice the + * error is usually smaller. + * A more intuitive presentation of the idea is available in + * \cite cavanna15geometric, and in a video \cite cavanna15visualizing. + * + * The interface of `Sparse_rips_complex` is similar to the one for the usual + * `Rips_complex`, except that one has to specify the approximation factor, and + * there is no option to limit the maximum filtration value (the way the + * approximation is done means that larger filtration values are much cheaper + * to handle than low filtration values, so the gain would be too small). + * + * Theoretical guarantees are only available for \f$\epsilon<1\f$. The + * construction accepts larger values of ε, and the size of the complex + * keeps decreasing, but there is no guarantee on the quality of the result. * * \section ripspointsdistance Point cloud and distance function * @@ -104,6 +146,24 @@ namespace rips_complex { * * \include Rips_complex/full_skeleton_rips_for_doc.txt * + * + * \subsection sparseripspointscloudexample Example of a sparse Rips from a point cloud + * + * This example builds the full sparse Rips of a set of 2D Euclidean points, then prints some minimal + * information about the complex. + * + * \include Rips_complex/example_sparse_rips.cpp + * + * When launching: + * + * \code $> ./Rips_complex_example_sparse + * \endcode + * + * the program output may be (the exact output varies from one run to the next): + * + * \code Sparse Rips complex is of dimension 2 - 19 simplices - 7 vertices. + * \endcode + * * * * \section ripsdistancematrix Distance matrix @@ -146,8 +206,33 @@ namespace rips_complex { * * \include Rips_complex/full_skeleton_rips_for_doc.txt * - * \copyright GNU General Public License v3. - * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim + * + * \section ripscorrelationematrix Correlation matrix + * + * Analogously to the case of distance matrix, Rips complexes can be also constructed based on correlation matrix. + * Given a correlation matrix M, comportment-wise 1-M is a distance matrix. + * This example builds the one skeleton graph from the given corelation matrix and threshold value. + * Then it creates a `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Rips_complex/example_one_skeleton_rips_from_correlation_matrix.cpp + * + * When launching: + * + * \code $> ./example_one_skeleton_from_correlation_matrix + * \endcode + * + * the program output is: + * + * \include Rips_complex/one_skeleton_rips_from_correlation_matrix_for_doc.txt + * + * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes + * construction from correlation matrices. + * + * @warning As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical + * tool will not work properly, this is a known issue. + * */ /** @} */ // end defgroup rips_complex diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index b854b1c9..e7772bdb 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -1,11 +1,8 @@ -cmake_minimum_required(VERSION 2.6) project(Rips_complex_examples) # Point cloud add_executable ( Rips_complex_example_from_off example_rips_complex_from_off_file.cpp ) -add_executable ( Rips_complex_example_from_fvecs example_rips_complex_from_fvecs.cpp ) - add_executable ( Rips_complex_example_one_skeleton_from_points example_one_skeleton_rips_from_points.cpp ) # Distance matrix @@ -13,17 +10,29 @@ add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_ add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp ) +# Sparse rips from points +add_executable ( Rips_complex_example_sparse example_sparse_rips.cpp ) + +# Correlation matrix +add_executable ( Rips_complex_example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp ) + if (TBB_FOUND) target_link_libraries(Rips_complex_example_from_off ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_sparse ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) endif() add_test(NAME Rips_complex_example_one_skeleton_from_points COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_points>) add_test(NAME Rips_complex_example_one_skeleton_from_distance_matrix COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_distance_matrix>) +add_test(NAME Rips_complex_example_sparse + COMMAND $<TARGET_FILE:Rips_complex_example_sparse>) +add_test(NAME Rips_complex_example_one_skeleton_rips_from_correlation_matrix + COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_rips_from_correlation_matrix>) add_test(NAME Rips_complex_example_from_off_doc_12_1 COMMAND $<TARGET_FILE:Rips_complex_example_from_off> "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "12.0" "1" "${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt") @@ -59,3 +68,4 @@ install(TARGETS Rips_complex_example_from_off DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_points DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_distance_matrix DESTINATION bin) install(TARGETS Rips_complex_example_from_csv_distance_matrix DESTINATION bin) +install(TARGETS Rips_complex_example_one_skeleton_rips_from_correlation_matrix DESTINATION bin) diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp new file mode 100644 index 00000000..05bacb9f --- /dev/null +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp @@ -0,0 +1,81 @@ +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for std::numeric_limits + +int main() { + // Type definitions + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + using Distance_matrix = std::vector<std::vector<Filtration_value>>; + + // User defined correlation matrix is: + // |1 0.06 0.23 0.01 0.89| + // |0.06 1 0.74 0.01 0.61| + // |0.23 0.74 1 0.72 0.03| + // |0.01 0.01 0.72 1 0.7 | + // |0.89 0.61 0.03 0.7 1 | + + Distance_matrix correlations; + correlations.push_back({}); + correlations.push_back({0.06}); + correlations.push_back({0.23, 0.74}); + correlations.push_back({0.01, 0.01, 0.72}); + correlations.push_back({0.89, 0.61, 0.03, 0.7}); + + // ---------------------------------------------------------------------------- + // Convert correlation matrix to a distance matrix: + // ---------------------------------------------------------------------------- + double threshold = 0; + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + // Here we check if our data comes from corelation matrix. + if ((correlations[i][j] < -1) || (correlations[i][j] > 1)) { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate.\n"; + throw "The input matrix is not a correlation matrix. The program will now terminate.\n"; + } + correlations[i][j] = 1 - correlations[i][j]; + // Here we make sure that we will get the treshold value equal to maximal + // distance in the matrix. + if (correlations[i][j] > threshold) threshold = correlations[i][j]; + } + } + + //----------------------------------------------------------------------------- + // Now the correlation matrix is a distance matrix and can be processed further. + //----------------------------------------------------------------------------- + Distance_matrix distances = correlations; + + Rips_complex rips_complex_from_points(distances, threshold); + + Simplex_tree stree; + rips_complex_from_points.create_complex(stree, 1); + // ---------------------------------------------------------------------------- + // Display information about the one skeleton Rips complex. Note that + // the filtration displayed here comes from the distance matrix computed + // above, which is 1 - initial correlation matrix. Only this way, we obtain + // a complex with filtration. If a correlation matrix is used instead, we would + // have a reverse filtration (i.e. filtration of boundary of each simplex S + // is greater or equal to the filtration of S). + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + + return 0; +} diff --git a/src/Rips_complex/example/example_rips_complex_from_fvecs.cpp b/src/Rips_complex/example/example_rips_complex_from_fvecs.cpp deleted file mode 100644 index 40aab5dc..00000000 --- a/src/Rips_complex/example/example_rips_complex_from_fvecs.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include <gudhi/Rips_complex.h> -#include <gudhi/Simplex_tree.h> -//#include <gudhi/Fake_simplex_tree.h> -//#include <gudhi/Sb_wrapper.h> -#include <gudhi/distance_functions.h> -#include <gudhi/Points_fvecs_reader.h> - - -#include <CGAL/Epick_d.h> - -#include <iostream> -#include <string> -#include <vector> -#include <ctime> - -void usage(int nbArgs, char * const progName) { - std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; - std::cerr << "Usage: " << progName << " filename.fvecs threshold dim_max [ouput_file.txt]\n"; - std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.fvecs 60.0\n"; - exit(-1); // ----- >> -} - -int main(int argc, char **argv) { - if (argc != 4) usage(argc, (argv[0] - 1)); - - std::string file_name(argv[1]); - double threshold = atof(argv[2]); - int dim_max = atoi(argv[3]); - - // Type definitions - using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; - using Point = typename K::Point_d; - using Simplex_tree = Gudhi::Simplex_tree<>; - //using Simplex_tree = Gudhi::Fake_simplex_tree; - //using Simplex_tree = Gudhi::Sb_wrapper; - using Filtration_value = Simplex_tree::Filtration_value; - using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; - using Point_vector = std::vector<Point>; - - // ---------------------------------------------------------------------------- - // Init of a Rips complex from an fvecs file - // ---------------------------------------------------------------------------- - Point_vector point_vector; - Gudhi::load_points_from_fvecs_file<K>(file_name, std::back_insert_iterator< Point_vector >(point_vector)); - - Rips_complex rips_complex_from_file(point_vector, threshold, Gudhi::Euclidean_distance()); - - - Simplex_tree stree; - - clock_t start, end; - start = clock(); - rips_complex_from_file.create_complex(stree, dim_max); - end = clock(); - - std::cout << "Strong witness complex took "<< static_cast<double>(end - start) / CLOCKS_PER_SEC << " s." << std::endl; - //std::cout << "Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices." << std::endl; - - return 0; -} diff --git a/src/Rips_complex/example/example_sparse_rips.cpp b/src/Rips_complex/example/example_sparse_rips.cpp new file mode 100644 index 00000000..1c95b48c --- /dev/null +++ b/src/Rips_complex/example/example_sparse_rips.cpp @@ -0,0 +1,30 @@ +#include <gudhi/Sparse_rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <vector> + +int main() { + using Point = std::vector<double>; + using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; + using Filtration_value = Simplex_tree::Filtration_value; + using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>; + + Point points[] = {{1.0, 1.0}, {7.0, 0.0}, {4.0, 6.0}, {9.0, 6.0}, {0.0, 14.0}, {2.0, 19.0}, {9.0, 17.0}}; + + // ---------------------------------------------------------------------------- + // Init from Euclidean points + // ---------------------------------------------------------------------------- + double epsilon = 2; // very rough, no guarantees + Sparse_rips sparse_rips(points, Gudhi::Euclidean_distance(), epsilon); + + Simplex_tree stree; + sparse_rips.create_complex(stree, 10); + + // ---------------------------------------------------------------------------- + // Display information about the complex + // ---------------------------------------------------------------------------- + std::cout << "Sparse Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; +} diff --git a/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt new file mode 100644 index 00000000..640d7083 --- /dev/null +++ b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt @@ -0,0 +1,17 @@ +Rips complex is of dimension 1 - 15 simplices - 5 vertices. +Iterator on Rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 4 0 ) -> [0.11] + ( 2 1 ) -> [0.26] + ( 3 2 ) -> [0.28] + ( 4 3 ) -> [0.3] + ( 4 1 ) -> [0.39] + ( 2 0 ) -> [0.77] + ( 1 0 ) -> [0.94] + ( 4 2 ) -> [0.97] + ( 3 0 ) -> [0.99] + ( 3 1 ) -> [0.99] diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h index 1e4b76a7..f0fe57f4 100644 --- a/src/Rips_complex/include/gudhi/Rips_complex.h +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau * - * Copyright (C) 2016 INRIA + * 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 diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h new file mode 100644 index 00000000..4dcc08ed --- /dev/null +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -0,0 +1,175 @@ +/* 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): Marc Glisse + * + * 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SPARSE_RIPS_COMPLEX_H_ +#define SPARSE_RIPS_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/choose_n_farthest_points.h> + +#include <boost/graph/adjacency_list.hpp> +#include <boost/range/metafunctions.hpp> + +#include <vector> + +namespace Gudhi { + +namespace rips_complex { + +// The whole interface is copied on Rips_complex. A redesign should be discussed with all complex creation classes in +// mind. + +/** + * \class Sparse_rips_complex + * \brief Sparse Rips complex data structure. + * + * \ingroup rips_complex + * + * \details + * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered + * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * + * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. + */ +template <typename Filtration_value> +class Sparse_rips_complex { + private: + // TODO(MG): use a different graph where we know we can safely insert in parallel. + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, + boost::property<vertex_filtration_t, Filtration_value>, + boost::property<edge_filtration_t, Filtration_value>> + Graph; + + typedef int Vertex_handle; + + public: + /** \brief Sparse_rips_complex constructor from a list of points. + * + * @param[in] points Range of points. + * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. + * @param[in] epsilon Approximation parameter. epsilon must be positive. + * + */ + template <typename RandomAccessPointRange, typename Distance> + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); + std::vector<Vertex_handle> sorted_points; + std::vector<Filtration_value> params; + auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; + Ker<decltype(dist_fun)> kernel(dist_fun); + subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, + std::back_inserter(sorted_points), std::back_inserter(params)); + compute_sparse_graph(sorted_points, params, dist_fun, epsilon); + } + + /** \brief Sparse_rips_complex constructor from a distance matrix. + * + * @param[in] distance_matrix Range of range of distances. + * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and + * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * distance\_matrix.size().\f$ + * @param[in] epsilon Approximation parameter. epsilon must be positive. + */ + template <typename DistanceMatrix> + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), + [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + + /** \brief Fills the simplicial complex with the sparse Rips graph and + * expands it with all the cliques, stopping at a given maximal dimension. + * + * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. + * + * @param[in] complex the complex to fill + * @param[in] dim_max maximal dimension of the simplicial complex. + * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. + * + */ + template <typename SimplicialComplexForRips> + void create_complex(SimplicialComplexForRips& complex, int dim_max) { + GUDHI_CHECK(complex.num_vertices() == 0, + std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); + + complex.insert_graph(graph_); + complex.expansion(dim_max); + } + + private: + // choose_n_farthest_points wants the distance function in this form... + template <class Distance> + struct Ker { + typedef std::size_t Point_d; // index into point range + Ker(Distance& d) : dist(d) {} + // Despite the name, this is not squared... + typedef Distance Squared_distance_d; + Squared_distance_d& squared_distance_d_object() const { return dist; } + Distance& dist; + }; + + // PointRange must be random access. + template <typename PointRange, typename ParamRange, typename Distance> + void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + const int n = boost::size(points); + graph_.~Graph(); + new (&graph_) Graph(n); + // for(auto v : vertices(g)) // doesn't work :-( + typename boost::graph_traits<Graph>::vertex_iterator v_i, v_e; + for (std::tie(v_i, v_e) = vertices(graph_); v_i != v_e; ++v_i) { + auto v = *v_i; + // This whole loop might not be necessary, leave it until someone investigates if it is safe to remove. + put(vertex_filtration_t(), graph_, v, 0); + } + + // TODO(MG): + // - make it parallel + // - only test near-enough neighbors + for (int i = 0; i < n; ++i) + for (int j = i + 1; j < n; ++j) { + auto&& pi = points[i]; + auto&& pj = points[j]; + auto d = dist(pi, pj); + auto li = params[i]; + auto lj = params[j]; + GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); + Filtration_value alpha; + + // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips + if (d * epsilon <= 2 * lj) + alpha = d; + else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) + alpha = (d - lj / epsilon) * 2; + else + continue; + + add_edge(pi, pj, alpha, graph_); + } + } + + Graph graph_; +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // SPARSE_RIPS_COMPLEX_H_ diff --git a/src/Rips_complex/test/CMakeLists.txt b/src/Rips_complex/test/CMakeLists.txt index 3da9c90d..745d953c 100644 --- a/src/Rips_complex/test/CMakeLists.txt +++ b/src/Rips_complex/test/CMakeLists.txt @@ -1,4 +1,3 @@ -cmake_minimum_required(VERSION 2.6) project(Rips_complex_tests) include(GUDHI_test_coverage) diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp index fc83f5f7..b8b444c9 100644 --- a/src/Rips_complex/test/test_rips_complex.cpp +++ b/src/Rips_complex/test/test_rips_complex.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2016 INRIA Saclay (France) + * 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 @@ -31,23 +31,22 @@ #include <algorithm> // std::max #include <gudhi/Rips_complex.h> +#include <gudhi/Sparse_rips_complex.h> // to construct Rips_complex from a OFF file of points #include <gudhi/Points_off_io.h> #include <gudhi/Simplex_tree.h> #include <gudhi/distance_functions.h> #include <gudhi/reader_utils.h> +#include <gudhi/Unitary_tests_utils.h> // Type definitions using Point = std::vector<double>; using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex<Simplex_tree::Filtration_value>; +using Sparse_rips_complex = Gudhi::rips_complex::Sparse_rips_complex<Simplex_tree::Filtration_value>; using Distance_matrix = std::vector<std::vector<Filtration_value>>; -bool are_almost_the_same(float a, float b) { - return std::fabs(a - b) < std::numeric_limits<float>::epsilon(); -} - BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { // ---------------------------------------------------------------------------- // @@ -92,7 +91,7 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { std::cout << ") - distance =" << Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); - BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), Gudhi::Euclidean_distance()(vp.at(0), vp.at(1))); } } @@ -113,14 +112,14 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { Simplex_tree::Filtration_value f12 = st2.filtration(st2.find({1, 2})); Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); std::cout << "f012= " << f012 << " | f01= " << f01 << " - f02= " << f02 << " - f12= " << f12 << std::endl; - BOOST_CHECK(are_almost_the_same(f012, std::max(f01, std::max(f02,f12)))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::max(f01, std::max(f02,f12))); Simplex_tree::Filtration_value f45 = st2.filtration(st2.find({4, 5})); Simplex_tree::Filtration_value f56 = st2.filtration(st2.find({5, 6})); Simplex_tree::Filtration_value f46 = st2.filtration(st2.find({4, 6})); Simplex_tree::Filtration_value f456 = st2.filtration(st2.find({4, 5, 6})); std::cout << "f456= " << f456 << " | f45= " << f45 << " - f56= " << f56 << " - f46= " << f46 << std::endl; - BOOST_CHECK(are_almost_the_same(f456, std::max(f45, std::max(f56,f46)))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f456, std::max(f45, std::max(f56,f46))); const int DIMENSION_3 = 3; Simplex_tree st3; @@ -140,7 +139,7 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { Simplex_tree::Filtration_value f0123 = st3.filtration(st3.find({0, 1, 2, 3})); std::cout << "f0123= " << f0123 << " | f012= " << f012 << " - f123= " << f123 << " - f013= " << f013 << " - f023= " << f023 << std::endl; - BOOST_CHECK(are_almost_the_same(f0123, std::max(f012, std::max(f123, std::max(f013, f023))))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f0123, std::max(f012, std::max(f123, std::max(f013, f023)))); } @@ -219,12 +218,77 @@ BOOST_AUTO_TEST_CASE(Rips_complex_from_points) { std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; switch (st.dimension(f_simplex)) { case 0: - BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), 0.0)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0); + break; + case 1: + case 2: + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 2.0); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} + +BOOST_AUTO_TEST_CASE(Sparse_rips_complex_from_points) { + // This is a clone of the test above + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_of_points points; + std::vector<double> 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 Rips complex from the list of points + // ---------------------------------------------------------------------------- + // .001 is small enough that we get a deterministic result matching the exact Rips + Sparse_rips_complex sparse_rips(points, Custom_square_euclidean_distance(), .001); + + std::cout << "========== Sparse_rips_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + sparse_rips.create_complex(st, DIMENSION); + + // Another way to check num_simplices + std::cout << "Iterator on Rips 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: case 2: case 3: - BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), 2.0)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 2.0); break; default: BOOST_CHECK(false); // Shall not happen @@ -276,7 +340,7 @@ BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) { } std::cout << ") - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vvh.size() == 2); - BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), distances[vvh.at(0)][vvh.at(1)])); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), distances[vvh.at(0)][vvh.at(1)]); } } @@ -297,14 +361,14 @@ BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) { Simplex_tree::Filtration_value f12 = st2.filtration(st2.find({1, 2})); Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); std::cout << "f012= " << f012 << " | f01= " << f01 << " - f02= " << f02 << " - f12= " << f12 << std::endl; - BOOST_CHECK(are_almost_the_same(f012, std::max(f01, std::max(f02,f12)))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::max(f01, std::max(f02,f12))); Simplex_tree::Filtration_value f45 = st2.filtration(st2.find({4, 5})); Simplex_tree::Filtration_value f56 = st2.filtration(st2.find({5, 6})); Simplex_tree::Filtration_value f46 = st2.filtration(st2.find({4, 6})); Simplex_tree::Filtration_value f456 = st2.filtration(st2.find({4, 5, 6})); std::cout << "f456= " << f456 << " | f45= " << f45 << " - f56= " << f56 << " - f46= " << f46 << std::endl; - BOOST_CHECK(are_almost_the_same(f456, std::max(f45, std::max(f56,f46)))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f456, std::max(f45, std::max(f56,f46))); const int DIMENSION_3 = 3; Simplex_tree st3; @@ -324,7 +388,7 @@ BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) { Simplex_tree::Filtration_value f0123 = st3.filtration(st3.find({0, 1, 2, 3})); std::cout << "f0123= " << f0123 << " | f012= " << f012 << " - f123= " << f123 << " - f013= " << f013 << " - f023= " << f023 << std::endl; - BOOST_CHECK(are_almost_the_same(f0123, std::max(f012, std::max(f123, std::max(f013, f023))))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f0123, std::max(f012, std::max(f123, std::max(f013, f023)))); } diff --git a/src/Rips_complex/utilities/CMakeLists.txt b/src/Rips_complex/utilities/CMakeLists.txt new file mode 100644 index 00000000..4b565628 --- /dev/null +++ b/src/Rips_complex/utilities/CMakeLists.txt @@ -0,0 +1,34 @@ +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}) + +add_executable(rips_correlation_matrix_persistence rips_correlation_matrix_persistence.cpp) +target_link_libraries(rips_correlation_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(sparse_rips_persistence sparse_rips_persistence.cpp) +target_link_libraries(sparse_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}) + target_link_libraries(rips_correlation_matrix_persistence ${TBB_LIBRARIES}) + target_link_libraries(sparse_rips_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence> + "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") +add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +add_test(NAME Rips_complex_utility_from_rips_correlation_matrix COMMAND $<TARGET_FILE:rips_correlation_matrix_persistence> + "${CMAKE_SOURCE_DIR}/data/correlation_matrix/lower_triangular_correlation_matrix.csv" "-c" "0.3" "-d" "3" "-p" "3" "-m" "0") +add_test(NAME Sparse_rips_complex_utility_on_tore_3D COMMAND $<TARGET_FILE:sparse_rips_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-e" "0.5" "-m" "0.2" "-d" "3" "-p" "2") + +install(TARGETS rips_distance_matrix_persistence DESTINATION bin) +install(TARGETS rips_persistence DESTINATION bin) +install(TARGETS rips_correlation_matrix_persistence DESTINATION bin) +install(TARGETS sparse_rips_persistence DESTINATION bin) diff --git a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp new file mode 100644 index 00000000..287e8915 --- /dev/null +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -0,0 +1,171 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Pawel Dlotko, Vincent Rouvreau + * + * Copyright (C) 2016 Inria + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/reader_utils.h> +#include <gudhi/writing_persistence_to_file.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <algorithm> // for sort + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; +using Correlation_matrix = std::vector<std::vector<Filtration_value>>; +using intervals_common = Gudhi::Persistence_interval_common<double, int>; + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& correlation_min, 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 correlation_min; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, correlation_min, dim_max, p, min_persistence); + + Correlation_matrix correlations = + Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file); + + Filtration_value threshold = 0; + + // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + correlations[i][j] = 1 - correlations[i][j]; + // Here we make sure that the values of corelations lie between -1 and 1. + // If not, we throw an exception. + if ((correlations[i][j] < -1) || (correlations[i][j] > 1)) { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate. \n"; + throw "The input matrix is not a correlation matrix. The program will now terminate. \n"; + } + if (correlations[i][j] > threshold) threshold = correlations[i][j]; + } + } + + Rips_complex rips_complex_from_file(correlations, 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); + // compute persistence + pcoh.compute_persistent_cohomology(min_persistence); + + // invert the persistence diagram. The reason for this procedure is the following: + // The input to the program is a corelation matrix M. When processing it, it is + // turned into 1-M and the obtained persistence intervals are in '1-M' units. + // Below we reverse every (birth,death) pair into (1-birth, 1-death) pair + // so that the input and the output to the program is expressed in the same + // units. + auto pairs = pcoh.get_persistent_pairs(); + std::vector<intervals_common> processed_persistence_intervals; + processed_persistence_intervals.reserve(pairs.size()); + for (auto pair : pairs) { + double birth = 1 - simplex_tree.filtration(get<0>(pair)); + double death = 1 - simplex_tree.filtration(get<1>(pair)); + unsigned dimension = (unsigned)simplex_tree.dimension(get<0>(pair)); + int field = get<2>(pair); + processed_persistence_intervals.push_back(intervals_common(birth, death, dimension, field)); + } + + // sort the processed intervals: + std::sort(processed_persistence_intervals.begin(), processed_persistence_intervals.end()); + + // and write them to a file + if (filediag.empty()) { + write_persistence_intervals_to_stream(processed_persistence_intervals); + } else { + std::ofstream out(filediag); + write_persistence_intervals_to_stream(processed_persistence_intervals, out); + } + return 0; +} + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()( + "input-file", po::value<std::string>(&csv_matrix_file), + "Name of file containing a corelation matrix. Can be square or lower triangular matrix. Separator is ';'."); + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "min-edge-corelation,c", po::value<Filtration_value>(&correlation_min)->default_value(0), + "Minimal corelation of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a corelation 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; + exit(-1); + } +} diff --git a/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp new file mode 100644 index 00000000..c73152cf --- /dev/null +++ b/src/Rips_complex/utilities/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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/reader_utils.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; +using Distance_matrix = std::vector<std::vector<Filtration_value>>; + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); + + Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file); + Rips_complex rips_complex_from_file(distances, threshold); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + return 0; +} + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()( + "input-file", po::value<std::string>(&csv_matrix_file), + "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'."); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of distance matrix.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Rips_complex/utilities/rips_persistence.cpp b/src/Rips_complex/utilities/rips_persistence.cpp new file mode 100644 index 00000000..9410b9c2 --- /dev/null +++ b/src/Rips_complex/utilities/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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree 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<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md new file mode 100644 index 00000000..6df49310 --- /dev/null +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -0,0 +1,109 @@ +--- +layout: page +title: "Rips complex" +meta_title: "Rips complex" +teaser: "" +permalink: /ripscomplex/ +--- +{::comment} +Leave the lines above as it is required by the web site generator 'Jekyll' +{:/comment} + + +## 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] <OFF input file>` + +**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_distance_matrix_persistence [options] <CSV input file>` + +where +`<CSV input file>` is the path to the file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. + +**Example** + +`rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + + +## rips_correlation_matrix_persistence ## + +Same as `rips_distance_matrix_persistence` but taking a correlation matrix as input. + +**Usage** + +`rips_correlation_matrix_persistence [options] <CSV input file>` + +where +`<CSV input file>` is the path to the file containing a correlation matrix. Can be square or lower triangular matrix. Separator is ';'. +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + +**Example** + +`rips_correlation_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + +**Warning** + +As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical tool will not work +properly, this is a known issue. + + +## sparse_rips_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* +of a sparse (1+epsilon)-approximation of the Rips complex defined on a set of input Euclidean 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 (`p` must be a prime number). + +**Usage** + +`sparse_rips_persistence [options] <OFF input file>` + +**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. +* `-e [ --approximation ]` (default = .5) Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex. +* `-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. + +**Example with Z/2Z coefficients** + +`sparse_rips_persistence ../../data/points/tore3D_1307.off -e .5 -m .2 -d 3 -p 2` diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp new file mode 100644 index 00000000..6d4d86fd --- /dev/null +++ b/src/Rips_complex/utilities/sparse_rips_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): Marc Glisse, Clément Maria + * + * 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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Sparse_rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; + +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon, + int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string off_file_points; + std::string filediag; + double epsilon; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, epsilon, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Sparse_rips sparse_rips(off_reader.get_point_cloud(), Gudhi::Euclidean_distance(), epsilon); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + sparse_rips.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, double& epsilon, + int& dim_max, int& p, Filtration_value& min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "approximation,e", po::value<double>(&epsilon)->default_value(.5), + "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a sparse (1+epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} |