diff options
author | Hind-M <hind.montassif@gmail.com> | 2022-05-24 16:55:02 +0200 |
---|---|---|
committer | Hind-M <hind.montassif@gmail.com> | 2022-05-24 16:55:02 +0200 |
commit | 050081c554416a90f2b4aee359f15c020af66303 (patch) | |
tree | 2b27233f55e4568ad010526e73be9b3b67246297 /src | |
parent | a6a68c11455a554619d8a5b5d2f92c1ddbf45e99 (diff) | |
parent | e415d500fe8739414803b8cc127484562079d27e (diff) |
Merge remote-tracking branch 'upstream/master' into cech_optimization
Diffstat (limited to 'src')
50 files changed, 1425 insertions, 996 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index 5ab23720..7f14f423 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -164,11 +164,11 @@ Table of Contents * <b>Requires:</b> \ref eigen ≥ 3.1.0 and \ref cgal ≥ 5.1.0. * * A weighted version for Alpha complex is available (cf. Alpha_complex). It is like a usual Alpha complex, but based - * on a <a href="https://doc.cgal.org/latest/Triangulation/index.html#title20">CGAL regular triangulation</a> instead + * on a <a href="https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT">CGAL regular triangulation</a> instead * of Delaunay. * * This example builds the CGAL weighted alpha shapes from a small molecule, and initializes the alpha complex with - * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13">CGAL 3d + * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes">CGAL 3d * weighted alpha shapes</a>. * * Then, it is asked to display information about the alpha complex. @@ -212,7 +212,7 @@ Table of Contents * Gudhi::alpha_complex::complexity::EXACT. * * This example builds the CGAL 3d weighted alpha shapes from a small molecule, and initializes the alpha complex with - * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13">CGAL 3d + * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes">CGAL 3d * weighted alpha shapes</a>. * * Then, it is asked to display information about the alpha complex. diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index ccc3d852..562ef139 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -12,7 +12,6 @@ #ifndef ALPHA_COMPLEX_3D_H_ #define ALPHA_COMPLEX_3D_H_ -#include <boost/version.hpp> #include <boost/variant.hpp> #include <boost/range/size.hpp> #include <boost/range/combine.hpp> @@ -58,8 +57,6 @@ namespace Gudhi { namespace alpha_complex { -thread_local double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001; - // Value_from_iterator returns the filtration value from an iterator on alpha shapes values // // FAST SAFE EXACT @@ -101,7 +98,7 @@ struct Value_from_iterator<complexity::EXACT> { * \tparam Periodic Boolean used to set/unset the periodic version of Alpha_complex_3d. Default value is false. * * For the weighted version, weights values are explained on CGAL - * <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0">Alpha shapes 3d</a> and + * <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#Alpha_shapes_3Definitions">Alpha shapes 3d</a> and * <a href="https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation">Regular * triangulation</a> documentation. * diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index 0d3c6027..1e3b8fab 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -64,7 +64,7 @@ N.B.: * Weights values are explained on CGAL
[dD Triangulations](https://doc.cgal.org/latest/Triangulation/index.html)
and
-[Regular triangulation](https://doc.cgal.org/latest/Triangulation/index.html#title20) documentation.
+[Regular triangulation](https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT) documentation.
## alpha_complex_3d_persistence ##
@@ -131,6 +131,6 @@ N.B.: * `alpha_complex_3d_persistence` only accepts OFF files in dimension 3.
* Filtration values are alpha square values.
* Weights values are explained on CGAL
-[Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0)
+[Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#Alpha_shapes_3Definitions)
and
[Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation.
diff --git a/src/Collapse/doc/intro_edge_collapse.h b/src/Collapse/doc/intro_edge_collapse.h index fde39707..12e909c8 100644 --- a/src/Collapse/doc/intro_edge_collapse.h +++ b/src/Collapse/doc/intro_edge_collapse.h @@ -17,68 +17,48 @@ namespace collapse { /** \defgroup edge_collapse Edge collapse * - * \author Siddharth Pritam + * \author Siddharth Pritam and Marc Glisse * * @{ * - * This module implements edge collapse of a filtered flag complex, in particular it reduces a filtration of - * Vietoris-Rips complex from its graph to another smaller flag filtration with same persistence. - * Where a filtration is a sequence of simplicial (here Rips) complexes connected with inclusions. + * This module implements edge collapse of a filtered flag complex as described in \cite edgecollapsearxiv, in + * particular it reduces a filtration of Vietoris-Rips complex represented by a graph to a smaller flag filtration with + * the same persistent homology. * * \section edge_collapse_definition Edge collapse definition * * An edge \f$e\f$ in a simplicial complex \f$K\f$ is called a <b>dominated edge</b> if the link of \f$e\f$ in * \f$K\f$, \f$lk_K(e)\f$ is a simplicial cone, that is, there exists a vertex \f$v^{\prime} \notin e\f$ and a - * subcomplex \f$L\f$ in \f$K\f$, such that \f$lk_K(e) = v^{\prime}L\f$. We say that the vertex \f$v^{\prime}\f$ is - * {dominating} \f$e\f$ and \f$e\f$ is {dominated} by \f$v^{\prime}\f$. - * An <b> elementary egde collapse </b> is the removal of a dominated edge \f$e\f$ from \f$K\f$, - * which we denote with \f$K\f$ \f${\searrow\searrow}^1 \f$ \f$K\setminus e\f$. - * The symbol \f$\mathbf{K\setminus e}\f$ (deletion of \f$e\f$ from \f$K\f$) refers to the subcomplex of \f$K\f$ which - * has all simplices of \f$K\f$ except \f$e\f$ and the ones containing \f$e\f$. - * There is an <b>edge collapse</b> from a simplicial complex \f$K\f$ to its subcomplex \f$L\f$, - * if there exists a series of elementary edge collapses from \f$K\f$ to \f$L\f$, denoted as \f$K\f$ - * \f${\searrow\searrow}\f$ \f$L\f$. - * - * An edge collapse is a homotopy preserving operation, and it can be further expressed as sequence of the classical - * elementary simple collapse. - * A complex without any dominated edge is called a \f$1\f$- minimal complex and the core \f$K^1\f$ of simplicial - * complex is a minimal complex such that \f$K\f$ \f${\searrow\searrow}\f$ \f$K^1\f$. - * Computation of a core (not unique) involves computation of dominated edges and the dominated edges can be easily - * characterized as follows: + * subcomplex \f$L\f$ in \f$K\f$, such that \f$lk_K(e) = v^{\prime}L\f$. We say that the vertex \f$v^{\prime}\f$ + * \e dominates \f$e\f$ and \f$e\f$ is \e dominated by \f$v^{\prime}\f$. + * An <b> elementary edge collapse </b> is the removal of a dominated edge \f$e\f$ from \f$K\f$ (the cofaces of \f$e\f$ + * are implicitly removed as well). + * Domination is used as a simple sufficient condition that ensures that this removal is a homotopy preserving + * operation. + * + * The dominated edges can be easily characterized as follows: * - * -- For general simplicial complex: An edge \f$e \in K\f$ is dominated by another vertex \f$v^{\prime} \in K\f$, - * <i>if and only if</i> all the maximal simplices of \f$K\f$ that contain \f$e\f$ also contain \f$v^{\prime}\f$ + * -- For a general simplicial complex: an edge \f$e \in K\f$ is dominated by another vertex \f$v^{\prime} \in K\f$, + * if and only if all the maximal simplices of \f$K\f$ that contain \f$e\f$ also contain \f$v^{\prime}\f$. * - * -- For a flag complex: An edge \f$e \in K\f$ is dominated by another vertex \f$v^{\prime} \in K\f$, <i>if and only - * if</i> all the vertices in \f$K\f$ that has an edge with both vertices of \f$e\f$ also has an edge with - * \f$v^{\prime}\f$. + * -- For a flag complex: an edge \f$e \in K\f$ is dominated by another vertex \f$v^{\prime} \in K\f$, if and only + * if all the vertices in \f$K\f$ that have an edge with both vertices of \f$e\f$ also have an edge with + * \f$v^{\prime}\f$. Notice that this only depends on the graph. * - * The algorithm to compute the smaller induced filtration is described in Section 5 \cite edgecollapsesocg2020. - * Edge collapse can be successfully employed to reduce any given filtration of flag complexes to a smaller induced + * In the context of a filtration, an edge collapse may translate into an increase of the filtration value of an edge, + * or its removal if it already had the largest filtration value. + * The algorithm to compute the smaller induced filtration is described in \cite edgecollapsearxiv. + * Edge collapse can be successfully employed to reduce any input filtration of flag complexes to a smaller induced * filtration which preserves the persistent homology of the original filtration and is a flag complex as well. * - * The general idea is that we consider edges in the filtered graph and sort them according to their filtration value - * giving them a total order. - * Each edge gets a unique index denoted as \f$i\f$ in this order. To reduce the filtration, we move forward with - * increasing filtration value - * in the graph and check if the current edge \f$e_i\f$ is dominated in the current graph \f$G_i := \{e_1, .. e_i\} \f$ - * or not. - * If the edge \f$e_i\f$ is dominated we remove it from the filtration and move forward to the next edge \f$e_{i+1}\f$. - * If \f$e_i\f$ is non-dominated then we keep it in the reduced filtration and then go backward in the current graph - * \f$G_i\f$ to look for new non-dominated edges that was dominated before but might become non-dominated at this - * point. - * If an edge \f$e_j, j < i \f$ during the backward search is found to be non-dominated, we include \f$e_j\f$ in to the - * reduced filtration and we set its new filtration value to be \f$i\f$ that is the index of \f$e_i\f$. - * The precise mechanism for this reduction has been described in Section 5 \cite edgecollapsesocg2020. - * Here we implement this mechanism for a filtration of Rips complex. - * After perfoming the reduction the filtration reduces to a flag-filtration with the same persistence as the original - * filtration. - * + * The algorithm implemented here does not produce a minimal filtration. Taking its output and applying the algorithm a + * second time may further simplify the filtration. + * * \subsection edgecollapseexample Basic edge collapse * * This example calls `Gudhi::collapse::flag_complex_collapse_edges()` from a proximity graph represented as a list of * `Filtered_edge`. - * Then it collapses edges and displays a new list of `Filtered_edge` (with less edges) + * Then it collapses edges and displays a new list of `Filtered_edge` (with fewer edges) * that will preserve the persistence homology computation. * * \include edge_collapse_basic_example.cpp @@ -88,7 +68,7 @@ namespace collapse { * \code $> ./Edge_collapse_example_basic * \endcode * - * the program output is: + * the program output could be: * * \include edge_collapse_example_basic.txt */ diff --git a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h index 713c6608..c823901f 100644 --- a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h +++ b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h @@ -1,11 +1,12 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Siddharth Pritam + * Author(s): Siddharth Pritam, Marc Glisse * * Copyright (C) 2020 Inria * * Modification(s): * - 2020/03 Vincent Rouvreau: integration to the gudhi library + * - 2021 Marc Glisse: complete rewrite * - YYYY/MM Author: Description of the modification */ @@ -14,367 +15,319 @@ #include <gudhi/Debug_utils.h> -#include <boost/functional/hash.hpp> -#include <boost/iterator/iterator_facade.hpp> - -#include <Eigen/Sparse> -#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST +#include <boost/container/flat_map.hpp> +#include <boost/container/flat_set.hpp> #ifdef GUDHI_USE_TBB #include <tbb/parallel_sort.h> #endif -#include <iostream> -#include <utility> // for std::pair +#include <utility> #include <vector> -#include <unordered_map> -#include <unordered_set> -#include <set> -#include <tuple> // for std::tie -#include <algorithm> // for std::includes -#include <iterator> // for std::inserter -#include <type_traits> // for std::decay - -// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 -#if !EIGEN_VERSION_AT_LEAST(3,1,0) -# error Edge Collapse is only available for Eigen3 >= 3.1.0 -#endif +#include <tuple> +#include <algorithm> +#include <limits> namespace Gudhi { namespace collapse { /** \private - * - * \brief Flag complex sparse matrix data structure. * - * \details - * This class stores a <a target="_blank" href="https://en.wikipedia.org/wiki/Clique_complex">Flag complex</a> - * in an <a target="_blank" href="https://eigen.tuxfamily.org/dox/group__TutorialSparse.html">Eigen sparse matrix</a>. + * \brief Flag complex sparse matrix data structure. * - * \tparam Vertex type must be a signed integer type. It admits a total order <. - * \tparam Filtration type for the value of the filtration function. Must be comparable with <. + * \tparam Vertex type must be an integer type. + * \tparam Filtration type for the value of the filtration function. */ -template<typename Vertex, typename Filtration> -class Flag_complex_edge_collapser { - public: - /** \brief Re-define Vertex as Vertex_handle type to ease the interface with `Gudhi::Proximity_graph`. */ - using Vertex_handle = Vertex; - /** \brief Re-define Filtration as Filtration_value type to ease the interface with `Gudhi::Proximity_graph`. */ - using Filtration_value = Filtration; - - private: - // internal numbering of vertices and edges - using IVertex = std::size_t; - using Edge_index = std::size_t; - using IEdge = std::pair<IVertex, IVertex>; - - // The sparse matrix data type - // (Eigen::SparseMatrix<Edge_index, Eigen::RowMajor> has slow insertions) - using Sparse_vector = Eigen::SparseVector<Edge_index>; - using Sparse_row_matrix = std::vector<Sparse_vector>; - - // Range of neighbors of a vertex - template<bool closed> - struct Neighbours { - class iterator : public boost::iterator_facade<iterator, - IVertex, /* value_type */ - std::input_iterator_tag, // or boost::single_pass_traversal_tag - IVertex /* reference */ > - { - public: - iterator():ptr(nullptr){} - iterator(Neighbours const*p):ptr(p){find_valid();} - private: - friend class boost::iterator_core_access; - Neighbours const*ptr; - void increment(){ - ++ptr->it; - find_valid(); - } - void find_valid(){ - auto& it = ptr->it; - do { - if(!it) { ptr=nullptr; break; } - if(IVertex(it.index()) == ptr->u) { - if(closed) break; - else continue; - } - Edge_index e = it.value(); - if(e <= ptr->ec->current_backward || ptr->ec->critical_edge_indicator_[e]) break; - } while(++it, true); - } - bool equal(iterator const& other) const { return ptr == other.ptr; } - IVertex dereference() const { return ptr->it.index(); } - }; - typedef iterator const_iterator; - mutable typename Sparse_vector::InnerIterator it; - Flag_complex_edge_collapser const*ec; - IVertex u; - iterator begin() const { return this; } - iterator end() const { return {}; } - explicit Neighbours(Flag_complex_edge_collapser const*p,IVertex u):it(p->sparse_row_adjacency_matrix_[u]),ec(p),u(u){} - }; - - // A range of row indices - using IVertex_vector = std::vector<IVertex>; - - public: - /** \brief Filtered_edge is a type to store an edge with its filtration value. */ - using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>; - - private: - // Map from row index to its vertex handle - std::vector<Vertex_handle> row_to_vertex_; - - // Index of the current edge in the backwards walk. Edges <= current_backward are part of the temporary graph, - // while edges > current_backward are removed unless critical_edge_indicator_. - Edge_index current_backward = -1; - - // Map from IEdge to its index - std::unordered_map<IEdge, Edge_index, boost::hash<IEdge>> iedge_to_index_map_; - - // Boolean vector to indicate if the edge is critical. - std::vector<bool> critical_edge_indicator_; - - // Map from vertex handle to its row index - std::unordered_map<Vertex_handle, IVertex> vertex_to_row_; - - // Stores the Sparse matrix of Filtration values representing the original graph. - // The matrix rows and columns are indexed by IVertex. - Sparse_row_matrix sparse_row_adjacency_matrix_; - - // The input, a vector of filtered edges. - std::vector<Filtered_edge> f_edge_vector_; - - // Edge is the actual edge (u,v), with Vertex_handle u and v, not IVertex. - bool edge_is_dominated(Vertex_handle u, Vertex_handle v) const - { - const IVertex rw_u = vertex_to_row_.at(u); - const IVertex rw_v = vertex_to_row_.at(v); -#ifdef DEBUG_TRACES - std::cout << "The edge {" << u << ", " << v << "} is going for domination check." << std::endl; -#endif // DEBUG_TRACES - auto common_neighbours = open_common_neighbours_row_index(rw_u, rw_v); -#ifdef DEBUG_TRACES - std::cout << "And its common neighbours are." << std::endl; - for (auto neighbour : common_neighbours) { - std::cout << row_to_vertex_[neighbour] << ", " ; - } - std::cout<< std::endl; -#endif // DEBUG_TRACES - if (common_neighbours.size() == 1) - return true; - else - for (auto rw_c : common_neighbours) { - auto neighbours_c = neighbours_row_index<true>(rw_c); - // If neighbours_c contains the common neighbours. - if (std::includes(neighbours_c.begin(), neighbours_c.end(), - common_neighbours.begin(), common_neighbours.end())) - return true; - } - return false; +template<typename Vertex, typename Filtration_value> +struct Flag_complex_edge_collapser { + using Filtered_edge = std::tuple<Vertex, Vertex, Filtration_value>; + typedef std::pair<Vertex,Vertex> Edge; + struct Cmpi { template<class T, class U> bool operator()(T const&a, U const&b)const{return b<a; } }; + typedef boost::container::flat_map<Vertex, Filtration_value> Ngb_list; + typedef std::vector<Ngb_list> Neighbors; + Neighbors neighbors; // closed neighborhood + std::size_t num_vertices; + std::vector<std::tuple<Vertex, Vertex, Filtration_value>> res; + +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + // Minimal matrix interface + // Using this matrix generally helps performance, but the memory use may be excessive for a very sparse graph + // (and in extreme cases the constant initialization of the matrix may start to dominate the runnning time). + // Are there cases where the matrix is too big but a hash table would help? + std::vector<Filtration_value> neighbors_data; + void init_neighbors_dense(){ + neighbors_data.clear(); + neighbors_data.resize(num_vertices*num_vertices, std::numeric_limits<Filtration_value>::infinity()); } + Filtration_value& neighbors_dense(Vertex i, Vertex j){return neighbors_data[num_vertices*j+i];} +#endif - // Returns the edges connecting u and v (extremities of crit) to their common neighbors (not themselves) - std::set<Edge_index> three_clique_indices(Edge_index crit) { - std::set<Edge_index> edge_indices; - - Vertex_handle u = std::get<0>(f_edge_vector_[crit]); - Vertex_handle v = std::get<1>(f_edge_vector_[crit]); - -#ifdef DEBUG_TRACES - std::cout << "The current critical edge to re-check criticality with filt value is : f {" << u << "," << v - << "} = " << std::get<2>(f_edge_vector_[crit]) << std::endl; -#endif // DEBUG_TRACES - auto rw_u = vertex_to_row_[u]; - auto rw_v = vertex_to_row_[v]; - - IVertex_vector common_neighbours = open_common_neighbours_row_index(rw_u, rw_v); + // This does not touch the events list, only the adjacency matrix(es) + void delay_neighbor(Vertex u, Vertex v, Filtration_value f) { + neighbors[u][v]=f; + neighbors[v][u]=f; +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + neighbors_dense(u,v)=f; + neighbors_dense(v,u)=f; +#endif + } + void remove_neighbor(Vertex u, Vertex v) { + neighbors[u].erase(v); + neighbors[v].erase(u); +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + neighbors_dense(u,v)=std::numeric_limits<Filtration_value>::infinity(); + neighbors_dense(v,u)=std::numeric_limits<Filtration_value>::infinity(); +#endif + } - for (auto rw_c : common_neighbours) { - IEdge e_with_new_nbhr_v = std::minmax(rw_u, rw_c); - IEdge e_with_new_nbhr_u = std::minmax(rw_v, rw_c); - edge_indices.emplace(iedge_to_index_map_[e_with_new_nbhr_v]); - edge_indices.emplace(iedge_to_index_map_[e_with_new_nbhr_u]); + template<class FilteredEdgeRange> + void read_edges(FilteredEdgeRange const&r){ + neighbors.resize(num_vertices); +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + init_neighbors_dense(); +#endif + // Use the raw sequence to avoid maintaining the order + std::vector<typename Ngb_list::sequence_type> neighbors_seq(num_vertices); + for(auto&&e : r){ + using std::get; + Vertex u = get<0>(e); + Vertex v = get<1>(e); + Filtration_value f = get<2>(e); + neighbors_seq[u].emplace_back(v, f); + neighbors_seq[v].emplace_back(u, f); +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + neighbors_dense(u,v)=f; + neighbors_dense(v,u)=f; +#endif + } + for(std::size_t i=0;i<neighbors_seq.size();++i){ + neighbors_seq[i].emplace_back(i, -std::numeric_limits<Filtration_value>::infinity()); + neighbors[i].adopt_sequence(std::move(neighbors_seq[i])); // calls sort +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + neighbors_dense(i,i)=-std::numeric_limits<Filtration_value>::infinity(); +#endif } - return edge_indices; } - // Detect and set all edges that are becoming critical - template<typename FilteredEdgeOutput> - void set_edge_critical(Edge_index indx, Filtration_value filt, FilteredEdgeOutput filtered_edge_output) { -#ifdef DEBUG_TRACES - std::cout << "The curent index with filtration value " << indx << ", " << filt << " is primary critical" << - std::endl; -#endif // DEBUG_TRACES - std::set<Edge_index> effected_indices = three_clique_indices(indx); - // Cannot use boost::adaptors::reverse in such dynamic cases apparently - for (auto it = effected_indices.rbegin(); it != effected_indices.rend(); ++it) { - current_backward = *it; - Vertex_handle u = std::get<0>(f_edge_vector_[current_backward]); - Vertex_handle v = std::get<1>(f_edge_vector_[current_backward]); - // If current_backward is not critical so it should be processed, otherwise it stays in the graph - if (!critical_edge_indicator_[current_backward]) { - if (!edge_is_dominated(u, v)) { -#ifdef DEBUG_TRACES - std::cout << "The curent index became critical " << current_backward << std::endl; -#endif // DEBUG_TRACES - critical_edge_indicator_[current_backward] = true; - filtered_edge_output(u, v, filt); - std::set<Edge_index> inner_effected_indcs = three_clique_indices(current_backward); - for (auto inr_idx : inner_effected_indcs) { - if(inr_idx < current_backward) // && !critical_edge_indicator_[inr_idx] - effected_indices.emplace(inr_idx); - } -#ifdef DEBUG_TRACES - std::cout << "The following edge is critical with filt value: {" << u << "," << v << "}; " - << filt << std::endl; -#endif // DEBUG_TRACES - } + // Open neighborhood + // At some point it helped gcc to add __attribute__((noinline)) here, otherwise we had +50% on the running time + // on one example. It looks ok now, or I forgot which example that was. + void common_neighbors(boost::container::flat_set<Vertex>& e_ngb, + std::vector<std::pair<Filtration_value, Vertex>>& e_ngb_later, + Vertex u, Vertex v, Filtration_value f_event){ + // Using neighbors_dense here seems to hurt, even if we loop on the smaller of nu and nv. + Ngb_list const&nu = neighbors[u]; + Ngb_list const&nv = neighbors[v]; + auto ui = nu.begin(); + auto ue = nu.end(); + auto vi = nv.begin(); + auto ve = nv.end(); + assert(ui != ue && vi != ve); + while(ui != ue && vi != ve){ + Vertex w = ui->first; + if(w < vi->first) { ++ui; continue; } + if(w > vi->first) { ++vi; continue; } + // nu and nv are closed, so we need to exclude e here. + if(w != u && w != v) { + Filtration_value f = std::max(ui->second, vi->second); + if(f > f_event) + e_ngb_later.emplace_back(f, w); + else + e_ngb.insert(e_ngb.end(), w); } + ++ui; ++vi; } - // Clear the implicit "removed from graph" data structure - current_backward = -1; } - // Returns list of neighbors of a particular vertex. - template<bool closed> - auto neighbours_row_index(IVertex rw_u) const - { - return Neighbours<closed>(this, rw_u); + // Test if the neighborhood of e is included in the closed neighborhood of c + template<class Ngb> + bool is_dominated_by(Ngb const& e_ngb, Vertex c, Filtration_value f){ + // The best strategy probably depends on the distribution, how sparse / dense the adjacency matrix is, + // how (un)balanced the sizes of e_ngb and nc are. + // Some efficient operations on sets work best with bitsets, although the need for a map complicates things. +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + for(auto v : e_ngb) { + // if(v==c)continue; + if(neighbors_dense(v,c) > f) return false; + } + return true; +#else + auto&&nc = neighbors[c]; + // if few neighbors, use dichotomy? Seems slower. + // I tried storing a copy of neighbors as a vector<absl::flat_hash_map> and using it for nc, but it was + // a bit slower here. It did help with neighbors[dominator].find(w) in the main function though, + // sometimes enough, sometimes not. + auto ci = nc.begin(); + auto ce = nc.end(); + auto eni = e_ngb.begin(); + auto ene = e_ngb.end(); + assert(eni != ene); + assert(ci != ce); + // if(*eni == c && ++eni == ene) return true; + for(;;){ + Vertex ve = *eni; + Vertex vc = ci->first; + while(ve > vc) { + // try a gallop strategy (exponential search)? Seems slower + if(++ci == ce) return false; + vc = ci->first; + } + if(ve < vc) return false; + // ve == vc + if(ci->second > f) return false; + if(++eni == ene)return true; + // If we stored an open neighborhood of c (excluding c), we would need to test for c here and before the loop + // if(*eni == c && ++eni == ene)return true; + if(++ci == ce) return false; + } +#endif } - // Returns the list of open neighbours of the edge :{u,v}. - IVertex_vector open_common_neighbours_row_index(IVertex rw_u, IVertex rw_v) const - { - auto non_zero_indices_u = neighbours_row_index<false>(rw_u); - auto non_zero_indices_v = neighbours_row_index<false>(rw_v); - IVertex_vector common; - std::set_intersection(non_zero_indices_u.begin(), non_zero_indices_u.end(), non_zero_indices_v.begin(), - non_zero_indices_v.end(), std::back_inserter(common)); - - return common; - } + template<class FilteredEdgeRange, class Delay> + void process_edges(FilteredEdgeRange const& edges, Delay&& delay) { + { + Vertex maxi = 0, maxj = 0; + for(auto& fe : edges) { + Vertex i = std::get<0>(fe); + Vertex j = std::get<1>(fe); + if (i > maxi) maxi = i; + if (j > maxj) maxj = j; + } + num_vertices = std::max(maxi, maxj) + 1; + } - // Insert a vertex in the data structure - IVertex insert_vertex(Vertex_handle vertex) { - auto n = row_to_vertex_.size(); - auto result = vertex_to_row_.emplace(vertex, n); - // If it was not already inserted - Value won't be updated by emplace if it is already present - if (result.second) { - // Expand the matrix. The size of rows is irrelevant. - sparse_row_adjacency_matrix_.emplace_back((std::numeric_limits<Eigen::Index>::max)()); - // Initializing the diagonal element of the adjency matrix corresponding to rw_b. - sparse_row_adjacency_matrix_[n].insert(n) = -1; // not an edge - // Must be done after reading its size() - row_to_vertex_.push_back(vertex); + read_edges(edges); + + boost::container::flat_set<Vertex> e_ngb; + e_ngb.reserve(num_vertices); + std::vector<std::pair<Filtration_value, Vertex>> e_ngb_later; + for(auto&e:edges) { + { + Vertex u = std::get<0>(e); + Vertex v = std::get<1>(e); + Filtration_value input_time = std::get<2>(e); + auto time = delay(input_time); + auto start_time = time; + e_ngb.clear(); + e_ngb_later.clear(); + common_neighbors(e_ngb, e_ngb_later, u, v, time); + // If we identify a good candidate (the first common neighbor) for being a dominator of e until infinity, + // we could check that a bit more cheaply. It does not seem to help though. + auto cmp1=[](auto const&a, auto const&b){return a.first > b.first;}; + auto e_ngb_later_begin=e_ngb_later.begin(); + auto e_ngb_later_end=e_ngb_later.end(); + bool heapified = false; + + bool dead = false; + while(true) { + Vertex dominator = -1; + // special case for size 1 + // if(e_ngb.size()==1){dominator=*e_ngb.begin();}else + // It is tempting to test the dominators in increasing order of filtration value, which is likely to reduce + // the number of calls to is_dominated_by before finding a dominator, but sorting, even partially / lazily, + // is very expensive. + for(auto c : e_ngb){ + if(is_dominated_by(e_ngb, c, time)){ + dominator = c; + break; + } + } + if(dominator==-1) break; + // Push as long as dominator remains a dominator. + // Iterate on times where at least one neighbor appears. + for (bool still_dominated = true; still_dominated; ) { + if(e_ngb_later_begin == e_ngb_later_end) { + dead = true; goto end_move; + } + if(!heapified) { + // Eagerly sorting can be slow + std::make_heap(e_ngb_later_begin, e_ngb_later_end, cmp1); + heapified=true; + } + time = e_ngb_later_begin->first; // first place it may become critical + // Update the neighborhood for this new time, while checking if any of the new neighbors break domination. + while (e_ngb_later_begin != e_ngb_later_end && e_ngb_later_begin->first <= time) { + Vertex w = e_ngb_later_begin->second; +#ifdef GUDHI_COLLAPSE_USE_DENSE_ARRAY + if (neighbors_dense(dominator,w) > e_ngb_later_begin->first) + still_dominated = false; +#else + auto& ngb_dom = neighbors[dominator]; + auto wit = ngb_dom.find(w); // neighborhood may be open or closed, it does not matter + if (wit == ngb_dom.end() || wit->second > e_ngb_later_begin->first) + still_dominated = false; +#endif + e_ngb.insert(w); + std::pop_heap(e_ngb_later_begin, e_ngb_later_end--, cmp1); + } + } // this doesn't seem to help that much... + } +end_move: + if(dead) { + remove_neighbor(u, v); + } else if(start_time != time) { + delay_neighbor(u, v, time); + res.emplace_back(u, v, time); + } else { + res.emplace_back(u, v, input_time); + } + } } - return result.first->second; } - // Insert an edge in the data structure - // @exception std::invalid_argument In debug mode, if u == v - IEdge insert_new_edge(Vertex_handle u, Vertex_handle v, Edge_index idx) - { - GUDHI_CHECK((u != v), std::invalid_argument("Flag_complex_edge_collapser::insert_new_edge with u == v")); - // The edge must not be added before, it should be a new edge. - IVertex rw_u = insert_vertex(u); - IVertex rw_v = insert_vertex(v); -#ifdef DEBUG_TRACES - std::cout << "Inserting the edge " << u <<", " << v << std::endl; -#endif // DEBUG_TRACES - sparse_row_adjacency_matrix_[rw_u].insert(rw_v) = idx; - sparse_row_adjacency_matrix_[rw_v].insert(rw_u) = idx; - return std::minmax(rw_u, rw_v); + std::vector<Filtered_edge> output() { + return std::move(res); } - public: - /** \brief Flag_complex_edge_collapser constructor from a range of filtered edges. - * - * @param[in] edges Range of Filtered edges range.There is no need the range to be sorted, as it will be performed in - * `Flag_complex_edge_collapser::process_edges`. - * - * \tparam FilteredEdgeRange must be a range for which std::begin and std::end return iterators on a - * `Flag_complex_edge_collapser::Filtered_edge`. - */ - template<typename FilteredEdgeRange> - Flag_complex_edge_collapser(const FilteredEdgeRange& edges) - : f_edge_vector_(std::begin(edges), std::end(edges)) { } +}; - /** \brief Performs edge collapse in a increasing sequence of the filtration value. - * - * \tparam filtered_edge_output is a functor that is called on the output edges, in non-decreasing order of - * filtration, as filtered_edge_output(u, v, f) where u and v are Vertex_handle representing the extremities of the - * edge, and f is its new Filtration_value. - */ - template<typename FilteredEdgeOutput> - void process_edges(FilteredEdgeOutput filtered_edge_output) { - // Sort edges - auto sort_by_filtration = [](const Filtered_edge& edge_a, const Filtered_edge& edge_b) -> bool - { - return (std::get<2>(edge_a) < std::get<2>(edge_b)); - }; +template<class R> R to_range(R&& r) { return std::move(r); } +template<class R, class T> R to_range(T&& t) { R r; r.insert(r.end(), t.begin(), t.end()); return r; } +template<class FilteredEdgeRange, class Delay> +auto flag_complex_collapse_edges(FilteredEdgeRange&& edges, Delay&&delay) { + // Would it help to label the points according to some spatial sorting? + auto first_edge_itr = std::begin(edges); + using Vertex = std::decay_t<decltype(std::get<0>(*first_edge_itr))>; + using Filtration_value = std::decay_t<decltype(std::get<2>(*first_edge_itr))>; + using Edge_collapser = Flag_complex_edge_collapser<Vertex, Filtration_value>; + if (first_edge_itr != std::end(edges)) { + auto edges2 = to_range<std::vector<typename Edge_collapser::Filtered_edge>>(std::forward<FilteredEdgeRange>(edges)); #ifdef GUDHI_USE_TBB - tbb::parallel_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration); + // I think this sorting is always negligible compared to the collapse, but parallelizing it shouldn't hurt. + tbb::parallel_sort(edges2.begin(), edges2.end(), + [](auto const&a, auto const&b){return std::get<2>(a)>std::get<2>(b);}); #else - std::sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration); + std::sort(edges2.begin(), edges2.end(), [](auto const&a, auto const&b){return std::get<2>(a)>std::get<2>(b);}); #endif - - for (Edge_index endIdx = 0; endIdx < f_edge_vector_.size(); endIdx++) { - Filtered_edge fec = f_edge_vector_[endIdx]; - Vertex_handle u = std::get<0>(fec); - Vertex_handle v = std::get<1>(fec); - Filtration_value filt = std::get<2>(fec); - - // Inserts the edge in the sparse matrix to update the graph (G_i) - IEdge ie = insert_new_edge(u, v, endIdx); - - iedge_to_index_map_.emplace(ie, endIdx); - critical_edge_indicator_.push_back(false); - - if (!edge_is_dominated(u, v)) { - critical_edge_indicator_[endIdx] = true; - filtered_edge_output(u, v, filt); - if (endIdx > 1) - set_edge_critical(endIdx, filt, filtered_edge_output); - } - } + Edge_collapser edge_collapser; + edge_collapser.process_edges(edges2, std::forward<Delay>(delay)); + return edge_collapser.output(); } - -}; + return std::vector<typename Edge_collapser::Filtered_edge>(); +} /** \brief Implicitly constructs a flag complex from edges as an input, collapses edges while preserving the persistent - * homology and returns the remaining edges as a range. + * homology and returns the remaining edges as a range. The filtration value of vertices is irrelevant to this function. * - * \param[in] edges Range of Filtered edges.There is no need the range to be sorted, as it will be performed. + * \param[in] edges Range of Filtered edges. There is no need for the range to be sorted, as it will be done internally. * - * \tparam FilteredEdgeRange furnishes `std::begin` and `std::end` methods and returns an iterator on a - * FilteredEdge of type `std::tuple<Vertex_handle, Vertex_handle, Filtration_value>` where `Vertex_handle` is the type - * of a vertex index and `Filtration_value` is the type of an edge filtration value. + * \tparam FilteredEdgeRange Range of `std::tuple<Vertex_handle, Vertex_handle, Filtration_value>` + * where `Vertex_handle` is the type of a vertex index. * * \return Remaining edges after collapse as a range of * `std::tuple<Vertex_handle, Vertex_handle, Filtration_value>`. - * + * * \ingroup edge_collapse - * + * + * \note + * Advanced: Defining the macro GUDHI_COLLAPSE_USE_DENSE_ARRAY tells gudhi to allocate a square table of size the + * maximum vertex index. This usually speeds up the computation for dense graphs. However, for sparse graphs, the memory + * use may be problematic and initializing this large table may be slow. */ template<class FilteredEdgeRange> auto flag_complex_collapse_edges(const FilteredEdgeRange& edges) { - auto first_edge_itr = std::begin(edges); - using Vertex_handle = std::decay_t<decltype(std::get<0>(*first_edge_itr))>; - using Filtration_value = std::decay_t<decltype(std::get<2>(*first_edge_itr))>; - using Edge_collapser = Flag_complex_edge_collapser<Vertex_handle, Filtration_value>; - std::vector<typename Edge_collapser::Filtered_edge> remaining_edges; - if (first_edge_itr != std::end(edges)) { - Edge_collapser edge_collapser(edges); - edge_collapser.process_edges( - [&remaining_edges](Vertex_handle u, Vertex_handle v, Filtration_value filtration) { - // insert the edge - remaining_edges.emplace_back(u, v, filtration); - }); - } - return remaining_edges; + return flag_complex_collapse_edges(edges, [](auto const&d){return d;}); } } // namespace collapse diff --git a/src/Collapse/test/collapse_unit_test.cpp b/src/Collapse/test/collapse_unit_test.cpp index b8876246..f41dbedd 100644 --- a/src/Collapse/test/collapse_unit_test.cpp +++ b/src/Collapse/test/collapse_unit_test.cpp @@ -98,8 +98,8 @@ BOOST_AUTO_TEST_CASE(collapse) { // o---o // 0 3 edges.emplace_back(0, 2, 2.); - edges.emplace_back(1, 3, 2.); - trace_and_check_collapse(edges, {{1, 3, 2.}}); + edges.emplace_back(1, 3, 2.1); + trace_and_check_collapse(edges, {{1, 3, 2.1}}); // 1 2 4 // o---o---o @@ -121,8 +121,8 @@ BOOST_AUTO_TEST_CASE(collapse) { // o---o---o // 0 3 5 edges.emplace_back(2, 5, 4.); - edges.emplace_back(4, 3, 4.); - trace_and_check_collapse(edges, {{1, 3, 2.}, {4, 3, 4.}}); + edges.emplace_back(4, 3, 4.1); + trace_and_check_collapse(edges, {{1, 3, 2.}, {4, 3, 4.1}}); // 1 2 4 // o---o---o @@ -132,8 +132,8 @@ BOOST_AUTO_TEST_CASE(collapse) { // o---o---o // 0 3 5 edges.emplace_back(1, 5, 5.); - edges.emplace_back(0, 4, 5.); - trace_and_check_collapse(edges, {{1, 3, 2.}, {4, 3, 4.}, {0, 4, 5.}}); + edges.emplace_back(0, 4, 5.1); + trace_and_check_collapse(edges, {{1, 3, 2.}, {4, 3, 4.}, {0, 4, 5.1}}); } BOOST_AUTO_TEST_CASE(collapse_from_array) { @@ -150,8 +150,8 @@ BOOST_AUTO_TEST_CASE(collapse_from_array) { {2, 3, 1.}, {3, 0, 1.}, {0, 2, 2.}, - {1, 3, 2.}}}; - trace_and_check_collapse(f_edge_array, {{1, 3, 2.}}); + {1, 3, 2.1}}}; + trace_and_check_collapse(f_edge_array, {{1, 3, 2.1}}); } BOOST_AUTO_TEST_CASE(collapse_from_proximity_graph) { diff --git a/src/Contraction/include/gudhi/Edge_contraction.h b/src/Contraction/include/gudhi/Edge_contraction.h index 6c0f4c78..58d627c2 100644 --- a/src/Contraction/include/gudhi/Edge_contraction.h +++ b/src/Contraction/include/gudhi/Edge_contraction.h @@ -26,6 +26,7 @@ namespace contraction { /** \defgroup contr Edge contraction +@{ \author David Salinas @@ -195,7 +196,6 @@ int main (int argc, char *argv[]) return EXIT_SUCCESS; } -} \endcode \verbatim diff --git a/src/Doxyfile.in b/src/Doxyfile.in index ae8db1a3..b15d47a6 100644 --- a/src/Doxyfile.in +++ b/src/Doxyfile.in @@ -229,13 +229,7 @@ TAB_SIZE = 2 # "Side Effects:". You can put \n's in the value part of an alias to insert # newlines. -ALIASES = - -# This tag can be used to specify a number of word-keyword mappings (TCL only). -# A mapping has the form "name=value". For example adding "class=itcl::class" -# will allow you to use the command class in the itcl::class meaning. - -TCL_SUBST = +ALIASES = gudhi_example_link{2}="@ref \2 \"\1/\2\"" # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources # only. Doxygen will then generate output that is more tailored for C. For @@ -874,7 +868,7 @@ EXCLUDE_SYMBOLS = # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = @CMAKE_SOURCE_DIR@/biblio/ \ +EXAMPLE_PATH = @CMAKE_SOURCE_DIR@ \ @CMAKE_SOURCE_DIR@/data/ \ @GUDHI_DOXYGEN_EXAMPLE_PATH@ @@ -1040,25 +1034,6 @@ USE_HTAGS = NO VERBATIM_HEADERS = YES -# If the CLANG_ASSISTED_PARSING tag is set to YES then doxygen will use the -# clang parser (see: http://clang.llvm.org/) for more accurate parsing at the -# cost of reduced performance. This can be particularly helpful with template -# rich C++ code for which doxygen's built-in parser lacks the necessary type -# information. -# Note: The availability of this option depends on whether or not doxygen was -# generated with the -Duse-libclang=ON option for CMake. -# The default value is: NO. - -CLANG_ASSISTED_PARSING = NO - -# If clang assisted parsing is enabled you can provide the compiler with command -# line options that you would normally use when invoking the compiler. Note that -# the include paths will already be set by doxygen for the files and directories -# specified with INPUT and INCLUDE_PATH. -# This tag requires that the tag CLANG_ASSISTED_PARSING is set to YES. - -CLANG_OPTIONS = - #--------------------------------------------------------------------------- # Configuration options related to the alphabetical class index #--------------------------------------------------------------------------- @@ -1070,13 +1045,6 @@ CLANG_OPTIONS = ALPHABETICAL_INDEX = YES -# The COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns in -# which the alphabetical index list will be split. -# Minimum value: 1, maximum value: 20, default value: 5. -# This tag requires that the tag ALPHABETICAL_INDEX is set to YES. - -COLS_IN_ALPHA_INDEX = 5 - # In case all classes in a project start with a common prefix, all classes will # be put under the same header in the alphabetical index. The IGNORE_PREFIX tag # can be used to specify a prefix (or a list of prefixes) that should be ignored @@ -1775,16 +1743,6 @@ LATEX_BATCHMODE = NO LATEX_HIDE_INDICES = NO -# If the LATEX_SOURCE_CODE tag is set to YES then doxygen will include source -# code with syntax highlighting in the LaTeX output. -# -# Note that which sources are shown also depends on other settings such as -# SOURCE_BROWSER. -# The default value is: NO. -# This tag requires that the tag GENERATE_LATEX is set to YES. - -LATEX_SOURCE_CODE = NO - # The LATEX_BIB_STYLE tag can be used to specify the style to use for the # bibliography, e.g. plainnat, or ieeetr. See # http://en.wikipedia.org/wiki/BibTeX and \cite for more info. @@ -1857,16 +1815,6 @@ RTF_STYLESHEET_FILE = RTF_EXTENSIONS_FILE = -# If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code -# with syntax highlighting in the RTF output. -# -# Note that which sources are shown also depends on other settings such as -# SOURCE_BROWSER. -# The default value is: NO. -# This tag requires that the tag GENERATE_RTF is set to YES. - -RTF_SOURCE_CODE = NO - #--------------------------------------------------------------------------- # Configuration options related to the man page output #--------------------------------------------------------------------------- @@ -1956,15 +1904,6 @@ GENERATE_DOCBOOK = NO DOCBOOK_OUTPUT = docbook -# If the DOCBOOK_PROGRAMLISTING tag is set to YES, doxygen will include the -# program listings (including syntax highlighting and cross-referencing -# information) to the DOCBOOK output. Note that enabling this will significantly -# increase the size of the DOCBOOK output. -# The default value is: NO. -# This tag requires that the tag GENERATE_DOCBOOK is set to YES. - -DOCBOOK_PROGRAMLISTING = NO - #--------------------------------------------------------------------------- # Configuration options for the AutoGen Definitions output #--------------------------------------------------------------------------- @@ -2139,12 +2078,6 @@ EXTERNAL_GROUPS = YES EXTERNAL_PAGES = YES -# The PERL_PATH should be the absolute path and name of the perl script -# interpreter (i.e. the result of 'which perl'). -# The default file (with absolute path) is: /usr/bin/perl. - -PERL_PATH = /usr/bin/perl - #--------------------------------------------------------------------------- # Configuration options related to the dot tool #--------------------------------------------------------------------------- @@ -2158,15 +2091,6 @@ PERL_PATH = /usr/bin/perl CLASS_DIAGRAMS = NO -# You can define message sequence charts within doxygen comments using the \msc -# command. Doxygen will then run the mscgen tool (see: -# http://www.mcternan.me.uk/mscgen/)) to produce the chart and insert it in the -# documentation. The MSCGEN_PATH tag allows you to specify the directory where -# the mscgen tool resides. If left empty the tool is assumed to be found in the -# default search path. - -MSCGEN_PATH = - # You can include diagrams made with dia in doxygen documentation. Doxygen will # then run dia to produce the diagram and insert it in the documentation. The # DIA_PATH tag allows you to specify the directory where the dia binary resides. diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h index a6098860..e1ab7cb3 100644 --- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -24,7 +24,7 @@ namespace cover_complex { * Visualizations of the simplicial complexes can be done with either * neato (from <a target="_blank" href="http://www.graphviz.org/">graphviz</a>), * <a target="_blank" href="http://www.geomview.org/">geomview</a>, - * <a target="_blank" href="https://github.com/MLWave/kepler-mapper">KeplerMapper</a>. + * <a target="_blank" href="https://github.com/scikit-tda/kepler-mapper">KeplerMapper</a>. * Input point clouds are assumed to be \ref FileFormatsOFF "OFF files" * * \section covers Covers diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index a3613d0d..94579564 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -131,8 +131,7 @@ namespace persistent_cohomology { We provide several example files: run these examples with -h for details on their use, and read the README file. -\li <a href="rips_persistence_8cpp-example.html"> -Rips_complex/rips_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its persistence +\li \gudhi_example_link{Rips_complex,rips_persistence.cpp} computes the Rips complex of a point cloud and outputs its persistence diagram. \code $> ./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3 \endcode \code The complex contains 177838 simplices @@ -144,12 +143,10 @@ diagram. More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. -\li <a href="rips_multifield_persistence_8cpp-example.html"> -Persistent_cohomology/rips_multifield_persistence.cpp</a> computes the Rips complex of a point cloud and outputs its +\li \gudhi_example_link{Persistent_cohomology,rips_multifield_persistence.cpp} computes the Rips complex of a point cloud and outputs its persistence diagram with a family of field coefficients. -\li <a href="rips_distance_matrix_persistence_8cpp-example.html"> -Rips_complex/rips_distance_matrix_persistence.cpp</a> computes the Rips complex of a distance matrix and +\li \gudhi_example_link{Rips_complex,rips_distance_matrix_persistence.cpp} computes the Rips complex of a distance matrix and outputs its persistence diagram. The file should contain square or lower triangular distance matrix with semicolons as separators. @@ -158,8 +155,7 @@ Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. -\li <a href="rips_correlation_matrix_persistence_8cpp-example.html"> -Rips_complex/rips_correlation_matrix_persistence.cpp</a> +\li \gudhi_example_link{Rips_complex,rips_correlation_matrix_persistence.cpp} computes the Rips complex of a correlation matrix and outputs its persistence diagram. Note that no check is performed if the matrix given as the input is a correlation matrix. @@ -169,8 +165,7 @@ Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv More details on the <a href="../../ripscomplex/">Rips complex utilities</a> dedicated page. -\li <a href="alpha_complex_3d_persistence_8cpp-example.html"> -Alpha_complex/alpha_complex_3d_persistence.cpp</a> computes the persistent homology with +\li \gudhi_example_link{Alpha_complex,alpha_complex_3d_persistence.cpp} computes the persistent homology with \f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. \code $> ./alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45 \endcode \code Simplex_tree dim: 3 @@ -235,8 +230,7 @@ Note that the lengths of the sides of the periodic cuboid have to be the same.<b 3 2 36.8838 inf 3 3 58.6783 inf \endcode -\li <a href="alpha_complex_persistence_8cpp-example.html"> -Alpha_complex/alpha_complex_persistence.cpp</a> computes the persistent homology with +\li \gudhi_example_link{Alpha_complex,alpha_complex_persistence.cpp} computes the persistent homology with \f$\mathbb{Z}/p\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. \code $> ./alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off \endcode \code Alpha complex is of dimension 3 - 9273 simplices - 300 vertices. @@ -248,8 +242,7 @@ Simplex_tree dim: 3 More details on the <a href="../../alphacomplex/">Alpha complex utilities</a> dedicated page. -\li <a href="plain_homology_8cpp-example.html"> -Persistent_cohomology/plain_homology.cpp</a> computes the plain homology of a simple simplicial complex without +\li \gudhi_example_link{Persistent_cohomology,plain_homology.cpp} computes the plain homology of a simple simplicial complex without filtration values. */ diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 3888ec8f..cd77b327 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -63,9 +63,8 @@ namespace rips_complex { * value set with \f$max(filtration(4,5), filtration(4,6), filtration(5,6))\f$. * And so on for simplex (0,1,2,3). * - * If the Rips_complex interfaces are not detailed enough for your need, please refer to - * <a href="rips_persistence_step_by_step_8cpp-example.html"> - * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and + * If the Rips_complex interfaces are not detailed enough for your need, please refer to the example + * \gudhi_example_link{Persistent_cohomology,rips_persistence_step_by_step.cpp} , where the constructions of the graph and * the Simplex_tree are more detailed. * * \section sparserips Sparse Rips complex diff --git a/src/Simplex_tree/doc/Intro_simplex_tree.h b/src/Simplex_tree/doc/Intro_simplex_tree.h index ef8dec91..2d3ecdec 100644 --- a/src/Simplex_tree/doc/Intro_simplex_tree.h +++ b/src/Simplex_tree/doc/Intro_simplex_tree.h @@ -39,11 +39,9 @@ namespace Gudhi { * \subsubsection filteredcomplexessimplextreeexamples Examples * * Here is a list of simplex tree examples : - * \li <a href="simple_simplex_tree_8cpp-example.html"> - * Simplex_tree/simple_simplex_tree.cpp</a> - Simple simplex tree construction and basic function use. + * \li \gudhi_example_link{Simplex_tree,simple_simplex_tree.cpp} - Simple simplex tree construction and basic function use. * - * \li <a href="simplex_tree_from_cliques_of_graph_8cpp-example.html"> - * Simplex_tree/simplex_tree_from_cliques_of_graph.cpp</a> - Simplex tree construction from cliques of graph read in + * \li \gudhi_example_link{Simplex_tree,simplex_tree_from_cliques_of_graph.cpp} - Simplex tree construction from cliques of graph read in * a file. * * Simplex tree construction with \f$\mathbb{Z}/3\mathbb{Z}\f$ coefficients on weighted graph Klein bottle file: @@ -54,12 +52,10 @@ Expand the simplex tree in 3.8e-05 s. Information of the Simplex Tree: Number of vertices = 10 Number of simplices = 98 \endcode * - * \li <a href="example_alpha_shapes_3_simplex_tree_from_off_file_8cpp-example.html"> - * Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp</a> - Simplex tree is computed and displayed + * \li \gudhi_example_link{Simplex_tree,example_alpha_shapes_3_simplex_tree_from_off_file.cpp} - Simplex tree is computed and displayed * from a 3D alpha complex (Requires CGAL, GMP and GMPXX to be installed). * - * \li <a href="graph_expansion_with_blocker_8cpp-example.html"> - * Simplex_tree/graph_expansion_with_blocker.cpp</a> - Simple simplex tree construction from a one-skeleton graph with + * \li \gudhi_example_link{Simplex_tree,graph_expansion_with_blocker.cpp} - Simple simplex tree construction from a one-skeleton graph with * a simple blocker expansion method. * * \subsection filteredcomplexeshassecomplex Hasse complex diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 85790baf..6dce947c 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -187,6 +187,12 @@ class Simplex_tree { typedef Simplex_tree_boundary_simplex_iterator<Simplex_tree> Boundary_simplex_iterator; /** \brief Range over the simplices of the boundary of a simplex. */ typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range; + /** \brief Iterator over the simplices of the boundary of a simplex and their opposite vertices. + * + * 'value_type' is std::pair<Simplex_handle, Vertex_handle>. */ + typedef Simplex_tree_boundary_opposite_vertex_simplex_iterator<Simplex_tree> Boundary_opposite_vertex_simplex_iterator; + /** \brief Range over the simplices of the boundary of a simplex and their opposite vertices. */ + typedef boost::iterator_range<Boundary_opposite_vertex_simplex_iterator> Boundary_opposite_vertex_simplex_range; /** \brief Iterator over the simplices of the simplicial complex. * * 'value_type' is Simplex_handle. */ @@ -296,6 +302,23 @@ class Simplex_tree { Boundary_simplex_iterator(this)); } + /** \brief Given a simplex, returns a range over the simplices of its boundary and their opposite vertices. + * + * The boundary of a simplex is the set of codimension \f$1\f$ subsimplices of the simplex. + * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation induced by \f$ v_0 < \cdots < v_d \f$, the + * iterator enumerates the simplices of the boundary in the order: + * \f$[v_0,\cdots,\widehat{v_i},\cdots,v_d]\f$ for \f$i\f$ from \f$d\f$ to \f$0\f$, where \f$\widehat{v_i}\f$ means + * that the vertex \f$v_i\f$, known as the opposite vertex, is omitted from boundary, but returned as the second + * element of a pair. + * + * @param[in] sh Simplex for which the boundary is computed. + */ + template<class SimplexHandle> + Boundary_opposite_vertex_simplex_range boundary_opposite_vertex_simplex_range(SimplexHandle sh) { + return Boundary_opposite_vertex_simplex_range(Boundary_opposite_vertex_simplex_iterator(this, sh), + Boundary_opposite_vertex_simplex_iterator(this)); + } + /** @} */ // end range and iterator methods /** \name Constructor/Destructor * @{ */ @@ -1060,8 +1083,8 @@ class Simplex_tree { * * Inserts all vertices and edges given by a OneSkeletonGraph. * OneSkeletonGraph must be a model of - * <a href="http://www.boost.org/doc/libs/1_76_0/libs/graph/doc/VertexAndEdgeListGraph.html">boost::VertexAndEdgeListGraph</a> - * and <a href="http://www.boost.org/doc/libs/1_76_0/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>. + * <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/VertexAndEdgeListGraph.html">boost::VertexAndEdgeListGraph</a> + * and <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>. * * The vertex filtration value is accessible through the property tag * vertex_filtration_t. @@ -1283,6 +1306,7 @@ class Simplex_tree { Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t + simplex->second.assign_children(new_sib); std::vector<Vertex_handle> blocked_new_sib_vertex_list; // As all intersections are inserted, we can call the blocker function on all new_sib members for (auto new_sib_member = new_sib->members().begin(); @@ -1305,7 +1329,6 @@ class Simplex_tree { new_sib->members().erase(blocked_new_sib_member); } // ensure recursive call - simplex->second.assign_children(new_sib); siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex); } } else { diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h index ee64a277..394c6ee1 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h @@ -5,6 +5,7 @@ * Copyright (C) 2014 Inria * * Modification(s): + * - 2022/04 Vincent Rouvreau: Add Simplex_tree_boundary_opposite_vertex_simplex_iterator for alpha and cech purpose * - YYYY/MM Author: Description of the modification */ @@ -14,10 +15,10 @@ #include <gudhi/Debug_utils.h> #include <boost/iterator/iterator_facade.hpp> -#include <boost/version.hpp> #include <boost/container/static_vector.hpp> #include <vector> +#include <utility> // for std::pair namespace Gudhi { @@ -124,7 +125,7 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< private: friend class boost::iterator_core_access; -// valid when iterating along the SAME boundary. + // valid when iterating along the SAME boundary. bool equal(Simplex_tree_boundary_simplex_iterator const& other) const { return sh_ == other.sh_; } @@ -179,6 +180,116 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< Simplex_handle sh_; // current Simplex_handle in the boundary SimplexTree * st_; // simplex containing the simplicial complex }; + +/* \brief Iterator over the simplices of the boundary of a simplex and their opposite vertices. + * + * Forward iterator, value_type is std::pair<SimplexTree::Simplex_handle, SimplexTree::Vertex_handle>.*/ +template<class SimplexTree> +class Simplex_tree_boundary_opposite_vertex_simplex_iterator : public boost::iterator_facade< + Simplex_tree_boundary_opposite_vertex_simplex_iterator<SimplexTree>, + std::pair<typename SimplexTree::Simplex_handle, typename SimplexTree::Vertex_handle> const, boost::forward_traversal_tag> { + public: + using Simplex_handle = typename SimplexTree::Simplex_handle; + using Vertex_handle = typename SimplexTree::Vertex_handle; + using Siblings = typename SimplexTree::Siblings; + + // For cython purpose only. The object it initializes should be overwritten ASAP and never used before it is + // overwritten. + Simplex_tree_boundary_opposite_vertex_simplex_iterator() + : sib_(nullptr), + st_(nullptr) { + } + + // any end() iterator + explicit Simplex_tree_boundary_opposite_vertex_simplex_iterator(SimplexTree * st) + : last_(st->null_vertex()), + next_(st->null_vertex()), + sib_(nullptr), + baov_(st->null_simplex(), st->null_vertex()), + st_(st) { + } + + template<class SimplexHandle> + Simplex_tree_boundary_opposite_vertex_simplex_iterator(SimplexTree * st, SimplexHandle sh) + : last_(sh->first), + next_(st->null_vertex()), + sib_(nullptr), + baov_(st->null_simplex(), sh->first), + st_(st) { + // Only check once at the beginning instead of for every increment, as this is expensive. + if (SimplexTree::Options::contiguous_vertices) + GUDHI_CHECK(st_->contiguous_vertices(), "The set of vertices is not { 0, ..., n } without holes"); + Siblings * sib = st->self_siblings(sh); + next_ = sib->parent(); + sib_ = sib->oncles(); + if (sib_ != nullptr) { + if (SimplexTree::Options::contiguous_vertices && sib_->oncles() == nullptr) + // Only relevant for edges + baov_.first = sib_->members_.begin()+next_; + else + baov_.first = sib_->find(next_); + } + } + + private: + friend class boost::iterator_core_access; + + // valid when iterating along the SAME boundary. + bool equal(Simplex_tree_boundary_opposite_vertex_simplex_iterator const& other) const { + return (baov_.first == other.baov_.first); + } + + std::pair<Simplex_handle, Vertex_handle> const& dereference() const { + return baov_; + } + + void increment() { + if (sib_ == nullptr) { + baov_.first = st_->null_simplex(); + return; // ------>> + } + Siblings * for_sib = sib_; + Siblings * new_sib = sib_->oncles(); + auto rit = suffix_.rbegin(); + if (SimplexTree::Options::contiguous_vertices && new_sib == nullptr) { + // We reached the root, use a short-cut to find a vertex. + if (rit == suffix_.rend()) { + baov_.second = baov_.first->first; + // Segment, this vertex is the last boundary simplex + baov_.first = for_sib->members_.begin()+last_; + sib_ = nullptr; + return; + } else { + // Dim >= 2, initial step of the descent + baov_.first = for_sib->members_.begin()+*rit; + for_sib = baov_.first->second.children(); + ++rit; + } + } + for (; rit != suffix_.rend(); ++rit) { + baov_.first = for_sib->find(*rit); + for_sib = baov_.first->second.children(); + } + baov_.first = for_sib->find(last_); // baov_.first points to the right simplex now + suffix_.push_back(next_); + next_ = sib_->parent(); + sib_ = new_sib; + baov_.second = suffix_.back(); + } + + // Most of the storage should be moved to the range, iterators should be light. + Vertex_handle last_; // last vertex of the simplex + Vertex_handle next_; // next vertex to push in suffix_ + // 40 seems a conservative bound on the dimension of a Simplex_tree for now, + // as it would not fit on the biggest hard-drive. + boost::container::static_vector<Vertex_handle, 40> suffix_; + // static_vector still has some overhead compared to a trivial hand-made + // version using std::aligned_storage, or compared to making suffix_ static. + Siblings * sib_; // where the next search will start from + std::pair<Simplex_handle, Vertex_handle> baov_; // a pair containing the current Simplex_handle in the boundary and its opposite vertex + SimplexTree * st_; // simplex containing the simplicial complex +}; + /*---------------------------------------------------------------------------*/ /* \brief Iterator over the simplices of a simplicial complex. * diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h index b53bad29..ae25d290 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h @@ -36,6 +36,7 @@ class Simplex_tree_siblings { template<class T> friend class Simplex_tree_boundary_simplex_iterator; template<class T> friend class Simplex_tree_complex_simplex_iterator; template<class T> friend class Simplex_tree_skeleton_simplex_iterator; + template<class T> friend class Simplex_tree_boundary_opposite_vertex_simplex_iterator; typedef typename SimplexTree::Vertex_handle Vertex_handle; typedef typename SimplexTree::Filtration_value Filtration_value; diff --git a/src/Simplex_tree/test/CMakeLists.txt b/src/Simplex_tree/test/CMakeLists.txt index cf2b0153..25b562e0 100644 --- a/src/Simplex_tree/test/CMakeLists.txt +++ b/src/Simplex_tree/test/CMakeLists.txt @@ -34,3 +34,9 @@ if (TBB_FOUND) target_link_libraries(Simplex_tree_make_filtration_non_decreasing_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Simplex_tree_make_filtration_non_decreasing_test_unit) + +add_executable ( Simplex_tree_graph_expansion_test_unit simplex_tree_graph_expansion_unit_test.cpp ) +if (TBB_FOUND) + target_link_libraries(Simplex_tree_graph_expansion_test_unit ${TBB_LIBRARIES}) +endif() +gudhi_add_boost_test(Simplex_tree_graph_expansion_test_unit) diff --git a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp index 881a06ae..6d63d8ae 100644 --- a/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp @@ -9,33 +9,62 @@ */ #include <iostream> -#include <fstream> -#include <string> -#include <algorithm> -#include <utility> // std::pair, std::make_pair -#include <cmath> // float comparison -#include <limits> -#include <functional> // greater +#include <vector> #define BOOST_TEST_DYN_LINK -#define BOOST_TEST_MODULE "simplex_tree" +#define BOOST_TEST_MODULE "simplex_tree_graph_expansion" #include <boost/test/unit_test.hpp> #include <boost/mpl/list.hpp> -// ^ -// /!\ Nothing else from Simplex_tree shall be included to test includes are well defined. #include "gudhi/Simplex_tree.h" +#include <gudhi/Unitary_tests_utils.h> using namespace Gudhi; typedef boost::mpl::list<Simplex_tree<>, Simplex_tree<Simplex_tree_options_fast_persistence>> list_of_tested_variants; +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_all_is_blocked, typeST, list_of_tested_variants) { + std::clog << "********************************************************************\n"; + std::clog << "simplex_tree_expansion_all_is_blocked\n"; + std::clog << "********************************************************************\n"; + using Simplex_handle = typename typeST::Simplex_handle; + // Construct the Simplex Tree with a 1-skeleton graph example + typeST simplex_tree; + + simplex_tree.insert_simplex({0, 1}, 0.); + simplex_tree.insert_simplex({0, 2}, 1.); + simplex_tree.insert_simplex({0, 3}, 2.); + simplex_tree.insert_simplex({1, 2}, 3.); + simplex_tree.insert_simplex({1, 3}, 4.); + simplex_tree.insert_simplex({2, 3}, 5.); + simplex_tree.insert_simplex({2, 4}, 6.); + simplex_tree.insert_simplex({3, 6}, 7.); + simplex_tree.insert_simplex({4, 5}, 8.); + simplex_tree.insert_simplex({4, 6}, 9.); + simplex_tree.insert_simplex({5, 6}, 10.); + simplex_tree.insert_simplex({6}, 10.); -bool AreAlmostTheSame(float a, float b) { - return std::fabs(a - b) < std::numeric_limits<float>::epsilon(); + typeST stree_copy = simplex_tree; + + simplex_tree.expansion_with_blockers(3, [&](Simplex_handle sh){ return true; }); + + std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; + std::clog << " - dimension " << simplex_tree.dimension() << "\n"; + std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + std::clog << " " << "[" << simplex_tree.filtration(f_simplex) << "] "; + for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) + std::clog << "(" << vertex << ")"; + std::clog << std::endl; + } + + BOOST_CHECK(stree_copy == simplex_tree); } BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, list_of_tested_variants) { + std::clog << "********************************************************************\n"; + std::clog << "simplex_tree_expansion_with_blockers_3\n"; + std::clog << "********************************************************************\n"; using Simplex_handle = typename typeST::Simplex_handle; // Construct the Simplex Tree with a 1-skeleton graph example typeST simplex_tree; @@ -72,9 +101,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, li return result; }); - std::clog << "********************************************************************\n"; - std::clog << "simplex_tree_expansion_with_blockers_3\n"; - std::clog << "********************************************************************\n"; std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; std::clog << " - dimension " << simplex_tree.dimension() << "\n"; std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; @@ -89,15 +115,23 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_3, typeST, li BOOST_CHECK(simplex_tree.dimension() == 3); // {4, 5, 6} shall be blocked BOOST_CHECK(simplex_tree.find({4, 5, 6}) == simplex_tree.null_simplex()); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 4.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 5.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 6.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 6.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), 7.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})), + static_cast<typename typeST::Filtration_value>(4.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})), + static_cast<typename typeST::Filtration_value>(5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})), + static_cast<typename typeST::Filtration_value>(6.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})), + static_cast<typename typeST::Filtration_value>(6.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), + static_cast<typename typeST::Filtration_value>(7.)); } BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, list_of_tested_variants) { + std::clog << "********************************************************************\n"; + std::clog << "simplex_tree_expansion_with_blockers_2\n"; + std::clog << "********************************************************************\n"; using Simplex_handle = typename typeST::Simplex_handle; // Construct the Simplex Tree with a 1-skeleton graph example typeST simplex_tree; @@ -134,9 +168,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, li return result; }); - std::clog << "********************************************************************\n"; - std::clog << "simplex_tree_expansion_with_blockers_2\n"; - std::clog << "********************************************************************\n"; std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; std::clog << " - dimension " << simplex_tree.dimension() << "\n"; std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; @@ -151,14 +182,22 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_blockers_2, typeST, li BOOST_CHECK(simplex_tree.dimension() == 2); // {4, 5, 6} shall be blocked BOOST_CHECK(simplex_tree.find({4, 5, 6}) == simplex_tree.null_simplex()); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 4.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 5.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 6.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 6.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})), + static_cast<typename typeST::Filtration_value>(4.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})), + static_cast<typename typeST::Filtration_value>(5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})), + static_cast<typename typeST::Filtration_value>(6.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})), + static_cast<typename typeST::Filtration_value>(6.)); BOOST_CHECK(simplex_tree.find({0,1,2,3}) == simplex_tree.null_simplex()); } -BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_variants) { +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_with_find_simplex_blockers, typeST, list_of_tested_variants) { + std::clog << "********************************************************************\n"; + std::clog << "simplex_tree_expansion_with_find_simplex_blockers\n"; + std::clog << "********************************************************************\n"; + using Simplex_handle = typename typeST::Simplex_handle; // Construct the Simplex Tree with a 1-skeleton graph example typeST simplex_tree; @@ -175,10 +214,66 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_var simplex_tree.insert_simplex({5, 6}, 10.); simplex_tree.insert_simplex({6}, 10.); - simplex_tree.expansion(3); + simplex_tree.expansion_with_blockers(3, [&](Simplex_handle sh){ + bool result = false; + std::clog << "Blocker on ["; + std::vector<typename typeST::Vertex_handle> simplex; + // User can loop on the vertices from the given simplex_handle i.e. + for (auto vertex : simplex_tree.simplex_vertex_range(sh)) { + // We block the expansion, if the vertex '1' is in the given list of vertices + if (vertex == 1) + result = true; + std::clog << vertex << ", "; + simplex.push_back(vertex); + } + std::clog << "] => " << result << std::endl; + // Not efficient but test it works - required by the python interface + BOOST_CHECK(simplex_tree.find(simplex) == sh); + return result; + }); + + std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; + std::clog << " - dimension " << simplex_tree.dimension() << "\n"; + std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + std::clog << " " << "[" << simplex_tree.filtration(f_simplex) << "] "; + for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) + std::clog << "(" << vertex << ")"; + std::clog << std::endl; + } + + BOOST_CHECK(simplex_tree.num_simplices() == 20); + BOOST_CHECK(simplex_tree.dimension() == 2); + + // {1, 2, 3}, {0, 1, 2} and {0, 1, 3} shall be blocked as it contains vertex 1 + BOOST_CHECK(simplex_tree.find({4, 5, 6}) != simplex_tree.null_simplex()); + BOOST_CHECK(simplex_tree.find({1, 2, 3}) == simplex_tree.null_simplex()); + BOOST_CHECK(simplex_tree.find({0, 2, 3}) != simplex_tree.null_simplex()); + BOOST_CHECK(simplex_tree.find({0, 1, 2}) == simplex_tree.null_simplex()); + BOOST_CHECK(simplex_tree.find({0, 1, 3}) == simplex_tree.null_simplex()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_3, typeST, list_of_tested_variants) { std::clog << "********************************************************************\n"; std::clog << "simplex_tree_expansion_3\n"; std::clog << "********************************************************************\n"; + // Construct the Simplex Tree with a 1-skeleton graph example + typeST simplex_tree; + + simplex_tree.insert_simplex({0, 1}, 0.); + simplex_tree.insert_simplex({0, 2}, 1.); + simplex_tree.insert_simplex({0, 3}, 2.); + simplex_tree.insert_simplex({1, 2}, 3.); + simplex_tree.insert_simplex({1, 3}, 4.); + simplex_tree.insert_simplex({2, 3}, 5.); + simplex_tree.insert_simplex({2, 4}, 6.); + simplex_tree.insert_simplex({3, 6}, 7.); + simplex_tree.insert_simplex({4, 5}, 8.); + simplex_tree.insert_simplex({4, 6}, 9.); + simplex_tree.insert_simplex({5, 6}, 10.); + simplex_tree.insert_simplex({6}, 10.); + + simplex_tree.expansion(3); std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; std::clog << " - dimension " << simplex_tree.dimension() << "\n"; std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; @@ -192,16 +287,25 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion, typeST, list_of_tested_var BOOST_CHECK(simplex_tree.num_simplices() == 24); BOOST_CHECK(simplex_tree.dimension() == 3); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({4,5,6})), 10.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 3.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 4.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 5.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 5.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), 5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({4,5,6})), + static_cast<typename typeST::Filtration_value>(10.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})), + static_cast<typename typeST::Filtration_value>(3.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})), + static_cast<typename typeST::Filtration_value>(4.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})), + static_cast<typename typeST::Filtration_value>(5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})), + static_cast<typename typeST::Filtration_value>(5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2,3})), + static_cast<typename typeST::Filtration_value>(5.)); } BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_variants) { + std::clog << "********************************************************************\n"; + std::clog << "simplex_tree_expansion_2\n"; + std::clog << "********************************************************************\n"; // Construct the Simplex Tree with a 1-skeleton graph example typeST simplex_tree; @@ -220,9 +324,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_v simplex_tree.expansion(2); - std::clog << "********************************************************************\n"; - std::clog << "simplex_tree_expansion_2\n"; - std::clog << "********************************************************************\n"; std::clog << "* The complex contains " << simplex_tree.num_simplices() << " simplices"; std::clog << " - dimension " << simplex_tree.dimension() << "\n"; std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; @@ -236,10 +337,15 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_expansion_2, typeST, list_of_tested_v BOOST_CHECK(simplex_tree.num_simplices() == 23); BOOST_CHECK(simplex_tree.dimension() == 2); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({4,5,6})), 10.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,2})), 3.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,1,3})), 4.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({0,2,3})), 5.)); - BOOST_CHECK(AreAlmostTheSame(simplex_tree.filtration(simplex_tree.find({1,2,3})), 5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({4,5,6})), + static_cast<typename typeST::Filtration_value>(10.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,2})), + static_cast<typename typeST::Filtration_value>(3.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,1,3})), + static_cast<typename typeST::Filtration_value>(4.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({0,2,3})), + static_cast<typename typeST::Filtration_value>(5.)); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(simplex_tree.find({1,2,3})), + static_cast<typename typeST::Filtration_value>(5.)); BOOST_CHECK(simplex_tree.find({0,1,2,3}) == simplex_tree.null_simplex()); } diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index bdd41d34..b18e2ec4 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -17,6 +17,8 @@ #include <limits> #include <functional> // greater #include <tuple> // std::tie +#include <iterator> // for std::distance +#include <cstddef> // for std::size_t #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE "simplex_tree" @@ -991,3 +993,48 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_reset_filtration, typeST, list_of_tes } +BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_boundaries_and_opposite_vertex_iterator, typeST, list_of_tested_variants) { + std::clog << "********************************************************************" << std::endl; + std::clog << "TEST OF BOUNDARIES AND OPPOSITE VERTEX ITERATORS" << std::endl; + typeST st; + + st.insert_simplex_and_subfaces({2, 1, 0}, 3.); + st.insert_simplex_and_subfaces({3, 0}, 2.); + st.insert_simplex_and_subfaces({3, 4, 5}, 3.); + st.insert_simplex_and_subfaces({0, 1, 6, 7}, 4.); + + /* Inserted simplex: */ + /* 1 6 */ + /* o---o */ + /* /X\7/ */ + /* o---o---o---o */ + /* 2 0 3\X/4 */ + /* o */ + /* 5 */ + using Simplex = std::vector<typename typeST::Vertex_handle>; + // simplices must be kept sorted by vertex number for std::vector to use operator== - cf. last BOOST_CHECK + std::vector<Simplex> simplices = {{0, 1, 2}, {0, 3}, {0, 1, 6, 7}, {3, 4, 5}, {3, 5}, {2}}; + for (auto simplex : simplices) { + Simplex opposite_vertices; + for(auto boundary_and_opposite_vertex : st.boundary_opposite_vertex_simplex_range(st.find(simplex))) { + Simplex output; + for (auto vertex : st.simplex_vertex_range(boundary_and_opposite_vertex.first)) { + std::clog << vertex << " "; + output.emplace_back(vertex); + } + std::clog << " - opposite vertex = " << boundary_and_opposite_vertex.second << std::endl; + // Check that boundary simplex + opposite vertex = simplex given as input + output.emplace_back(boundary_and_opposite_vertex.second); + std::sort(output.begin(), output.end()); + BOOST_CHECK(simplex == output); + opposite_vertices.emplace_back(boundary_and_opposite_vertex.second); + } + // Check that the list of all opposite vertices = simplex given as input + // no opposite vertices if simplex given as input is of dimension 1 + std::sort(opposite_vertices.begin(), opposite_vertices.end()); + if (simplex.size() > 1) + BOOST_CHECK(simplex == opposite_vertices); + else + BOOST_CHECK(opposite_vertices.size() == 0); + } +} diff --git a/src/Skeleton_blocker/concept/SkeletonBlockerDS.h b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h index 0c2014bd..23eb3670 100644 --- a/src/Skeleton_blocker/concept/SkeletonBlockerDS.h +++ b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h @@ -29,7 +29,7 @@ struct SkeletonBlockerDS { /** * @brief Root_vertex_handle and Vertex_handle are similar to global and local vertex descriptor - * used in <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a> + * used in <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/subgraph.html">boost subgraphs</a> * and allow to localize a vertex of a subcomplex on its parent root complex. * * In gross, vertices are stored in a vector diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h index 0c0cc624..d091d7dd 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h @@ -28,7 +28,7 @@ namespace skeleton_blocker { */ struct Skeleton_blocker_simple_traits { /** - * @brief Global and local handle similar to <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a>. + * @brief Global and local handle similar to <a href="https://www.boost.org/doc/libs/release/libs/graph/doc/subgraph.html">boost subgraphs</a>. * Vertices are stored in a vector. * For the root simplicial complex, the local and global descriptors are the same. * For a subcomplex L and one of its vertices 'v', the local descriptor of 'v' is its position in diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h index 125c6387..031bcb9c 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h @@ -1071,7 +1071,6 @@ class Skeleton_blocker_complex { /** * Removes all the popable blockers of the complex and delete them. - * @returns the number of popable blockers deleted */ void remove_popable_blockers(); @@ -1103,7 +1102,6 @@ class Skeleton_blocker_complex { public: /** * Remove the star of the edge connecting vertices a and b. - * @returns the number of blocker that have been removed */ void remove_star(Vertex_handle a, Vertex_handle b); diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h index 404f04f9..5db9c2fd 100644..100755 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -39,7 +39,6 @@ bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_han /** * Removes all the popable blockers of the complex and delete them. - * @returns the number of popable blockers deleted */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() { @@ -160,7 +159,6 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_s /** * Remove the star of the edge connecting vertices a and b. - * @returns the number of blocker that have been removed */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) { diff --git a/src/Subsampling/include/gudhi/sparsify_point_set.h b/src/Subsampling/include/gudhi/sparsify_point_set.h index 4571b8f3..b325fe3c 100644 --- a/src/Subsampling/include/gudhi/sparsify_point_set.h +++ b/src/Subsampling/include/gudhi/sparsify_point_set.h @@ -11,12 +11,7 @@ #ifndef SPARSIFY_POINT_SET_H_ #define SPARSIFY_POINT_SET_H_ -#include <boost/version.hpp> -#if BOOST_VERSION < 106600 -# include <boost/function_output_iterator.hpp> -#else # include <boost/iterator/function_output_iterator.hpp> -#endif #include <gudhi/Kd_tree_search.h> #ifdef GUDHI_SUBSAMPLING_PROFILING diff --git a/src/Tangential_complex/benchmark/benchmark_tc.cpp b/src/Tangential_complex/benchmark/benchmark_tc.cpp index 6da1425f..8e7c72ff 100644 --- a/src/Tangential_complex/benchmark/benchmark_tc.cpp +++ b/src/Tangential_complex/benchmark/benchmark_tc.cpp @@ -33,6 +33,7 @@ const std::size_t ONLY_LOAD_THE_FIRST_N_POINTS = 20000000; #include <gudhi/sparsify_point_set.h> #include <gudhi/random_point_generators.h> #include <gudhi/Tangential_complex/utilities.h> +#include <gudhi/Simplex_tree.h> #include <CGAL/assertions_behaviour.h> #include <CGAL/Epick_d.h> diff --git a/src/Tangential_complex/example/example_basic.cpp b/src/Tangential_complex/example/example_basic.cpp index ab35edf0..c50b9b8c 100644 --- a/src/Tangential_complex/example/example_basic.cpp +++ b/src/Tangential_complex/example/example_basic.cpp @@ -1,7 +1,6 @@ #include <gudhi/Tangential_complex.h> #include <gudhi/sparsify_point_set.h> -//#include <gudhi/Fake_simplex_tree.h> - +#include <gudhi/Simplex_tree.h> #include <CGAL/Epick_d.h> #include <CGAL/Random.h> @@ -39,7 +38,6 @@ int main(void) { // Export the TC into a Simplex_tree Gudhi::Simplex_tree<> stree; - //Gudhi::Fake_simplex_tree stree; tc.create_complex(stree); // Display stats about inconsistencies diff --git a/src/Tangential_complex/example/example_with_perturb.cpp b/src/Tangential_complex/example/example_with_perturb.cpp index d0d877ea..e70e2980 100644 --- a/src/Tangential_complex/example/example_with_perturb.cpp +++ b/src/Tangential_complex/example/example_with_perturb.cpp @@ -1,5 +1,6 @@ #include <gudhi/Tangential_complex.h> #include <gudhi/sparsify_point_set.h> +#include <gudhi/Simplex_tree.h> #include <CGAL/Epick_d.h> #include <CGAL/Random.h> diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp index 023c1e1a..a24b9ae2 100644 --- a/src/Tangential_complex/test/test_tangential_complex.cpp +++ b/src/Tangential_complex/test/test_tangential_complex.cpp @@ -14,6 +14,7 @@ #include <gudhi/Tangential_complex.h> #include <gudhi/sparsify_point_set.h> +#include <gudhi/Simplex_tree.h> #include <CGAL/Epick_d.h> #include <CGAL/Random.h> diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index b316740d..6a94d1f5 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -1,12 +1,14 @@ # This files manage third party libraries required by GUDHI -find_package(Boost 1.56.0 QUIET OPTIONAL_COMPONENTS filesystem unit_test_framework program_options) +find_package(Boost 1.66.0 QUIET OPTIONAL_COMPONENTS filesystem unit_test_framework program_options) # Boost_FOUND is not reliable if(NOT Boost_VERSION) message(FATAL_ERROR "NOTICE: This program requires Boost and will not be compiled.") endif(NOT Boost_VERSION) include_directories(${Boost_INCLUDE_DIRS}) +message(STATUS "boost include dirs:" ${Boost_INCLUDE_DIRS}) +message(STATUS "boost library dirs:" ${Boost_LIBRARY_DIRS}) find_package(GMP) if(GMP_FOUND) @@ -17,6 +19,15 @@ if(GMP_FOUND) endif() endif() +# from windows vcpkg eigen 3.4.0#2 : build fails with +# error C2440: '<function-style-cast>': cannot convert from 'Eigen::EigenBase<Derived>::Index' to '__gmp_expr<mpq_t,mpq_t>' +# cf. https://gitlab.com/libeigen/eigen/-/issues/2476 +# Workaround is to compile with '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int' +if (FORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT) + message("++ User explicit demand to force EIGEN_DEFAULT_DENSE_INDEX_TYPE to int") + add_definitions(-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int) +endif() + # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CMAKE_CXX_FLAGS are overwritten. # cf. http://doc.cgal.org/latest/Manual/installation.html#title40 # A workaround is to include(${CGAL_USE_FILE}) before adding "-std=c++11". @@ -90,9 +101,6 @@ add_definitions( -DBOOST_ALL_DYN_LINK ) # problem on Mac with boost_system and boost_thread add_definitions( -DBOOST_SYSTEM_NO_DEPRECATED ) -message(STATUS "boost include dirs:" ${Boost_INCLUDE_DIRS}) -message(STATUS "boost library dirs:" ${Boost_LIBRARY_DIRS}) - if (WITH_GUDHI_PYTHON) # Find the correct Python interpreter. # Can be set with -DPYTHON_EXECUTABLE=/usr/bin/python3 or -DPython_ADDITIONAL_VERSIONS=3 for instance. @@ -181,4 +189,4 @@ if (WITH_GUDHI_PYTHON) endif(NOT SPHINX_PATH) endif(SPHINX_FOUND) endif(PYTHONINTERP_FOUND AND CYTHON_FOUND) -endif (WITH_GUDHI_PYTHON)
\ No newline at end of file +endif (WITH_GUDHI_PYTHON) diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h index 0c4320f6..1634b19e 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -1,6 +1,6 @@ // List of GUDHI examples and utils - Doxygen needs at least a file tag to analyse comments // Generated from scripts/cpp_examples_for_doxygen.py -/*! @file Examples +/*! @file * \section Witness_complex_example_section Witness_complex * @example strong_witness_persistence.cpp * @example weak_witness_persistence.cpp diff --git a/src/common/doc/installation.h b/src/common/doc/installation.h index ef668dfb..24a7fc7a 100644 --- a/src/common/doc/installation.h +++ b/src/common/doc/installation.h @@ -5,8 +5,8 @@ * Examples of GUDHI headers inclusion can be found in \ref utilities. * * \section compiling Compiling - * The library uses c++14 and requires <a target="_blank" href="http://www.boost.org/">Boost</a> ≥ 1.56.0 - * and <a target="_blank" href="https://www.cmake.org/">CMake</a> ≥ 3.5. + * The library uses c++14 and requires <a target="_blank" href="https://www.boost.org/">Boost</a> ≥ 1.66.0 + * and <a target="_blank" href="https://cmake.org/">CMake</a> ≥ 3.5. * It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2015. * * \subsection utilities Utilities and examples @@ -56,10 +56,9 @@ make \endverbatim * The multi-field persistent homology algorithm requires GMP which is a free library for arbitrary-precision * arithmetic, operating on signed integers, rational numbers, and floating point numbers. * - * The following example requires the <a target="_blank" href="http://gmplib.org/">GNU Multiple Precision Arithmetic + * The following example requires the <a target="_blank" href="https://gmplib.org/">GNU Multiple Precision Arithmetic * Library</a> (GMP) and will not be built if GMP is not installed: - * \li <a href="rips_multifield_persistence_8cpp-example.html"> - * Persistent_cohomology/rips_multifield_persistence.cpp</a> + * \li \gudhi_example_link{Persistent_cohomology,rips_multifield_persistence.cpp} * * Having GMP version 4.2 or higher installed is recommended. * @@ -76,179 +75,101 @@ make \endverbatim * * The following examples/utilities require the <a target="_blank" href="http://www.cgal.org/">Computational Geometry Algorithms * Library</a> (CGAL \cite cgal:eb-19b) and will not be built if CGAL version 4.11.0 or higher is not installed: - * \li <a href="example_alpha_shapes_3_simplex_tree_from_off_file_8cpp-example.html"> - * Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp</a> - * \li <a href="strong_witness_persistence_8cpp-example.html"> - * Witness_complex/strong_witness_persistence.cpp</a> - * \li <a href="weak_witness_persistence_8cpp-example.html"> - * Witness_complex/weak_witness_persistence.cpp</a> - * \li <a href="example_strong_witness_complex_off_8cpp-example.html"> - * Witness_complex/example_strong_witness_complex_off.cpp</a> - * \li <a href="example_witness_complex_off_8cpp-example.html"> - * Witness_complex/example_witness_complex_off.cpp</a> - * \li <a href="example_witness_complex_sphere_8cpp-example.html"> - * Witness_complex/example_witness_complex_sphere.cpp</a> - * \li <a href="_alpha_complex_from_off_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="_alpha_complex_from_points_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_points.cpp</a> - * \li <a href="alpha_complex_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_persistence.cpp</a> - * \li <a href="custom_persistence_sort_8cpp-example.html"> - * Persistent_cohomology/custom_persistence_sort.cpp</a> - * \li <a href="alpha_rips_persistence_bottleneck_distance_8cpp-example.html"> - * Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp.cpp</a> - * \li <a href="bottleneck_basic_example_8cpp-example.html"> - * Bottleneck_distance/bottleneck_basic_example.cpp</a> - * \li <a href="bottleneck_distance_8cpp-example.html"> - * Bottleneck_distance/bottleneck_distance.cpp</a> - * \li <a href="_coord_g_i_c_8cpp-example.html"> - * Nerve_GIC/CoordGIC.cpp</a> - * \li <a href="_func_g_i_c_8cpp-example.html"> - * Nerve_GIC/FuncGIC.cpp</a> - * \li <a href="_nerve_8cpp-example.html"> - * Nerve_GIC/Nerve.cpp</a> - * \li <a href="_voronoi_g_i_c_8cpp-example.html"> - * Nerve_GIC/VoronoiGIC.cpp</a> - * \li <a href="example_spatial_searching_8cpp-example.html"> - * Spatial_searching/example_spatial_searching.cpp</a> - * \li <a href="example_choose_n_farthest_points_8cpp-example.html"> - * Subsampling/example_choose_n_farthest_points.cpp</a> - * \li <a href="example_pick_n_random_points_8cpp-example.html"> - * Subsampling/example_pick_n_random_points.cpp</a> - * \li <a href="example_sparsify_point_set_8cpp-example.html"> - * Subsampling/example_sparsify_point_set.cpp</a> - * \li <a href="example_basic_8cpp-example.html"> - * Tangential_complex/example_basic.cpp</a> - * \li <a href="example_with_perturb_8cpp-example.html"> - * Tangential_complex/example_with_perturb.cpp</a> - * \li <a href="_weighted_alpha_complex_3d_from_points_8cpp-example.html"> - * Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp</a> - * \li <a href="alpha_complex_3d_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_3d_persistence.cpp</a> - * \li <a href="_coxeter_triangulation_2manifold_tracing_flat_torus_with_boundary_8cpp-example.html"> - * Coxeter_triangulation/manifold_tracing_flat_torus_with_boundary.cpp</a> + * \li \gudhi_example_link{Simplex_tree,example_alpha_shapes_3_simplex_tree_from_off_file.cpp} + * \li \gudhi_example_link{Witness_complex,strong_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,weak_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,example_strong_witness_complex_off.cpp} + * \li \gudhi_example_link{Witness_complex,example_witness_complex_off.cpp} + * \li \gudhi_example_link{Witness_complex,example_witness_complex_sphere.cpp} + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_off.cpp} + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_points.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_persistence.cpp} + * \li \gudhi_example_link{Persistent_cohomology,custom_persistence_sort.cpp} + * \li \gudhi_example_link{Bottleneck_distance,alpha_rips_persistence_bottleneck_distance.cpp} + * \li \gudhi_example_link{Bottleneck_distance,bottleneck_basic_example.cpp} + * \li \gudhi_example_link{Bottleneck_distance,bottleneck_distance.cpp} + * \li \gudhi_example_link{Nerve_GIC,CoordGIC.cpp} + * \li \gudhi_example_link{Nerve_GIC,FuncGIC.cpp} + * \li \gudhi_example_link{Nerve_GIC,Nerve.cpp} + * \li \gudhi_example_link{Nerve_GIC,VoronoiGIC.cpp} + * \li \gudhi_example_link{Spatial_searching,example_spatial_searching.cpp} + * \li \gudhi_example_link{Subsampling,example_choose_n_farthest_points.cpp} + * \li \gudhi_example_link{Subsampling,example_pick_n_random_points.cpp} + * \li \gudhi_example_link{Subsampling,example_sparsify_point_set.cpp} + * \li \gudhi_example_link{Tangential_complex,example_basic.cpp} + * \li \gudhi_example_link{Tangential_complex,example_with_perturb.cpp} + * \li \gudhi_example_link{Alpha_complex,Weighted_alpha_complex_3d_from_points.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_3d_persistence.cpp} + * \li \gudhi_example_link{Coxeter_triangulation,manifold_tracing_flat_torus_with_boundary.cpp} * * \subsection eigen Eigen * Some GUDHI modules (cf. \ref main_page "modules list"), and few examples require - * <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen</a> is a C++ template library for linear algebra: + * <a target="_blank" href="https://eigen.tuxfamily.org">Eigen</a> is a C++ template library for linear algebra: * matrices, vectors, numerical solvers, and related algorithms. * - * The following examples/utilities require the <a target="_blank" href="http://eigen.tuxfamily.org/">Eigen</a> and will not be + * The following examples/utilities require the <a target="_blank" href="https://eigen.tuxfamily.org">Eigen</a> and will not be * built if Eigen is not installed: - * \li <a href="_alpha_complex_from_off_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="_alpha_complex_from_points_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_points.cpp</a> - * \li <a href="alpha_complex_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_persistence.cpp</a> - * \li <a href="alpha_complex_3d_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_3d_persistence.cpp</a> - * \li <a href="_weighted_alpha_complex_3d_from_points_8cpp-example.html"> - * Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp</a> - * \li <a href="alpha_rips_persistence_bottleneck_distance_8cpp-example.html"> - * Bottleneck_distance/alpha_rips_persistence_bottleneck_distance.cpp.cpp</a> - * \li <a href="custom_persistence_sort_8cpp-example.html"> - * Persistent_cohomology/custom_persistence_sort.cpp</a> - * \li <a href="example_spatial_searching_8cpp-example.html"> - * Spatial_searching/example_spatial_searching.cpp</a> - * \li <a href="example_choose_n_farthest_points_8cpp-example.html"> - * Subsampling/example_choose_n_farthest_points.cpp</a> - * \li <a href="example_pick_n_random_points_8cpp-example.html"> - * Subsampling/example_pick_n_random_points.cpp</a> - * \li <a href="example_sparsify_point_set_8cpp-example.html"> - * Subsampling/example_sparsify_point_set.cpp</a> - * \li <a href="example_basic_8cpp-example.html"> - * Tangential_complex/example_basic.cpp</a> - * \li <a href="example_with_perturb_8cpp-example.html"> - * Tangential_complex/example_with_perturb.cpp</a> - * \li <a href="strong_witness_persistence_8cpp-example.html"> - * Witness_complex/strong_witness_persistence.cpp</a> - * \li <a href="weak_witness_persistence_8cpp-example.html"> - * Witness_complex/weak_witness_persistence.cpp</a> - * \li <a href="example_strong_witness_complex_off_8cpp-example.html"> - * Witness_complex/example_strong_witness_complex_off.cpp</a> - * \li <a href="example_witness_complex_off_8cpp-example.html"> - * Witness_complex/example_witness_complex_off.cpp</a> - * \li <a href="example_witness_complex_sphere_8cpp-example.html"> - * Witness_complex/example_witness_complex_sphere.cpp</a> - * \li <a href="_coxeter_triangulation_2cell_complex_from_basic_circle_manifold_8cpp-example.html"> - * Coxeter_triangulation/cell_complex_from_basic_circle_manifold.cpp</a> - * \li <a href="_coxeter_triangulation_2manifold_tracing_custom_function_8cpp-example.html"> - * Coxeter_triangulation/manifold_tracing_custom_function.cpp</a> - * \li <a href="_coxeter_triangulation_2manifold_tracing_flat_torus_with_boundary_8cpp-example.html"> - * Coxeter_triangulation/manifold_tracing_flat_torus_with_boundary.cpp</a> + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_off.cpp} + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_points.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_persistence.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_3d_persistence.cpp} + * \li \gudhi_example_link{Alpha_complex,Weighted_alpha_complex_3d_from_points.cpp} + * \li \gudhi_example_link{Bottleneck_distance,alpha_rips_persistence_bottleneck_distance.cpp} + * \li \gudhi_example_link{Persistent_cohomology,custom_persistence_sort.cpp} + * \li \gudhi_example_link{Spatial_searching,example_spatial_searching.cpp} + * \li \gudhi_example_link{Subsampling,example_choose_n_farthest_points.cpp} + * \li \gudhi_example_link{Subsampling,example_pick_n_random_points.cpp} + * \li \gudhi_example_link{Subsampling,example_sparsify_point_set.cpp} + * \li \gudhi_example_link{Tangential_complex,example_basic.cpp} + * \li \gudhi_example_link{Tangential_complex,example_with_perturb.cpp} + * \li \gudhi_example_link{Witness_complex,strong_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,weak_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,example_strong_witness_complex_off.cpp} + * \li \gudhi_example_link{Witness_complex,example_witness_complex_off.cpp} + * \li \gudhi_example_link{Witness_complex,example_witness_complex_sphere.cpp} + * \li \gudhi_example_link{Coxeter_triangulation,cell_complex_from_basic_circle_manifold.cpp} + * \li \gudhi_example_link{Coxeter_triangulation,manifold_tracing_custom_function.cpp} + * \li \gudhi_example_link{Coxeter_triangulation,manifold_tracing_flat_torus_with_boundary.cpp} * * \subsection tbb Threading Building Blocks - * <a target="_blank" href="https://www.threadingbuildingblocks.org/">Intel® TBB</a> lets you easily write parallel + * <a target="_blank" href="https://github.com/oneapi-src/oneTBB">Intel® TBB</a> lets you easily write parallel * C++ programs that take full advantage of multicore performance, that are portable and composable, and that have * future-proof scalability. * * Having Intel® TBB installed is recommended to parallelize and accelerate some GUDHI computations. * * The following examples/utilities are using Intel® TBB if installed: - * \li <a href="_alpha_complex_from_off_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_off.cpp</a> - * \li <a href="_alpha_complex_from_points_8cpp-example.html"> - * Alpha_complex/Alpha_complex_from_points.cpp</a> - * \li <a href="alpha_complex_3d_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_3d_persistence.cpp</a> - * \li <a href="alpha_complex_persistence_8cpp-example.html"> - * Alpha_complex/alpha_complex_persistence.cpp</a> - * \li <a href="cubical_complex_persistence_8cpp-example.html"> - * Bitmap_cubical_complex/cubical_complex_persistence.cpp</a> - * \li <a href="periodic_cubical_complex_persistence_8cpp-example.html"> - * Bitmap_cubical_complex/periodic_cubical_complex_persistence.cpp</a> - * \li <a href="_random_bitmap_cubical_complex_8cpp-example.html"> - * Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp</a> - * \li <a href="_coord_g_i_c_8cpp-example.html"> - * Nerve_GIC/CoordGIC.cpp</a> - * \li <a href="_func_g_i_c_8cpp-example.html"> - * Nerve_GIC/FuncGIC.cpp</a> - * \li <a href="_nerve_8cpp-example.html"> - * Nerve_GIC/Nerve.cpp</a> - * \li <a href="_voronoi_g_i_c_8cpp-example.html"> - * Nerve_GIC/VoronoiGIC.cpp</a> - * \li <a href="simple_simplex_tree_8cpp-example.html"> - * Simplex_tree/simple_simplex_tree.cpp</a> - * \li <a href="example_alpha_shapes_3_simplex_tree_from_off_file_8cpp-example.html"> - * Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp</a> - * \li <a href="simplex_tree_from_cliques_of_graph_8cpp-example.html"> - * Simplex_tree/simplex_tree_from_cliques_of_graph.cpp</a> - * \li <a href="graph_expansion_with_blocker_8cpp-example.html"> - * Simplex_tree/graph_expansion_with_blocker.cpp</a> - * \li <a href="alpha_complex_3d_persistence_8cpp-example.html"> - * Persistent_cohomology/alpha_complex_3d_persistence.cpp</a> - * \li <a href="alpha_complex_persistence_8cpp-example.html"> - * Persistent_cohomology/alpha_complex_persistence.cpp</a> - * \li <a href="rips_persistence_via_boundary_matrix_8cpp-example.html"> - * Persistent_cohomology/rips_persistence_via_boundary_matrix.cpp</a> - * \li <a href="persistence_from_file_8cpp-example.html"> - * Persistent_cohomology/persistence_from_file.cpp</a> - * \li <a href="persistence_from_simple_simplex_tree_8cpp-example.html"> - * Persistent_cohomology/persistence_from_simple_simplex_tree.cpp</a> - * \li <a href="plain_homology_8cpp-example.html"> - * Persistent_cohomology/plain_homology.cpp</a> - * \li <a href="rips_multifield_persistence_8cpp-example.html"> - * Persistent_cohomology/rips_multifield_persistence.cpp</a> - * \li <a href="rips_persistence_step_by_step_8cpp-example.html"> - * Persistent_cohomology/rips_persistence_step_by_step.cpp</a> - * \li <a href="custom_persistence_sort_8cpp-example.html"> - * Persistent_cohomology/custom_persistence_sort.cpp</a> - * \li <a href="example_one_skeleton_rips_from_points_8cpp-example.html"> - * Rips_complex/example_one_skeleton_rips_from_points.cpp</a> - * \li <a href="example_rips_complex_from_off_file_8cpp-example.html"> - * Rips_complex/example_rips_complex_from_off_file.cpp</a> - * \li <a href="rips_distance_matrix_persistence_8cpp-example.html"> - * Rips_complex/rips_distance_matrix_persistence.cpp</a> - * \li <a href="rips_persistence_8cpp-example.html"> - * Rips_complex/rips_persistence.cpp</a> - * \li <a href="strong_witness_persistence_8cpp-example.html"> - * Witness_complex/strong_witness_persistence.cpp</a> - * \li <a href="weak_witness_persistence_8cpp-example.html"> - * Witness_complex/weak_witness_persistence.cpp</a> - * \li <a href="example_nearest_landmark_table_8cpp-example.html"> - * Witness_complex/example_nearest_landmark_table.cpp</a> + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_off.cpp} + * \li \gudhi_example_link{Alpha_complex,Alpha_complex_from_points.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_3d_persistence.cpp} + * \li \gudhi_example_link{Alpha_complex,alpha_complex_persistence.cpp} + * \li \gudhi_example_link{Bitmap_cubical_complex,cubical_complex_persistence.cpp} + * \li \gudhi_example_link{Bitmap_cubical_complex,periodic_cubical_complex_persistence.cpp} + * \li \gudhi_example_link{Bitmap_cubical_complex,Random_bitmap_cubical_complex.cpp} + * \li \gudhi_example_link{Nerve_GIC,CoordGIC.cpp} + * \li \gudhi_example_link{Nerve_GIC,FuncGIC.cpp} + * \li \gudhi_example_link{Nerve_GIC,Nerve.cpp} + * \li \gudhi_example_link{Nerve_GIC,VoronoiGIC.cpp} + * \li \gudhi_example_link{Simplex_tree,simple_simplex_tree.cpp} + * \li \gudhi_example_link{Simplex_tree,example_alpha_shapes_3_simplex_tree_from_off_file.cpp} + * \li \gudhi_example_link{Simplex_tree,simplex_tree_from_cliques_of_graph.cpp} + * \li \gudhi_example_link{Simplex_tree,graph_expansion_with_blocker.cpp} + * \li \gudhi_example_link{Persistent_cohomology,alpha_complex_3d_persistence.cpp} + * \li \gudhi_example_link{Persistent_cohomology,alpha_complex_persistence.cpp} + * \li \gudhi_example_link{Persistent_cohomology,rips_persistence_via_boundary_matrix.cpp} + * \li \gudhi_example_link{Persistent_cohomology,persistence_from_file.cpp} + * \li \gudhi_example_link{Persistent_cohomology,persistence_from_simple_simplex_tree.cpp} + * \li \gudhi_example_link{Persistent_cohomology,plain_homology.cpp} + * \li \gudhi_example_link{Persistent_cohomology,rips_multifield_persistence.cpp} + * \li \gudhi_example_link{Persistent_cohomology,rips_persistence_step_by_step.cpp} + * \li \gudhi_example_link{Persistent_cohomology,custom_persistence_sort.cpp} + * \li \gudhi_example_link{Rips_complex,example_one_skeleton_rips_from_points.cpp} + * \li \gudhi_example_link{Rips_complex,example_rips_complex_from_off_file.cpp} + * \li \gudhi_example_link{Rips_complex,rips_distance_matrix_persistence.cpp} + * \li \gudhi_example_link{Rips_complex,rips_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,strong_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,weak_witness_persistence.cpp} + * \li \gudhi_example_link{Witness_complex,example_nearest_landmark_table.cpp} * * \section Contributions Bug reports and contributions * Please help us improving the quality of the GUDHI library. diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index 6f995fee..2cb02e3f 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -231,13 +231,12 @@ homology of the input sequence. The resulting method is simple and extremely efficient. Computation of edge collapse and persistent homology of a filtered flag complex via edge collapse as described in - \cite edgecollapsesocg2020. + \cite edgecollapsearxiv. </td> <td width="15%"> - <b>Author:</b> Siddharth Pritam<br> + <b>Author:</b> Siddharth Pritam, Marc Glisse<br> <b>Introduced in:</b> GUDHI 3.3.0<br> - <b>Copyright:</b> MIT<br> - <b>Requires:</b> \ref eigen + <b>Copyright:</b> MIT </td> </tr> <tr> diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h index a1b104e2..29d5423d 100644 --- a/src/common/include/gudhi/reader_utils.h +++ b/src/common/include/gudhi/reader_utils.h @@ -14,11 +14,7 @@ #include <gudhi/graph_simplicial_complex.h> #include <gudhi/Debug_utils.h> -#if BOOST_VERSION < 106600 -# include <boost/function_output_iterator.hpp> -#else # include <boost/iterator/function_output_iterator.hpp> -#endif #include <boost/graph/adjacency_list.hpp> #include <iostream> diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 8eb7478e..54221151 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -148,10 +148,6 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("Eigen3 version ${EIGEN3_VERSION}") # No problem, even if no CGAL found set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_EIGEN3', ") - set(GUDHI_USE_EIGEN3 "True") - else (EIGEN3_FOUND) - set(GUDHI_USE_EIGEN3 "False") endif (EIGEN3_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") @@ -180,6 +176,16 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") endif () + # from windows vcpkg eigen 3.4.0#2 : build fails with + # error C2440: '<function-style-cast>': cannot convert from 'Eigen::EigenBase<Derived>::Index' to '__gmp_expr<mpq_t,mpq_t>' + # cf. https://gitlab.com/libeigen/eigen/-/issues/2476 + # Workaround is to compile with '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int' + if (FORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int', ") + endif() + + + add_gudhi_debug_info("Boost version ${Boost_VERSION}") if(CGAL_FOUND) if(NOT CGAL_VERSION VERSION_LESS 5.3.0) # CGAL_HEADER_ONLY has been dropped for CGAL >= 5.3. Only the header-only version is supported. @@ -214,13 +220,14 @@ if(PYTHONINTERP_FOUND) endif(NOT GMP_LIBRARIES_DIR) add_GUDHI_PYTHON_lib_dir(${GMP_LIBRARIES_DIR}) message("** Add gmp ${GMP_LIBRARIES_DIR}") + # When FORCE_CGAL_NOT_TO_BUILD_WITH_GMPXX is set, not defining CGAL_USE_GMPXX is sufficient enough if(GMPXX_FOUND) add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ") add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}") add_GUDHI_PYTHON_lib_dir(${GMPXX_LIBRARIES_DIR}) message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}") - endif(GMPXX_FOUND) + endif() endif(GMP_FOUND) if(MPFR_FOUND) add_gudhi_debug_info("MPFR_LIBRARIES = ${MPFR_LIBRARIES}") @@ -264,6 +271,10 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_INCLUDE_DIRS "${GUDHI_PYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ") endif() + if(DEBUG_TRACES) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DDEBUG_TRACES', ") + endif(DEBUG_TRACES) + if(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) set( GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}") endif(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) @@ -580,6 +591,10 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_dtm_rips_complex) endif() + # persistence graphical tools + if(MATPLOTLIB_FOUND) + add_gudhi_py_test(test_persistence_graphical_tools) + endif() # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") @@ -590,4 +605,4 @@ if(PYTHONINTERP_FOUND) else(PYTHONINTERP_FOUND) message("++ Python module will not be compiled because no Python interpreter was found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") -endif(PYTHONINTERP_FOUND)
\ No newline at end of file +endif(PYTHONINTERP_FOUND) diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index cfd22742..9e67d38a 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -27,7 +27,8 @@ Remarks If you pass :code:`precision = 'exact'` to :func:`~gudhi.AlphaComplex.__init__`, the filtration values are the exact ones converted to float. This can be very slow. If you pass :code:`precision = 'safe'` (the default), the filtration values are only - guaranteed to have a small multiplicative error compared to the exact value. + guaranteed to have a small multiplicative error compared to the exact value, see + :func:`~gudhi.AlphaComplex.set_float_relative_precision` to modify the precision. A drawback, when computing persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval [10^12,10^12+10^6]. Using :code:`precision = 'fast'` makes the computations slightly faster, and the combinatorics are still exact, but @@ -177,11 +178,11 @@ Weighted version ^^^^^^^^^^^^^^^^ A weighted version for Alpha complex is available. It is like a usual Alpha complex, but based on a -`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#title20>`_. +`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT>`_. This example builds the weighted alpha-complex of a small molecule, where atoms have different sizes. It is taken from -`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13>`_. +`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes>`_. Then, it is asked to display information about the alpha complex. diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index 35c344e3..cff84691 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -33,25 +33,19 @@ Compiling These instructions are for people who want to compile gudhi from source, they are unnecessary if you installed a binary package of Gudhi as above. They assume that you have downloaded a `release <https://github.com/GUDHI/gudhi-devel/releases>`_, -with a name like `gudhi.3.2.0.tar.gz`, then run `tar xf gudhi.3.2.0.tar.gz`, which -created a directory `gudhi.3.2.0`, hereinafter referred to as `/path-to-gudhi/`. +with a name like `gudhi.3.X.Y.tar.gz`, then run `tar xf gudhi.3.X.Y.tar.gz`, which +created a directory `gudhi.3.X.Y`, hereinafter referred to as `/path-to-gudhi/`. If you are instead using a git checkout, beware that the paths are a bit different, and in particular the `python/` subdirectory is actually `src/python/` there. -The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.56.0, +The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.66.0, `CMake <https://www.cmake.org/>`_ :math:`\geq` 3.5 to generate makefiles, -`NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ and -`pybind11 <https://github.com/pybind/pybind11>`_ to compile -the GUDHI Python module. -It is a multi-platform library and compiles on Linux, Mac OSX and Visual -Studio 2017 or later. +Python :math:`\geq` 3.5, `NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ +:math:`\geq` 0.27 and `pybind11 <https://github.com/pybind/pybind11>`_ to compile the GUDHI Python module. +It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2017 or later. -On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python -:math:`\geq` 3.5 are available because of the required Visual Studio version. - -On other systems, if you have several Python/python installed, the version 2.X -will be used by default, but you can force it by adding +If you have several Python/python installed, the version 2.X may be used by default, but you can force it by adding :code:`-DPython_ADDITIONAL_VERSIONS=3` to the cmake command. GUDHI Python module compilation @@ -142,54 +136,63 @@ If :code:`import gudhi` succeeds, please have a look to debug information: .. code-block:: python - import gudhi - print(gudhi.__debug_info__) + import gudhi as gd + print(gd.__debug_info__) + print("+ Installed modules are: " + gd.__available_modules) + print("+ Missing modules are: " + gd.__missing_modules) You shall have something like: .. code-block:: none - Python version 2.7.15 - Cython version 0.26.1 - Numpy version 1.14.1 - Eigen3 version 3.1.1 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;reader_utils;witness_complex; - strong_witness_complex;alpha_complex; - Missing modules are: bottleneck_distance;nerve_gic;subsampling; - tangential_complex;persistence_graphical_tools; - euclidean_witness_complex;euclidean_strong_witness_complex; - CGAL version 4.7.1000 - GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so - GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so - TBB version 9107 found and used + Pybind11 version 2.8.1 + Python version 3.7.12 + Cython version 0.29.25 + Numpy version 1.21.4 + Boost version 1.77.0 + + Installed modules are: off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex; + + Missing modules are: bottleneck;nerve_gic;subsampling;tangential_complex;alpha_complex;euclidean_witness_complex; + euclidean_strong_witness_complex; -Here, you can see that bottleneck_distance, nerve_gic, subsampling and -tangential_complex are missing because of the CGAL version. -persistence_graphical_tools is not available as matplotlib is not -available. +Here, you can see that the modules that need CGAL are missing, because CGAL is not installed. +:code:`persistence_graphical_tools` is installed, but +`its functions <https://gudhi.inria.fr/python/latest/persistence_graphical_tools_ref.html>`_ will produce an error as +matplotlib is not available. Unitary tests cannot be run as pytest is missing. A complete configuration would be : .. code-block:: none - Python version 3.6.5 - Cython version 0.28.2 - Pytest version 3.3.2 - Matplotlib version 2.2.2 - Numpy version 1.14.5 - Eigen3 version 3.3.4 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;persistence_graphical_tools; - reader_utils;witness_complex;strong_witness_complex; - persistence_graphical_tools;bottleneck_distance;nerve_gic;subsampling; - tangential_complex;alpha_complex;euclidean_witness_complex; - euclidean_strong_witness_complex; - CGAL header only version 4.11.0 + Pybind11 version 2.8.1 + Python version 3.9.7 + Cython version 0.29.24 + Pytest version 6.2.5 + Matplotlib version 3.5.0 + Numpy version 1.21.4 + Scipy version 1.7.3 + Scikit-learn version 1.0.1 + POT version 0.8.0 + HNSWlib found + PyKeOps version [pyKeOps]: 1.5 + EagerPy version 0.30.0 + TensorFlow version 2.7.0 + Sphinx version 4.3.0 + Sphinx-paramlinks version 0.5.2 + python_docs_theme found + Eigen3 version 3.4.0 + Boost version 1.74.0 + CGAL version 5.3 GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so + MPFR_LIBRARIES = /usr/lib/x86_64-linux-gnu/libmpfr.so TBB version 9107 found and used + + Installed modules are: bottleneck;off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex;nerve_gic;subsampling; + tangential_complex;alpha_complex;euclidean_witness_complex;euclidean_strong_witness_complex; + + Missing modules are: + Documentation ============= @@ -345,8 +348,8 @@ You can still deactivate LaTeX rendering by saying: .. code-block:: python - import gudhi - gudhi.persistence_graphical_tools._gudhi_matplotlib_use_tex=False + import gudhi as gd + gd.persistence_graphical_tools._gudhi_matplotlib_use_tex=False PyKeOps ------- diff --git a/src/python/doc/nerve_gic_complex_user.rst b/src/python/doc/nerve_gic_complex_user.rst index 0b820abf..8633cadb 100644 --- a/src/python/doc/nerve_gic_complex_user.rst +++ b/src/python/doc/nerve_gic_complex_user.rst @@ -12,7 +12,7 @@ Definition Visualizations of the simplicial complexes can be done with either neato (from `graphviz <http://www.graphviz.org/>`_), `geomview <http://www.geomview.org/>`_, -`KeplerMapper <https://github.com/MLWave/kepler-mapper>`_. +`KeplerMapper <https://github.com/scikit-tda/kepler-mapper>`_. Input point clouds are assumed to be OFF files (cf. `OFF file format <fileformats.html#off-file-format>`_). Covers diff --git a/src/python/doc/persistence_graphical_tools_user.rst b/src/python/doc/persistence_graphical_tools_user.rst index d95b9d2b..e1d28c71 100644 --- a/src/python/doc/persistence_graphical_tools_user.rst +++ b/src/python/doc/persistence_graphical_tools_user.rst @@ -60,7 +60,7 @@ of shape (N x 2) encoding a persistence diagram (in a given dimension). import matplotlib.pyplot as plt import gudhi import numpy as np - d = np.array([[0, 1], [1, 2], [1, np.inf]]) + d = np.array([[0., 1.], [1., 2.], [1., np.inf]]) gudhi.plot_persistence_diagram(d) plt.show() diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in index 3043201a..79e12fbc 100644 --- a/src/python/gudhi/__init__.py.in +++ b/src/python/gudhi/__init__.py.in @@ -23,10 +23,6 @@ __all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@] __available_modules = '' __missing_modules = '' -# For unitary tests purpose -# could use "if 'collapse_edges' in gudhi.__all__" when collapse edges will have a python module -__GUDHI_USE_EIGEN3 = @GUDHI_USE_EIGEN3@ - # Try to import * from gudhi.__module_name for default modules. # Extra modules require an explicit import by the user (mostly because of # unusual dependencies, but also to avoid cluttering namespace gudhi and diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index a4888914..375e1561 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -31,6 +31,10 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + + @staticmethod + void set_float_relative_precision(double precision) nogil + @staticmethod + double get_float_relative_precision() nogil # AlphaComplex python interface cdef class AlphaComplex: @@ -133,3 +137,28 @@ cdef class AlphaComplex: self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, mas, compute_filtration) return stree + + @staticmethod + def set_float_relative_precision(precision): + """ + :param precision: When the AlphaComplex is constructed with :code:`precision = 'safe'` (the default), + one can set the float relative precision of filtration values computed in + :func:`~gudhi.AlphaComplex.create_simplex_tree`. + Default is :code:`1e-5` (cf. :func:`~gudhi.AlphaComplex.get_float_relative_precision`). + For more details, please refer to + `CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double <https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html>`_ + :type precision: float + """ + if precision <=0. or precision >= 1.: + raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1") + Alpha_complex_interface.set_float_relative_precision(precision) + + @staticmethod + def get_float_relative_precision(): + """ + :returns: The float relative precision of filtration values computation in + :func:`~gudhi.AlphaComplex.create_simplex_tree` when the AlphaComplex is constructed with + :code:`precision = 'safe'` (the default). + :rtype: float + """ + return Alpha_complex_interface.get_float_relative_precision() diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py index 848dc03e..7ed11360 100644 --- a/src/python/gudhi/persistence_graphical_tools.py +++ b/src/python/gudhi/persistence_graphical_tools.py @@ -12,6 +12,9 @@ from os import path from math import isfinite import numpy as np from functools import lru_cache +import warnings +import errno +import os from gudhi.reader_utils import read_persistence_intervals_in_dimension from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension @@ -22,6 +25,7 @@ __license__ = "MIT" _gudhi_matplotlib_use_tex = True + def __min_birth_max_death(persistence, band=0.0): """This function returns (min_birth, max_death) from the persistence. @@ -44,20 +48,46 @@ def __min_birth_max_death(persistence, band=0.0): min_birth = float(interval[1][0]) if band > 0.0: max_death += band + # can happen if only points at inf death + if min_birth == max_death: + max_death = max_death + 1.0 return (min_birth, max_death) def _array_handler(a): - ''' + """ :param a: if array, assumes it is a (n x 2) np.array and return a persistence-compatible list (padding with 0), so that the plot can be performed seamlessly. - ''' - if isinstance(a[0][1], np.float64) or isinstance(a[0][1], float): + """ + if isinstance(a[0][1], (np.floating, float)): return [[0, x] for x in a] else: return a + +def _limit_to_max_intervals(persistence, max_intervals, key): + """This function returns truncated persistence if length is bigger than max_intervals. + :param persistence: Persistence intervals values list. Can be grouped by dimension or not. + :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death). + :param max_intervals: maximal number of intervals to display. + Selected intervals are those with the longest life time. Set it + to 0 to see all. Default value is 1000. + :type max_intervals: int. + :param key: key function for sort algorithm. + :type key: function or lambda. + """ + if max_intervals > 0 and max_intervals < len(persistence): + warnings.warn( + "There are %s intervals given as input, whereas max_intervals is set to %s." + % (len(persistence), max_intervals) + ) + # Sort by life time, then takes only the max_intervals elements + return sorted(persistence, key=key, reverse=True)[:max_intervals] + else: + return persistence + + @lru_cache(maxsize=1) def _matplotlib_can_use_tex(): """This function returns True if matplotlib can deal with LaTeX, False otherwise. @@ -65,17 +95,17 @@ def _matplotlib_can_use_tex(): """ try: from matplotlib import checkdep_usetex + return checkdep_usetex(True) - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_barcode( persistence=[], persistence_file="", alpha=0.6, - max_intervals=1000, - max_barcodes=1000, + max_intervals=20000, inf_delta=0.1, legend=False, colormap=None, @@ -97,7 +127,7 @@ def plot_persistence_barcode( :type alpha: float. :param max_intervals: maximal number of intervals to display. Selected intervals are those with the longest life time. Set it - to 0 to see all. Default value is 1000. + to 0 to see all. Default value is 20000. :type max_intervals: int. :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x inf_delta)` above :code:`max_death` value. A reasonable value is @@ -119,99 +149,68 @@ def plot_persistence_barcode( import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if path.isfile(persistence_file): # Reset persistence persistence = [] - diag = read_persistence_intervals_grouped_by_dimension( - persistence_file=persistence_file - ) + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) for key in diag.keys(): for persistence_interval in diag[key]: persistence.append((key, persistence_interval)) else: - print("file " + persistence_file + " not found.") - return None - - persistence = _array_handler(persistence) - - if max_barcodes != 1000: - print("Deprecated parameter. It has been replaced by max_intervals") - max_intervals = max_barcodes - - if max_intervals > 0 and max_intervals < len(persistence): - # Sort by life time, then takes only the max_intervals elements - persistence = sorted( - persistence, - key=lambda life_time: life_time[1][1] - life_time[1][0], - reverse=True, - )[:max_intervals] - - if colormap == None: - colormap = plt.cm.Set1.colors - if axes == None: - fig, axes = plt.subplots(1, 1) + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) - persistence = sorted(persistence, key=lambda birth: birth[1][0]) + try: + persistence = _array_handler(persistence) + persistence = _limit_to_max_intervals( + persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + (min_birth, max_death) = __min_birth_max_death(persistence) + persistence = sorted(persistence, key=lambda birth: birth[1][0]) + except IndexError: + min_birth, max_death = 0.0, 1.0 + pass - (min_birth, max_death) = __min_birth_max_death(persistence) - ind = 0 delta = (max_death - min_birth) * inf_delta # Replace infinity values with max_death + delta for bar code to be more # readable infinity = max_death + delta axis_start = min_birth - delta - # Draw horizontal bars in loop - for interval in reversed(persistence): - if float(interval[1][1]) != float("inf"): - # Finite death case - axes.barh( - ind, - (interval[1][1] - interval[1][0]), - height=0.8, - left=interval[1][0], - alpha=alpha, - color=colormap[interval[0]], - linewidth=0, - ) - else: - # Infinite death case for diagram to be nicer - axes.barh( - ind, - (infinity - interval[1][0]), - height=0.8, - left=interval[1][0], - alpha=alpha, - color=colormap[interval[0]], - linewidth=0, - ) - ind = ind + 1 + + if axes == None: + _, axes = plt.subplots(1, 1) + if colormap == None: + colormap = plt.cm.Set1.colors + + x=[birth for (dim,(birth,death)) in persistence] + y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence] + c=[colormap[dim] for (dim,(birth,death)) in persistence] + + axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0) if legend: dimensions = list(set(item[0] for item in persistence)) axes.legend( - handles=[ - mpatches.Patch(color=colormap[dim], label=str(dim)) - for dim in dimensions - ], - loc="lower right", + handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right", ) axes.set_title("Persistence barcode", fontsize=fontsize) # Ends plot on infinity value and starts a little bit before min_birth - axes.axis([axis_start, infinity, 0, ind]) + if len(x) != 0: + axes.axis([axis_start, infinity, 0, len(x)]) return axes - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_diagram( @@ -219,14 +218,13 @@ def plot_persistence_diagram( persistence_file="", alpha=0.6, band=0.0, - max_intervals=1000, - max_plots=1000, + max_intervals=1000000, inf_delta=0.1, legend=False, colormap=None, axes=None, fontsize=16, - greyblock=True + greyblock=True, ): """This function plots the persistence diagram from persistence values list, a np.array of shape (N x 2) representing a diagram in a single @@ -244,7 +242,7 @@ def plot_persistence_diagram( :type band: float. :param max_intervals: maximal number of intervals to display. Selected intervals are those with the longest life time. Set it - to 0 to see all. Default value is 1000. + to 0 to see all. Default value is 1000000. :type max_intervals: int. :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x inf_delta)` above :code:`max_death` value. A reasonable value is @@ -268,47 +266,35 @@ def plot_persistence_diagram( import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if path.isfile(persistence_file): # Reset persistence persistence = [] - diag = read_persistence_intervals_grouped_by_dimension( - persistence_file=persistence_file - ) + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) for key in diag.keys(): for persistence_interval in diag[key]: persistence.append((key, persistence_interval)) else: - print("file " + persistence_file + " not found.") - return None - - persistence = _array_handler(persistence) + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) - if max_plots != 1000: - print("Deprecated parameter. It has been replaced by max_intervals") - max_intervals = max_plots - - if max_intervals > 0 and max_intervals < len(persistence): - # Sort by life time, then takes only the max_intervals elements - persistence = sorted( - persistence, - key=lambda life_time: life_time[1][1] - life_time[1][0], - reverse=True, - )[:max_intervals] - - if colormap == None: - colormap = plt.cm.Set1.colors - if axes == None: - fig, axes = plt.subplots(1, 1) + try: + persistence = _array_handler(persistence) + persistence = _limit_to_max_intervals( + persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + min_birth, max_death = __min_birth_max_death(persistence, band) + except IndexError: + min_birth, max_death = 0.0, 1.0 + pass - (min_birth, max_death) = __min_birth_max_death(persistence, band) delta = (max_death - min_birth) * inf_delta # Replace infinity values with max_death + delta for diagram to be more # readable @@ -316,61 +302,56 @@ def plot_persistence_diagram( axis_end = max_death + delta / 2 axis_start = min_birth - delta + if axes == None: + _, axes = plt.subplots(1, 1) + if colormap == None: + colormap = plt.cm.Set1.colors # bootstrap band if band > 0.0: x = np.linspace(axis_start, infinity, 1000) axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red") # lower diag patch if greyblock: - axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], fill=True, color='lightgrey')) - # Draw points in loop - pts_at_infty = False # Records presence of pts at infty - for interval in reversed(persistence): - if float(interval[1][1]) != float("inf"): - # Finite death case - axes.scatter( - interval[1][0], - interval[1][1], - alpha=alpha, - color=colormap[interval[0]], + axes.add_patch( + mpatches.Polygon( + [[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], + fill=True, + color="lightgrey", ) - else: - pts_at_infty = True - # Infinite death case for diagram to be nicer - axes.scatter( - interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]] - ) - if pts_at_infty: + ) + # line display of equation : birth = death + axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k") + + x=[birth for (dim,(birth,death)) in persistence] + y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence] + c=[colormap[dim] for (dim,(birth,death)) in persistence] + + axes.scatter(x,y,alpha=alpha,color=c) + if float("inf") in (death for (dim,(birth,death)) in persistence): # infinity line and text - axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k") axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha) # Infinity label yt = axes.get_yticks() - yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity + yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity yt = np.append(yt, infinity) ytl = ["%.3f" % e for e in yt] # to avoid float precision error - ytl[-1] = r'$+\infty$' + ytl[-1] = r"$+\infty$" axes.set_yticks(yt) axes.set_yticklabels(ytl) if legend: dimensions = list(set(item[0] for item in persistence)) - axes.legend( - handles=[ - mpatches.Patch(color=colormap[dim], label=str(dim)) - for dim in dimensions - ] - ) + axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions]) axes.set_xlabel("Birth", fontsize=fontsize) axes.set_ylabel("Death", fontsize=fontsize) axes.set_title("Persistence diagram", fontsize=fontsize) # Ends plot on infinity value and starts a little bit before min_birth - axes.axis([axis_start, axis_end, axis_start, infinity + delta/2]) + axes.axis([axis_start, axis_end, axis_start, infinity + delta / 2]) return axes - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_density( @@ -384,7 +365,7 @@ def plot_persistence_density( legend=False, axes=None, fontsize=16, - greyblock=False + greyblock=False, ): """This function plots the persistence density from persistence values list, np.array of shape (N x 2) representing a diagram @@ -444,12 +425,13 @@ def plot_persistence_density( import matplotlib.patches as mpatches from scipy.stats import kde from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if dimension is None: @@ -460,10 +442,16 @@ def plot_persistence_density( persistence_file=persistence_file, only_this_dim=dimension ) else: - print("file " + persistence_file + " not found.") - return None + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) + + # default cmap value cannot be done at argument definition level as matplotlib is not yet defined. + if cmap is None: + cmap = plt.cm.hot_r + if axes == None: + _, axes = plt.subplots(1, 1) - if len(persistence) > 0: + try: + # if not read from file but given by an argument persistence = _array_handler(persistence) persistence_dim = np.array( [ @@ -472,47 +460,54 @@ def plot_persistence_density( if (dim_interval[0] == dimension) or (dimension is None) ] ) - - persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])] - if max_intervals > 0 and max_intervals < len(persistence_dim): - # Sort by life time, then takes only the max_intervals elements + persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])] persistence_dim = np.array( - sorted( - persistence_dim, - key=lambda life_time: life_time[1] - life_time[0], - reverse=True, - )[:max_intervals] + _limit_to_max_intervals( + persistence_dim, max_intervals, key=lambda life_time: life_time[1] - life_time[0] + ) ) - # Set as numpy array birth and death (remove undefined values - inf and NaN) - birth = persistence_dim[:, 0] - death = persistence_dim[:, 1] - - # default cmap value cannot be done at argument definition level as matplotlib is not yet defined. - if cmap is None: - cmap = plt.cm.hot_r - if axes == None: - fig, axes = plt.subplots(1, 1) + # Set as numpy array birth and death (remove undefined values - inf and NaN) + birth = persistence_dim[:, 0] + death = persistence_dim[:, 1] + birth_min = birth.min() + birth_max = birth.max() + death_min = death.min() + death_max = death.max() + + # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents + k = kde.gaussian_kde([birth, death], bw_method=bw_method) + xi, yi = np.mgrid[ + birth_min : birth_max : nbins * 1j, death_min : death_max : nbins * 1j, + ] + zi = k(np.vstack([xi.flatten(), yi.flatten()])) + # Make the plot + img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto") + plot_success = True + + # IndexError on empty diagrams, ValueError on only inf death values + except (IndexError, ValueError): + birth_min = 0.0 + birth_max = 1.0 + death_min = 0.0 + death_max = 1.0 + plot_success = False + pass # line display of equation : birth = death - x = np.linspace(death.min(), birth.max(), 1000) + x = np.linspace(death_min, birth_max, 1000) axes.plot(x, x, color="k", linewidth=1.0) - # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents - k = kde.gaussian_kde([birth, death], bw_method=bw_method) - xi, yi = np.mgrid[ - birth.min() : birth.max() : nbins * 1j, - death.min() : death.max() : nbins * 1j, - ] - zi = k(np.vstack([xi.flatten(), yi.flatten()])) - - # Make the plot - img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap) - if greyblock: - axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey')) + axes.add_patch( + mpatches.Polygon( + [[birth_min, birth_min], [death_max, birth_min], [death_max, death_max]], + fill=True, + color="lightgrey", + ) + ) - if legend: + if plot_success and legend: plt.colorbar(img, ax=axes) axes.set_xlabel("Birth", fontsize=fontsize) @@ -521,7 +516,5 @@ def plot_persistence_density( return axes - except ImportError: - print( - "This function is not available, you may be missing matplotlib and/or scipy." - ) + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 006a24ed..4f229663 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -45,6 +45,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>": Simplex_tree_interface_full_featured() nogil + Simplex_tree_interface_full_featured(Simplex_tree_interface_full_featured&) nogil double simplex_filtration(vector[int] simplex) nogil void assign_simplex_filtration(vector[int] simplex, double filtration) nogil void initialize_filtration() nogil @@ -65,6 +66,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except + void reset_filtration(double filtration, int dimension) nogil + bint operator==(Simplex_tree_interface_full_featured) nogil # Iterators over Simplex tree pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil @@ -74,6 +76,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + + # Expansion with blockers + ctypedef bool (*blocker_func_t)(vector[int], void *user_data) + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index c3720936..2c53a872 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -16,6 +16,9 @@ __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +cdef bool callback(vector[int] simplex, void *blocker_func): + return (<object>blocker_func)(simplex) + # SimplexTree python interface cdef class SimplexTree: """The simplex tree is an efficient and flexible data structure for @@ -38,13 +41,29 @@ cdef class SimplexTree: cdef Simplex_tree_persistence_interface * pcohptr # Fake constructor that does nothing but documenting the constructor - def __init__(self): + def __init__(self, other = None): """SimplexTree constructor. + + :param other: If `other` is `None` (default value), an empty `SimplexTree` is created. + If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`. + :type other: SimplexTree (Optional) + :returns: An empty or a copy simplex tree. + :rtype: SimplexTree + + :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`. + :note: If the `SimplexTree` is a copy, the persistence information is not copied. If you need it in the clone, + you have to call :func:`compute_persistence` on it even if you had already computed it in the original. """ # The real cython constructor - def __cinit__(self): - self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) + def __cinit__(self, other = None): + if other: + if isinstance(other, SimplexTree): + self.thisptr = _get_copy_intptr(other) + else: + raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.") + else: + self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) def __dealloc__(self): cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() @@ -63,6 +82,21 @@ cdef class SimplexTree: """ return self.pcohptr != NULL + def copy(self): + """ + :returns: A simplex tree that is a deep copy of itself. + :rtype: SimplexTree + + :note: The persistence information is not copied. If you need it in the clone, you have to call + :func:`compute_persistence` on it even if you had already computed it in the original. + """ + stree = SimplexTree() + stree.thisptr = _get_copy_intptr(self) + return stree + + def __deepcopy__(self): + return self.copy() + def filtration(self, simplex): """This function returns the filtration value for a given N-simplex in this simplicial complex, or +infinity if it is not in the complex. @@ -443,6 +477,27 @@ cdef class SimplexTree: persistence_result = self.pcohptr.get_persistence() return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence) + def expansion_with_blocker(self, max_dim, blocker_func): + """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added + incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func` + returns `True` for this simplex. + + The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a + filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func` + with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed, + otherwise it is kept. The algorithm then proceeds with the next candidate. + + .. warning:: + Several candidates of the same dimension may be inserted simultaneously before calling `block_simplex`, so + if you examine the complex in `block_simplex`, you may hit a few simplices of the same dimension that have + not been vetted by `block_simplex` yet, or have already been rejected but not yet removed. + + :param max_dim: Expansion maximal dimension value. + :type max_dim: int + :param blocker_func: Blocker oracle. + :type blocker_func: Callable[[List[int]], bool] + """ + self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func) def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """This function computes and returns the persistence of the simplicial complex. @@ -621,18 +676,17 @@ cdef class SimplexTree: return (normal0, normals, infinite0, infinites) def collapse_edges(self, nb_iterations = 1): - """Assuming the simplex tree is a 1-skeleton graph, this method collapse edges (simplices of higher dimension - are ignored) and resets the simplex tree from the remaining edges. - A good candidate is to build a simplex tree on top of a :class:`~gudhi.RipsComplex` of dimension 1 before - collapsing edges + """Assuming the complex is a graph (simplices of higher dimension are ignored), this method implicitly + interprets it as the 1-skeleton of a flag complex, and replaces it with another (smaller) graph whose + expansion has the same persistent homology, using a technique known as edge collapses + (see :cite:`edgecollapsearxiv`). + + A natural application is to get a simplex tree of dimension 1 from :class:`~gudhi.RipsComplex`, + then collapse edges, perform :meth:`expansion()` and finally compute persistence (cf. :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`). - For implementation details, please refer to :cite:`edgecollapsesocg2020`. :param nb_iterations: The number of edge collapse iterations to perform. Default is 1. :type nb_iterations: int - - :note: collapse_edges method requires `Eigen <installation.html#eigen>`_ >= 3.1.0 and an exception is thrown - if this method is not available. """ # Backup old pointer cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() @@ -642,3 +696,13 @@ cdef class SimplexTree: self.thisptr = <intptr_t>(ptr.collapse_edges(nb_iter)) # Delete old pointer del ptr + + def __eq__(self, other:SimplexTree): + """Test for structural equality + :returns: True if the 2 simplex trees are equal, False otherwise. + :rtype: bool + """ + return dereference(self.get_ptr()) == dereference(other.get_ptr()) + +cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil: + return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr()))) diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 0541572b..16f63c3d 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -12,9 +12,11 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ Class to generate a weighted Rips complex from a distance matrix and weights on vertices, - in the way described in :cite:`dtmfiltrations`. + in the way described in :cite:`dtmfiltrations` with `p=1`. The filtration value of vertex `i` is `2*weights[i]`, + and the filtration value of edge `ij` is `distance_matrix[i][j]+weights[i]+weights[j]`, + or the maximum of the filtrations of its extremities, whichever is largest. Remark that all the filtration values are doubled compared to the definition in the paper - for the consistency with RipsComplex. + for consistency with RipsComplex. """ def __init__(self, distance_matrix, diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 671af4a4..469b91ce 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -57,6 +57,16 @@ class Alpha_complex_interface { alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } + static void set_float_relative_precision(double precision) { + // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h + CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::set_relative_precision_of_to_double(precision); + } + + static double get_float_relative_precision() { + // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h + return CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::get_relative_precision_of_to_double(); + } + private: std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; }; diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 629f6083..7f9b0067 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -15,9 +15,7 @@ #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> -#ifdef GUDHI_USE_EIGEN3 #include <gudhi/Flag_complex_edge_collapser.h> -#endif #include <iostream> #include <vector> @@ -42,6 +40,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator; + typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); public: @@ -164,7 +163,6 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { } Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) { -#ifdef GUDHI_USE_EIGEN3 using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>; std::vector<Filtered_edge> edges; for (Simplex_handle sh : Base::skeleton_simplex_range(1)) { @@ -178,7 +176,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { } for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) { - edges = Gudhi::collapse::flag_complex_collapse_edges(edges); + edges = Gudhi::collapse::flag_complex_collapse_edges(std::move(edges)); } Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(); // Copy the original 0-skeleton @@ -190,9 +188,13 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge)); } return collapsed_stree_ptr; -#else - throw std::runtime_error("Unable to collapse edges as it requires Eigen3 >= 3.1.0."); -#endif + } + + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) { + Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){ + Simplex simplex(Base::simplex_vertex_range(sh).begin(), Base::simplex_vertex_range(sh).end()); + return user_func(simplex, user_data); + }); } // Iterator over the simplex tree diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml index a9fb4985..55b64466 100644 --- a/src/python/pyproject.toml +++ b/src/python/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools", "wheel", "numpy>=1.15.0", "cython", "pybind11"] +requires = ["setuptools>=24.2.0", "wheel", "numpy>=1.15.0", "cython>=0.27", "pybind11"] build-backend = "setuptools.build_meta" diff --git a/src/python/setup.py.in b/src/python/setup.py.in index 23746998..2c67c2c5 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -5,6 +5,7 @@ Copyright (C) 2019 Inria Modification(s): + - 2021/12 Vincent Rouvreau: Python 3.5 as minimal version - YYYY/MM Author: Description of the modification """ @@ -43,7 +44,7 @@ for module in cython_modules: include_dirs=include_dirs, runtime_library_dirs=runtime_library_dirs,)) -ext_modules = cythonize(ext_modules, compiler_directives={'language_level': str(sys.version_info[0])}) +ext_modules = cythonize(ext_modules, compiler_directives={'language_level': '3'}) for module in pybind11_modules: my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)] @@ -86,6 +87,7 @@ setup( long_description_content_type='text/x-rst', long_description=long_description, ext_modules = ext_modules, + python_requires='>=3.5.0', install_requires = ['numpy >= 1.15.0',], package_data={"": ["*.dll"], }, ) diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index f15284f3..f81e6137 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -286,3 +286,30 @@ def _weighted_doc_example(precision): def test_weighted_doc_example(): for precision in ['fast', 'safe', 'exact']: _weighted_doc_example(precision) + +def test_float_relative_precision(): + assert AlphaComplex.get_float_relative_precision() == 1e-5 + # Must be > 0. + with pytest.raises(ValueError): + AlphaComplex.set_float_relative_precision(0.) + # Must be < 1. + with pytest.raises(ValueError): + AlphaComplex.set_float_relative_precision(1.) + + points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + st = AlphaComplex(points=points).create_simplex_tree() + filtrations = list(st.get_filtration()) + + # Get a better precision + AlphaComplex.set_float_relative_precision(1e-15) + assert AlphaComplex.get_float_relative_precision() == 1e-15 + + st = AlphaComplex(points=points).create_simplex_tree() + filtrations_better_resolution = list(st.get_filtration()) + + assert len(filtrations) == len(filtrations_better_resolution) + for idx in range(len(filtrations)): + # check simplex is the same + assert filtrations[idx][0] == filtrations_better_resolution[idx][0] + # check filtration is about the same with a relative precision of the worst case + assert filtrations[idx][1] == pytest.approx(filtrations_better_resolution[idx][1], rel=1e-5) diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py new file mode 100644 index 00000000..c19836b7 --- /dev/null +++ b/src/python/test/test_persistence_graphical_tools.py @@ -0,0 +1,121 @@ +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Vincent Rouvreau + + Copyright (C) 2021 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +import gudhi as gd +import numpy as np +import matplotlib as plt +import pytest + + +def test_array_handler(): + diags = np.array([[1, 2], [3, 4], [5, 6]], float) + arr_diags = gd.persistence_graphical_tools._array_handler(diags) + for idx in range(len(diags)): + assert arr_diags[idx][0] == 0 + np.testing.assert_array_equal(arr_diags[idx][1], diags[idx]) + + diags = [(1.0, 2.0), (3.0, 4.0), (5.0, 6.0)] + arr_diags = gd.persistence_graphical_tools._array_handler(diags) + for idx in range(len(diags)): + assert arr_diags[idx][0] == 0 + assert arr_diags[idx][1] == diags[idx] + + diags = [(0, (1.0, 2.0)), (0, (3.0, 4.0)), (0, (5.0, 6.0))] + assert gd.persistence_graphical_tools._array_handler(diags) == diags + + +def test_min_birth_max_death(): + diags = [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + (0, (0.0, 0.118398)), + (0, (0.118398, 1.0)), + (0, (0.0, 0.117908)), + (0, (0.0, 0.112307)), + (0, (0.0, 0.107535)), + (0, (0.0, 0.106382)), + ] + assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (0.0, 1.0) + assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (0.0, 5.0) + + +def test_limit_min_birth_max_death(): + diags = [ + (0, (2.0, float("inf"))), + (0, (2.0, float("inf"))), + ] + assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (2.0, 3.0) + assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (2.0, 6.0) + + +def test_limit_to_max_intervals(): + diags = [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + (0, (0.0, 0.118398)), + (0, (0.118398, 1.0)), + (0, (0.0, 0.117908)), + (0, (0.0, 0.112307)), + (0, (0.0, 0.107535)), + (0, (0.0, 0.106382)), + ] + # check no warnings if max_intervals equals to the diagrams number + with pytest.warns(None) as record: + truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals( + diags, 10, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + # check diagrams are not sorted + assert truncated_diags == diags + assert len(record) == 0 + + # check warning if max_intervals lower than the diagrams number + with pytest.warns(UserWarning) as record: + truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals( + diags, 5, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + # check diagrams are truncated and sorted by life time + assert truncated_diags == [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.118398, 1.0)), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + ] + assert len(record) == 1 + + +def _limit_plot_persistence(function): + pplot = function(persistence=[]) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[], legend=True) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[(0, float("inf"))]) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[(0, float("inf"))], legend=True) + assert isinstance(pplot, plt.axes.SubplotBase) + + +def test_limit_plot_persistence(): + for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]: + _limit_plot_persistence(function) + + +def _non_existing_persistence_file(function): + with pytest.raises(FileNotFoundError): + function(persistence_file="pouetpouettralala.toubiloubabdou") + + +def test_non_existing_persistence_file(): + for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]: + _non_existing_persistence_file(function) diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 31c46213..688f4fd6 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import SimplexTree, __GUDHI_USE_EIGEN3 +from gudhi import SimplexTree import numpy as np import pytest @@ -354,16 +354,11 @@ def test_collapse_edges(): assert st.num_simplices() == 10 - if __GUDHI_USE_EIGEN3: - st.collapse_edges() - assert st.num_simplices() == 9 - assert st.find([1, 3]) == False - for simplex in st.get_skeleton(0): - assert simplex[1] == 1. - else: - # If no Eigen3, collapse_edges throws an exception - with pytest.raises(RuntimeError): - st.collapse_edges() + st.collapse_edges() + assert st.num_simplices() == 9 + assert st.find([0, 2]) == False # [1, 3] would be fine as well + for simplex in st.get_skeleton(0): + assert simplex[1] == 1. def test_reset_filtration(): st = SimplexTree() @@ -447,4 +442,116 @@ def test_persistence_intervals_in_dimension(): assert np.array_equal(H2, np.array([[ 0., float("inf")]])) # Test empty case assert st.persistence_intervals_in_dimension(3).shape == (0, 2) -
\ No newline at end of file + +def test_equality_operator(): + st1 = SimplexTree() + st2 = SimplexTree() + + assert st1 == st2 + + st1.insert([1,2,3], 4.) + assert st1 != st2 + + st2.insert([1,2,3], 4.) + assert st1 == st2 + +def test_simplex_tree_deep_copy(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = st.copy() + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_deep_copy_constructor(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = SimplexTree(st) + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_constructor_exception(): + with pytest.raises(TypeError): + st = SimplexTree(other = "Construction from a string shall raise an exception") + +def test_expansion_with_blocker(): + st=SimplexTree() + st.insert([0,1],0) + st.insert([0,2],1) + st.insert([0,3],2) + st.insert([1,2],3) + st.insert([1,3],4) + st.insert([2,3],5) + st.insert([2,4],6) + st.insert([3,6],7) + st.insert([4,5],8) + st.insert([4,6],9) + st.insert([5,6],10) + st.insert([6],10) + + def blocker(simplex): + try: + # Block all simplices that countains vertex 6 + simplex.index(6) + print(simplex, ' is blocked') + return True + except ValueError: + print(simplex, ' is accepted') + st.assign_filtration(simplex, st.filtration(simplex) + 1.) + return False + + st.expansion_with_blocker(2, blocker) + assert st.num_simplices() == 22 + assert st.dimension() == 2 + assert st.find([4,5,6]) == False + assert st.filtration([0,1,2]) == 4. + assert st.filtration([0,1,3]) == 5. + assert st.filtration([0,2,3]) == 6. + assert st.filtration([1,2,3]) == 6. + + st.expansion_with_blocker(3, blocker) + assert st.num_simplices() == 23 + assert st.dimension() == 3 + assert st.find([4,5,6]) == False + assert st.filtration([0,1,2]) == 4. + assert st.filtration([0,1,3]) == 5. + assert st.filtration([0,2,3]) == 6. + assert st.filtration([1,2,3]) == 6. + assert st.filtration([0,1,2,3]) == 7. |