diff options
Diffstat (limited to 'src')
38 files changed, 1605 insertions, 118 deletions
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index b12c9690..e76edc5f 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -23,7 +23,7 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "alpha.pers") + "-p" "2" "-m" "0.45" "-o" "safe.pers") add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" @@ -31,13 +31,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "safe.pers" "-s") + "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") if (DIFF_PATH) add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH} - "exact.pers" "alpha.pers") + "exact.pers" "safe.pers") add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH} - "safe.pers" "alpha.pers") + "fast.pers" "safe.pers") endif() add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index 19e608ad..09c84eb3 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -94,11 +94,11 @@ int main(int argc, char **argv) { int coeff_field_characteristic = 0; Filtration_value min_persistence = 0.; bool exact_version = false; - bool safe_version = false; + bool fast_version = false; bool weighted_version = false; bool periodic_version = false; - program_options(argc, argv, off_file_points, exact_version, safe_version, weight_file, cuboid_file, output_file_diag, + program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, cuboid_file, output_file_diag, alpha_square_max_value, coeff_field_characteristic, min_persistence); std::vector<double> weights; @@ -120,16 +120,16 @@ int main(int argc, char **argv) { periodic_version = true; } - Gudhi::alpha_complex::complexity complexity = Gudhi::alpha_complex::complexity::FAST; + Gudhi::alpha_complex::complexity complexity = Gudhi::alpha_complex::complexity::SAFE; if (exact_version) { - if (safe_version) { - std::cerr << "You cannot set the exact and the safe version." << std::endl; + if (fast_version) { + std::cerr << "You cannot set the exact and the fast version." << std::endl; exit(-1); } complexity = Gudhi::alpha_complex::complexity::EXACT; } - if (safe_version) { - complexity = Gudhi::alpha_complex::complexity::SAFE; + if (fast_version) { + complexity = Gudhi::alpha_complex::complexity::FAST; } Simplex_tree simplex_tree; @@ -258,7 +258,7 @@ int main(int argc, char **argv) { return 0; } -void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &safe, +void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, std::string &weight_file, std::string &cuboid_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence) { @@ -270,9 +270,9 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "exact,e", po::bool_switch(&exact), - "To activate exact version of Alpha complex 3d (default is false, not available if safe is set)")( - "safe,s", po::bool_switch(&safe), - "To activate safe version of Alpha complex 3d (default is false, not available if exact is set)")( + "To activate exact version of Alpha complex 3d (default is false, not available if fast is set)")( + "fast,f", po::bool_switch(&fast), + "To activate fast version of Alpha complex 3d (default is false, not available if exact is set)")( "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()), "Name of file containing a point weights. Format is one weight per line:\n W1\n ...\n Wn ")( "cuboid-file,c", po::value<std::string>(&cuboid_file), @@ -303,7 +303,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool std::cout << std::endl; std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; std::cout << "of a 3D Alpha complex defined on a set of input points.\n"; - std::cout << "3D Alpha complex can be exact or safe, weighted and/or periodic\n\n"; + std::cout << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\n\n"; std::cout << "The output diagram contains one bar per line, written with the convention: \n"; std::cout << " p dim b d \n"; std::cout << "where dim is the dimension of the homological feature,\n"; diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index 98f56802..50a39d32 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -107,6 +107,7 @@ It must be in the format described points (one value per line).
* `-e [ --exact ]` for the exact computation version (not compatible with
weight and periodic version).
+* `-f [ --fast ]` for the fast computation version.
**Example**
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b40d506a..40fdcf2b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -35,6 +35,7 @@ add_gudhi_module(Skeleton_blocker) add_gudhi_module(Spatial_searching) add_gudhi_module(Subsampling) add_gudhi_module(Tangential_complex) +add_gudhi_module(Toplex_map) add_gudhi_module(Witness_complex) add_gudhi_module(Nerve_GIC) diff --git a/src/Doxyfile.in b/src/Doxyfile.in index 858a9299..1c293d1c 100644 --- a/src/Doxyfile.in +++ b/src/Doxyfile.in @@ -785,7 +785,8 @@ EXCLUDE = data/ \ GudhUI/ \ cmake/ \ src/cython/ \ - include/gudhi_patches/ + include/gudhi_patches/ \ + README.md # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded diff --git a/src/Rips_complex/concept/SimplicialComplexForRips.h b/src/Rips_complex/concept/SimplicialComplexForRips.h index 3c5acecf..36ab1b0c 100644 --- a/src/Rips_complex/concept/SimplicialComplexForRips.h +++ b/src/Rips_complex/concept/SimplicialComplexForRips.h @@ -34,6 +34,9 @@ struct SimplicialComplexForRips { /** \brief Type used to store the filtration values of the simplicial complex. */ typedef unspecified Filtration_value; + /** \brief Handle type to a simplex contained in the simplicial complex. */ + typedef unspecified Simplex_handle; + /** \brief Inserts a given `Gudhi::rips_complex::Rips_complex::OneSkeletonGraph` in the simplicial complex. */ template<class OneSkeletonGraph> void insert_graph(const OneSkeletonGraph& skel_graph); @@ -42,6 +45,24 @@ struct SimplicialComplexForRips { * explained in \ref ripsdefinition. */ void expansion(int max_dim); + /** \brief Expands a simplicial complex 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 `block_simplex` + * returns true for this simplex. + * + * @param[in] max_dim Expansion maximal dimension value. + * @param[in] block_simplex Blocker oracle. Its concept is <CODE>bool block_simplex(Simplex_handle sh)</CODE> + * + * 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 + * `block_simplex` on a `Simplex_handle` for this new simplex. If `block_simplex` returns true, the simplex is + * removed, otherwise it is kept. + */ + template< typename Blocker > + void expansion_with_blockers(int max_dim, Blocker block_simplex); + + /** \brief Returns a range over the vertices of a simplex. */ + unspecified simplex_vertex_range(Simplex_handle sh); + /** \brief Returns the number of vertices in the simplicial complex. */ std::size_t num_vertices(); diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index a2537036..97d66fbd 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -92,21 +92,27 @@ namespace rips_complex { * The sparse Rips filtration was introduced by Don Sheehy * \cite sheehy13linear. We are using the version described in * \cite buchet16efficient (except that we multiply all filtration values - * by 2, to match the usual Rips complex), which proves a - * \f$\frac{1+\epsilon}{1-\epsilon}\f$-interleaving, although in practice the + * by 2, to match the usual Rips complex), for which \cite cavanna15geometric proves a + * \f$(1,\frac{1}{1-\epsilon})\f$-interleaving, although in practice the * error is usually smaller. * A more intuitive presentation of the idea is available in * \cite cavanna15geometric, and in a video \cite cavanna15visualizing. * * The interface of `Sparse_rips_complex` is similar to the one for the usual - * `Rips_complex`, except that one has to specify the approximation factor, and - * there is no option to limit the maximum filtration value (the way the - * approximation is done means that larger filtration values are much cheaper - * to handle than low filtration values, so the gain would be too small). + * `Rips_complex`, except that one has to specify the approximation factor. + * There is an option to limit the minimum and maximum filtration values, but + * they are not recommended: the way the approximation is done means that + * larger filtration values are much cheaper to handle than low filtration + * values, so the gain in ignoring the large ones is small, and + * `Gudhi::subsampling::sparsify_point_set()` is a more efficient way of + * ignoring small filtration values. * * Theoretical guarantees are only available for \f$\epsilon<1\f$. The * construction accepts larger values of ε, and the size of the complex * keeps decreasing, but there is no guarantee on the quality of the result. + * Note that while the number of edges decreases when ε increases, the + * number of higher-dimensional simplices may not be monotonous when + * \f$\frac12\leq\epsilon\leq 1\f$. * * \section ripspointsdistance Point cloud and distance function * @@ -119,8 +125,8 @@ namespace rips_complex { * * \include Rips_complex/example_one_skeleton_rips_from_points.cpp * - * When launching (Rips maximal distance between 2 points is 12.0, is expanded until dimension 1 - one skeleton graph - * in other words): + * When launching (Rips maximal distance between 2 points is 12.0, is expanded + * until dimension 1 - one skeleton graph in other words): * * \code $> ./Rips_complex_example_one_skeleton_from_points * \endcode diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h index e902e52c..ee100867 100644 --- a/src/Rips_complex/include/gudhi/Rips_complex.h +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -90,7 +90,7 @@ class Rips_complex { * @param[in] threshold Rips value. * * \tparam DistanceMatrix must have a `size()` method and on which `distance_matrix[i][j]` returns - * the distance between points \f$i\f$ and \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * the distance between points \f$i\f$ and \f$j\f$ as long as \f$ 0 \leqslant j < i \leqslant * distance\_matrix.size().\f$ */ template<typename DistanceMatrix> diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 00da148f..8df6e387 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -47,7 +47,9 @@ namespace rips_complex { * * \details * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered - * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * simplicial complex that is multiplicatively + * \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. More precisely, + * this is a \f$(1,\frac{1}{1-\epsilon}\f$-interleaving. * * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. */ @@ -68,32 +70,36 @@ class Sparse_rips_complex { * @param[in] points Range of points. * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. * */ template <typename RandomAccessPointRange, typename Distance> - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) + : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); - std::vector<Vertex_handle> sorted_points; - std::vector<Filtration_value> params; auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; Ker<decltype(dist_fun)> kernel(dist_fun); subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, std::back_inserter(sorted_points), std::back_inserter(params)); - compute_sparse_graph(sorted_points, params, dist_fun, epsilon); + compute_sparse_graph(dist_fun, epsilon, mini, maxi); } /** \brief Sparse_rips_complex constructor from a distance matrix. * * @param[in] distance_matrix Range of range of distances. * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and - * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * \f$j\f$ as long as \f$ 0 \leqslant j < i \leqslant * distance\_matrix.size().\f$ * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. */ template <typename DistanceMatrix> - Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), - [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + [&](Vertex_handle i, Vertex_handle j) { return (i==j) ? 0 : (i<j) ? distance_matrix[j][i] : distance_matrix[i][j]; }, + epsilon, mini, maxi) {} /** \brief Fills the simplicial complex with the sparse Rips graph and * expands it with all the cliques, stopping at a given maximal dimension. @@ -111,7 +117,26 @@ class Sparse_rips_complex { std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); complex.insert_graph(graph_); - complex.expansion(dim_max); + if(epsilon_ >= 1) { + complex.expansion(dim_max); + return; + } + const int n = boost::size(params); + std::vector<Filtration_value> lambda(n); + // lambda[original_order]=params[sorted_order] + for(int i=0;i<n;++i) + lambda[sorted_points[i]] = params[i]; + double cst = epsilon_ * (1 - epsilon_) / 2; + auto block = [cst,&complex,&lambda](typename SimplicialComplexForRips::Simplex_handle sh){ + auto filt = complex.filtration(sh); + auto mini = filt * cst; + for(auto v : complex.simplex_vertex_range(sh)){ + if(lambda[v] < mini) + return true; // v died before this simplex could be born + } + return false; + }; + complex.expansion_with_blockers(dim_max, block); } private: @@ -127,9 +152,11 @@ class Sparse_rips_complex { }; // PointRange must be random access. - template <typename PointRange, typename ParamRange, typename Distance> - void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + template <typename Distance> + void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { + const auto& points = sorted_points; // convenience alias const int n = boost::size(points); + double cst = epsilon * (1 - epsilon) / 2; graph_.~Graph(); new (&graph_) Graph(n); // for(auto v : vertices(g)) // doesn't work :-( @@ -143,29 +170,43 @@ class Sparse_rips_complex { // TODO(MG): // - make it parallel // - only test near-enough neighbors - for (int i = 0; i < n; ++i) + for (int i = 0; i < n; ++i) { + auto&& pi = points[i]; + auto li = params[i]; + if (li < mini) break; for (int j = i + 1; j < n; ++j) { - auto&& pi = points[i]; auto&& pj = points[j]; auto d = dist(pi, pj); - auto li = params[i]; auto lj = params[j]; + if (lj < mini) break; GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); Filtration_value alpha; // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips if (d * epsilon <= 2 * lj) alpha = d; - else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) - alpha = (d - lj / epsilon) * 2; - else + else if (d * epsilon > li + lj) continue; + else { + alpha = (d - lj / epsilon) * 2; + // Keep the test exactly the same as in block to avoid inconsistencies + if (epsilon < 1 && alpha * cst > lj) + continue; + } - add_edge(pi, pj, alpha, graph_); + if (alpha <= maxi) + add_edge(pi, pj, alpha, graph_); } + } } Graph graph_; + double epsilon_; + // Because of the arbitrary split between constructor and create_complex + // sorted_points[sorted_order]=original_order + std::vector<Vertex_handle> sorted_points; + // params[sorted_order]=distance to previous points + std::vector<Filtration_value> params; }; } // namespace rips_complex diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 108cdd50..03838085 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -85,7 +85,7 @@ properly, this is a known issue. ## sparse_rips_persistence ## This program computes the persistent homology with coefficient field *Z/pZ* -of a sparse (1+epsilon)/(1-epsilon)-approximation of the Rips complex defined on a set of input Euclidean points. The output diagram contains one bar per line, written with the convention: +of a sparse 1/(1-epsilon)-approximation of the Rips complex defined on a set of input Euclidean points. The output diagram contains one bar per line, written with the convention: `p dim birth death` @@ -100,7 +100,7 @@ where `dim` is the dimension of the homological feature, `birth` and `death` are * `-h [ --help ]` Produce help message * `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. * `-e [ --approximation ]` (default = .5) Epsilon, where the sparse Rips complex is a (1+epsilon)/(1-epsilon)-approximation of the Rips complex. -* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute. +* `-d [ --cpx-dimension ]` (default = INT_MAX) Maximal dimension of the Rips complex we want to compute. * `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp index 6d4d86fd..3840d9f7 100644 --- a/src/Rips_complex/utilities/sparse_rips_persistence.cpp +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -98,7 +98,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: "Name of file in which the persistence diagram is written. Default print in std::cout")( "approximation,e", po::value<double>(&epsilon)->default_value(.5), "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")( - "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()), "Maximal dimension of the Rips complex we want to compute.")( "field-charac,p", po::value<int>(&p)->default_value(11), "Characteristic p of the coefficient field Z/pZ for computing homology.")( @@ -119,7 +119,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: if (vm.count("help") || !vm.count("input-file")) { std::cout << std::endl; std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; - std::cout << "of a sparse (1+epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; + std::cout << "of a sparse 1/(1-epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; std::cout << "The output diagram contains one bar per line, written with the convention: \n"; std::cout << " p dim b d \n"; std::cout << "where dim is the dimension of the homological feature,\n"; diff --git a/src/Tangential_complex/doc/Intro_tangential_complex.h b/src/Tangential_complex/doc/Intro_tangential_complex.h index f4fc8ac7..501f4a8b 100644 --- a/src/Tangential_complex/doc/Intro_tangential_complex.h +++ b/src/Tangential_complex/doc/Intro_tangential_complex.h @@ -35,9 +35,11 @@ namespace tangential_complex { \section tangentialdefinition Definition -A Tangential Delaunay complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> +A Tangential Delaunay complex is a +<a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> designed to reconstruct a \f$k\f$-dimensional smooth manifold embedded in \f$d\f$-dimensional Euclidean space. -The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of "small" intrinsic dimension. +The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of +"small" intrinsic dimension. The running time depends only linearly on the extrinsic dimension \f$ d \f$ and exponentially on the intrinsic dimension \f$ k \f$. @@ -46,17 +48,19 @@ An extensive description of the Tangential complex can be found in \cite tangent \subsection whatisthetc What is a Tangential Complex? Let us start with the description of the Tangential complex of a simple example, with \f$ k=1 \f$ and \f$ d=2 \f$. -The input data is 4 points \f$ P \f$ located on a curve embedded in 2D. +The point set \f$ \mathscr P \f$ is located on a closed curve embedded in 2D. +Only 4 points will be displayed (more are required for PCA) to simplify the figures. \image html "tc_example_01.png" "The input" -For each point \f$ p \f$, estimate its tangent subspace \f$ T_p \f$ (e.g. using PCA). +For each point \f$ P \f$, estimate its tangent subspace \f$ T_P \f$ using PCA. \image html "tc_example_02.png" "The estimated normals" -Let us add the Voronoi diagram of the points in orange. For each point \f$ p \f$, construct its star in the Delaunay triangulation of \f$ P \f$ restricted to \f$ T_p \f$. +Let us add the Voronoi diagram of the points in orange. For each point \f$ P \f$, construct its star in the Delaunay +triangulation of \f$ \mathscr P \f$ restricted to \f$ T_P \f$. \image html "tc_example_03.png" "The Voronoi diagram" The Tangential Delaunay complex is the union of those stars. In practice, neither the ambient Voronoi diagram nor the ambient Delaunay triangulation is computed. -Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only need the star of each point. -More details can be found in \cite tangentialcomplex2014. +Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only +need the star of each point. More details can be found in \cite tangentialcomplex2014. \subsection inconsistencies Inconsistencies @@ -65,7 +69,7 @@ An inconsistency occurs when a simplex is not in the star of all its vertices. Let us take the same example. \image html "tc_example_07_before.png" "Before" -Let us slightly move the tangent subspace \f$ T_q \f$ +Let us slightly move the tangent subspace \f$ T_Q \f$ \image html "tc_example_07_after.png" "After" Now, the star of \f$ Q \f$ contains \f$ QP \f$, but the star of \f$ P \f$ does not contain \f$ QP \f$. We have an inconsistency. \image html "tc_example_08.png" "After" diff --git a/src/Tangential_complex/example/example_basic.cpp b/src/Tangential_complex/example/example_basic.cpp index 4f2b859e..ab35edf0 100644 --- a/src/Tangential_complex/example/example_basic.cpp +++ b/src/Tangential_complex/example/example_basic.cpp @@ -1,5 +1,7 @@ #include <gudhi/Tangential_complex.h> #include <gudhi/sparsify_point_set.h> +//#include <gudhi/Fake_simplex_tree.h> + #include <CGAL/Epick_d.h> #include <CGAL/Random.h> @@ -20,7 +22,7 @@ CGAL::Parallel_tag> TC; int main(void) { const int INTRINSIC_DIM = 2; const int AMBIENT_DIM = 3; - const int NUM_POINTS = 1000; + const int NUM_POINTS = 100; Kernel k; @@ -37,6 +39,7 @@ 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/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index 37cdf1b4..4a78127c 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -322,7 +322,11 @@ class Tangential_complex { for (std::size_t i = 0; i < m_points.size(); ++i) m_are_tangent_spaces_computed[i] = true; } - /// Computes the tangential complex. + /** \brief Computes the tangential complex. + * \exception std::invalid_argument In debug mode, if the computed star dimension is too low. Try to set a bigger + * maximal edge length value with `Tangential_complex::set_max_squared_edge_length` if + * this happens. + */ void compute_tangential_complex() { #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS std::cerr << red << "WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. " @@ -1984,6 +1988,13 @@ class Tangential_complex { return os; } + /** \brief Sets the maximal possible squared edge length for the edges in the triangulations. + * + * @param[in] max_squared_edge_length Maximal possible squared edge length. + * + * If the maximal edge length value is too low `Tangential_complex::compute_tangential_complex` will throw an + * exception in debug mode. + */ void set_max_squared_edge_length(FT max_squared_edge_length) { m_max_squared_edge_length = max_squared_edge_length; } private: diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp index 4e2d4f65..103b8b30 100644 --- a/src/Tangential_complex/test/test_tangential_complex.cpp +++ b/src/Tangential_complex/test/test_tangential_complex.cpp @@ -126,3 +126,33 @@ BOOST_AUTO_TEST_CASE(test_mini_tangential) { BOOST_CHECK(stree.num_vertices() == 4); BOOST_CHECK(stree.num_simplices() == 6); } + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE(test_basic_example_throw) { + typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel; + typedef Kernel::FT FT; + typedef Kernel::Point_d Point; + typedef Kernel::Vector_d Vector; + typedef tc::Tangential_complex<Kernel, CGAL::Dynamic_dimension_tag,CGAL::Parallel_tag> TC; + + const int INTRINSIC_DIM = 2; + const int AMBIENT_DIM = 3; + const int NUM_POINTS = 1000; + + Kernel k; + + // Generate points on a 2-sphere + CGAL::Random_points_on_sphere_d<Point> generator(AMBIENT_DIM, 3.); + std::vector<Point> points; + points.reserve(NUM_POINTS); + for (int i = 0; i < NUM_POINTS; ++i) + points.push_back(*generator++); + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.set_max_squared_edge_length(0.01); + std::cout << "test_basic_example_throw - set_max_squared_edge_length(0.01) to make GUDHI_CHECK fail" << std::endl; + BOOST_CHECK_THROW(tc.compute_tangential_complex(), std::invalid_argument); + +} +#endif diff --git a/src/Toplex_map/benchmark/CMakeLists.txt b/src/Toplex_map/benchmark/CMakeLists.txt new file mode 100644 index 00000000..2d58a156 --- /dev/null +++ b/src/Toplex_map/benchmark/CMakeLists.txt @@ -0,0 +1,3 @@ +project(Toplex_map_benchmark) + +add_executable(Toplex_map_benchmark benchmark_tm.cpp) diff --git a/src/Toplex_map/benchmark/benchmark_tm.cpp b/src/Toplex_map/benchmark/benchmark_tm.cpp new file mode 100644 index 00000000..eedb442b --- /dev/null +++ b/src/Toplex_map/benchmark/benchmark_tm.cpp @@ -0,0 +1,137 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2018 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <random> +#include <chrono> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Toplex_map.h> +#include <gudhi/Lazy_toplex_map.h> + +using namespace Gudhi; + +typedef Toplex_map::Simplex Simplex; +typedef Toplex_map::Vertex Vertex; +typedef std::pair<Simplex_tree<>::Simplex_handle, bool> typePairSimplexBool; + +class ST_wrapper { + public: + void insert_simplex(const Simplex& tau) { + /*std::cout << "insert_simplex - " << simplexTree.num_simplices() << " - "; + for (auto v : tau) + std::cout << v << ", "; + std::cout << std::endl; + */ + simplexTree.insert_simplex_and_subfaces(tau); + } + + bool membership(const Simplex& tau) { return simplexTree.find(tau) != simplexTree.null_simplex(); } + + Vertex contraction(const Vertex x, const Vertex y) { + // TODO (VR): edge contraction is not yet available for Simplex_tree + return y; + } + + std::size_t num_maximal_simplices() { return simplexTree.num_simplices(); } + + private: + Simplex_tree<> simplexTree; + void erase_max(const Simplex& sigma) { + if (membership(sigma)) simplexTree.remove_maximal_simplex(simplexTree.find(sigma)); + } +}; + +int n = 300; + +int nb_insert_simplex1 = 3000; +int nb_membership1 = 4000; +int nb_contraction = 300; +int nb_insert_simplex2 = 3000; +int nb_membership2 = 400000; + +Simplex random_simplex(int n, std::size_t d) { + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<std::size_t> dis(1, n); + Simplex s; + while (s.size() < d) s.insert(dis(gen)); + return s; +} + +std::vector<Simplex> r_vector_simplices(int n, int max_d, int m) { + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<std::size_t> dis(1, max_d); + std::vector<Simplex> v; + for (int i = 0; i < m; i++) v.push_back(random_simplex(n, dis(gen))); + return v; +} + +template <typename complex_type> +void chrono(int n, int d) { + complex_type K; + std::vector<Simplex> simplices_insert_simplex1 = r_vector_simplices(n, d, nb_insert_simplex1); + std::vector<Simplex> simplices_membership1 = r_vector_simplices(n, d, nb_membership1); + std::vector<Simplex> simplices_insert_simplex2 = r_vector_simplices(n - 2 * nb_contraction, d, nb_insert_simplex2); + std::vector<Simplex> simplices_membership2 = r_vector_simplices(n - 2 * nb_contraction, d, nb_membership2); + std::chrono::time_point<std::chrono::system_clock> start, end; + + for (const Simplex& s : simplices_insert_simplex1) K.insert_simplex(s); + + for (const Simplex& s : simplices_membership1) K.membership(s); + + start = std::chrono::system_clock::now(); + for (int i = 1; i <= nb_contraction; i++) K.contraction(n - 2 * i, n - 2 * i - 1); + end = std::chrono::system_clock::now(); + auto c3 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + + start = std::chrono::system_clock::now(); + for (const Simplex& s : simplices_insert_simplex2) K.insert_simplex(s); + end = std::chrono::system_clock::now(); + auto c1 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + + start = std::chrono::system_clock::now(); + for (const Simplex& s : simplices_membership2) K.membership(s); + end = std::chrono::system_clock::now(); + auto c2 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); + + if (c3 > 0) + std::cout << c1 << "\t \t" << c2 << "\t \t" << c3 << "\t \t" << K.num_maximal_simplices() << std::endl; + else + std::cout << c1 << "\t \t" << c2 << "\t \tN/A\t \t" << K.num_maximal_simplices() << std::endl; +} + +int main() { + for (int d = 5; d <= 40; d += 5) { + std::cout << "d=" << d << " \t Insertions \t Membership \t Contractions \t Size" << std::endl; + std::cout << "T Map \t \t"; + chrono<Toplex_map>(n, d); + std::cout << "Lazy \t \t"; + chrono<Lazy_toplex_map>(n, d); + if (d <= 15) { + std::cout << "ST \t \t"; + chrono<ST_wrapper>(n, d); + } + std::cout << std::endl; + } +} diff --git a/src/Toplex_map/doc/Intro_Toplex_map.h b/src/Toplex_map/doc/Intro_Toplex_map.h new file mode 100644 index 00000000..a925dc2b --- /dev/null +++ b/src/Toplex_map/doc/Intro_Toplex_map.h @@ -0,0 +1,61 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2017 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_ +#define DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { + +/** \defgroup toplex_map Toplex Map + * + * \author François Godi + * @{ + * + * \section toplexmapdefinition Definition + * + * A Toplex_map is a data structure to represent and store a simplicial complex. A "toplex" is the contraction of + * "top-simplex", also known as a maximal simplex. The plural form of "toplex" will be called "toplices". + * + * Let's consider a simplicial complex, denote by \f$d\f$ its dimension and by \f$k\f$ its number of maximal simplices. + * Furthermore, denote by \f$\Gamma_0\f$ the maximal number of toplices, i.e. maximal simplices, that contain a same + * vertex. + * + * The goal of the Toplex Map is both to represent the complex in optimal O(kd) space and to provide fast standard + * operations such as : insertion, removal, contraction of an edge, collapses and membership of a simplex. The time + * needed for these operation is linear or quadratic in \f$\Gamma_0\f$ and \f$d\f$. + * + * Toplex map is composed firstly of a raw storage of toplices and secondly of a map which associate any vertex to a + * set of pointers toward all toplices containing this vertex. The data structure is described in + * \cite boissonnat_et_al:LIPIcs:2015:5098 (aka. Simplex Array List or SAL). + * + * \image html map.png + * + * The performances are a lot better than the `Simplex_tree` as soon you use maximal simplices and not simplices + * (cf. \cite DBLP:journals/corr/BoissonnatS16 ). + * + */ +/** @} */ // end defgroup toplex_map + +} // namespace Gudhi + +#endif // DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_ diff --git a/src/Toplex_map/doc/map.png b/src/Toplex_map/doc/map.png Binary files differnew file mode 100644 index 00000000..0f9cde2b --- /dev/null +++ b/src/Toplex_map/doc/map.png diff --git a/src/Toplex_map/example/CMakeLists.txt b/src/Toplex_map/example/CMakeLists.txt new file mode 100644 index 00000000..1d54fe87 --- /dev/null +++ b/src/Toplex_map/example/CMakeLists.txt @@ -0,0 +1,4 @@ +project(Toplex_map_examples) + +add_executable(Toplex_map_example_simple_toplex_map simple_toplex_map.cpp) +add_test(NAME Toplex_map_example_simple_toplex_map COMMAND $<TARGET_FILE:Toplex_map_example_simple_toplex_map>) diff --git a/src/Toplex_map/example/simple_toplex_map.cpp b/src/Toplex_map/example/simple_toplex_map.cpp new file mode 100644 index 00000000..e1c12ed6 --- /dev/null +++ b/src/Toplex_map/example/simple_toplex_map.cpp @@ -0,0 +1,165 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Toplex_map.h> + +#include <iostream> +#include <utility> // for pair +#include <vector> +#include <cassert> + +int main(int argc, char* const argv[]) { + using Simplex = Gudhi::Toplex_map::Simplex; + Simplex sigma1 = {1, 2, 3}; + Simplex sigma2 = {2, 3, 4, 5}; + + Gudhi::Toplex_map tm; + tm.insert_simplex(sigma1); + tm.insert_simplex(sigma2); + + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* /X\5/ */ + /* o---o */ + /* 1 3 */ + + std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices() + << std::endl; + + // Browse maximal cofaces + Simplex sigma3 = {2, 3}; + std::cout << "Maximal cofaces of {2, 3} are :" << std::endl; + for (auto simplex_ptr : tm.maximal_cofaces(sigma3, 2)) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + Simplex sigma4 = {1, 3}; + assert(tm.membership(sigma4)); + + Gudhi::Toplex_map::Vertex v = tm.contraction(1, 3); + std::cout << "After contraction(1, 3) - " << v << std::endl; + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* \5/ */ + /* o */ + /* 3 */ + std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices() + << std::endl; + + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + Simplex sigma5 = {3, 4}; + assert(tm.membership(sigma5)); + + v = tm.contraction(3, 4); + std::cout << "After contraction(3, 4) - " << v << std::endl; + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* \X/ */ + /* o */ + /* 5 */ + std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices() + << std::endl; + + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + tm.insert_simplex(sigma1); + tm.insert_simplex(sigma2); + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* /X\5/ */ + /* o---o */ + /* 1 3 */ + tm.remove_simplex(sigma1); + + std::cout << "After remove_simplex(1, 2, 3)" << std::endl; + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* / \5/ */ + /* o---o */ + /* 1 3 */ + std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices() + << std::endl; + + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + tm.remove_vertex(1); + + std::cout << "After remove_vertex(1)" << std::endl; + /* Simplex is: */ + /* 2 4 */ + /* o---o */ + /* \5/ */ + /* o */ + /* 3 */ + std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices() + << std::endl; + + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + } + + return 0; +} diff --git a/src/Toplex_map/include/gudhi/Lazy_toplex_map.h b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h new file mode 100644 index 00000000..b0b3706e --- /dev/null +++ b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h @@ -0,0 +1,271 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2018 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef LAZY_TOPLEX_MAP_H +#define LAZY_TOPLEX_MAP_H + +#include <gudhi/Toplex_map.h> +#include <boost/heap/fibonacci_heap.hpp> + +namespace Gudhi { + +/** + * \brief Lazy toplex map data structure for representing unfiltered simplicial complexes. + * + * \details A Toplex_map is an unordered map from vertices to maximal simplices (aka. toplices). + * The lazy version is not always up to date as it requires clean operation in order to be. + * + * \ingroup toplex_map */ +class Lazy_toplex_map { + public: + /** Vertex is the type of vertices. */ + using Vertex = Toplex_map::Vertex; + + /** Simplex is the type of simplices. */ + using Simplex = Toplex_map::Simplex; + + /** The type of the pointers to maximal simplices. */ + using Simplex_ptr = Toplex_map::Simplex_ptr; + + /** The type of the sets of Simplex_ptr. */ + using Simplex_ptr_set = Toplex_map::Simplex_ptr_set; + + /** Adds the given simplex to the complex. + * The simplex must not be in the complex already, and it must not contain one of the current toplices. */ + template <typename Input_vertex_range> + void insert_independent_simplex(const Input_vertex_range &vertex_range); + + /** \brief Adds the given simplex to the complex. + * Nothing happens if the simplex is already in the complex (i.e. it is a face of one of the toplices). */ + template <typename Input_vertex_range> + bool insert_simplex(const Input_vertex_range &vertex_range); + + /** \brief Removes the given simplex and its cofaces from the complex. + * Its faces are kept inside. */ + template <typename Input_vertex_range> + void remove_simplex(const Input_vertex_range &vertex_range); + + /** Does a simplex belong to the complex ? */ + template <typename Input_vertex_range> + bool membership(const Input_vertex_range &vertex_range); + + /** Do all the facets of a simplex belong to the complex ? */ + template <typename Input_vertex_range> + bool all_facets_inside(const Input_vertex_range &vertex_range); + + /** Contracts one edge in the complex. + * The edge has to verify the link condition if you want to preserve topology. + * Returns the remaining vertex. */ + Vertex contraction(const Vertex x, const Vertex y); + + /** \brief Number of maximal simplices. */ + std::size_t num_maximal_simplices() const { return size; } + + /** \brief Number of vertices. */ + std::size_t num_vertices() const { return t0.size(); } + + private: + template <typename Input_vertex_range> + void erase_max(const Input_vertex_range &vertex_range); + template <typename Input_vertex_range> + Vertex best_index(const Input_vertex_range &vertex_range); + void clean(const Vertex v); + + std::unordered_map<Vertex, std::size_t> gamma0_lbounds; + + std::unordered_map<Vertex, Simplex_ptr_set> t0; + bool empty_toplex; // Is the empty simplex a toplex ? + + typedef boost::heap::fibonacci_heap<std::pair<std::size_t, Vertex>> PriorityQueue; + PriorityQueue cleaning_priority; + std::unordered_map<Vertex, PriorityQueue::handle_type> cp_handles; + + std::size_t get_gamma0_lbound(const Vertex v) const; + + std::size_t size_lbound = 0; + std::size_t size = 0; + + const double ALPHA = 4; // time + const double BETTA = 8; // memory +}; + +template <typename Input_vertex_range> +void Lazy_toplex_map::insert_independent_simplex(const Input_vertex_range &vertex_range) { + for (const Vertex &v : vertex_range) + if (!gamma0_lbounds.count(v)) + gamma0_lbounds.emplace(v, 1); + else + gamma0_lbounds[v]++; + size_lbound++; + insert_simplex(vertex_range); +} + +template <typename Input_vertex_range> +bool Lazy_toplex_map::insert_simplex(const Input_vertex_range &vertex_range) { + Simplex sigma(vertex_range.begin(), vertex_range.end()); + // Check empty face management + empty_toplex = (sigma.size() == 0); + Simplex_ptr sptr = std::make_shared<Simplex>(sigma); + bool inserted = false; + for (const Vertex &v : sigma) { + if (!t0.count(v)) { + t0.emplace(v, Simplex_ptr_set()); + auto v_handle = cleaning_priority.push(std::make_pair(0, v)); + cp_handles.emplace(v, v_handle); + } + inserted = t0.at(v).emplace(sptr).second; + cleaning_priority.update(cp_handles.at(v), std::make_pair(t0.at(v).size() - get_gamma0_lbound(v), v)); + } + if (inserted) size++; + if (size > (size_lbound + 1) * BETTA) clean(cleaning_priority.top().second); + return inserted; +} + +template <typename Input_vertex_range> +void Lazy_toplex_map::remove_simplex(const Input_vertex_range &vertex_range) { + if (vertex_range.begin() == vertex_range.end()) { + t0.clear(); + gamma0_lbounds.clear(); + cleaning_priority.clear(); + size_lbound = 0; + size = 0; + empty_toplex = false; + } else { + const Vertex &v = best_index(vertex_range); + // Copy constructor needed because the set is modified + if (t0.count(v)) + for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(v))) + if (included(vertex_range, *sptr)) { + erase_max(*sptr); + for (const Simplex &f : facets(vertex_range)) insert_independent_simplex(f); + } + } +} + +template <typename Input_vertex_range> +bool Lazy_toplex_map::membership(const Input_vertex_range &vertex_range) { + if (t0.size() == 0 && !empty_toplex) return false; // empty complex + if (vertex_range.begin() == vertex_range.end()) return true; // empty query simplex + Vertex v = best_index(vertex_range); + if (!t0.count(v)) return false; + for (const Simplex_ptr &sptr : t0.at(v)) + if (included(vertex_range, *sptr)) return true; + return false; +} + +template <typename Input_vertex_range> +bool Lazy_toplex_map::all_facets_inside(const Input_vertex_range &vertex_range) { + Simplex sigma(vertex_range.begin(), vertex_range.end()); + Vertex v = best_index(sigma); + if (!t0.count(v)) return false; + Simplex f = sigma; + f.erase(v); + if (!membership(f)) return false; + std::unordered_set<Vertex> facets_inside; + for (const Simplex_ptr &sptr : t0.at(v)) + for (const Vertex &w : sigma) { + f = sigma; + f.erase(w); + if (included(f, *sptr)) facets_inside.insert(w); + } + return facets_inside.size() == sigma.size() - 1; +} + +/* Returns the remaining vertex */ +Lazy_toplex_map::Vertex Lazy_toplex_map::contraction(const Vertex x, const Vertex y) { + if (!t0.count(x)) return y; + if (!t0.count(y)) return x; + Vertex k, d; + if (t0.at(x).size() > t0.at(y).size()) + k = x, d = y; + else + k = y, d = x; + // Copy constructor needed because the set is modified + for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(d))) { + Simplex sigma(*sptr); + erase_max(sigma); + sigma.erase(d); + sigma.insert(k); + insert_simplex(sigma); + } + t0.erase(d); + return k; +} + +/* No facets insert_simplexed */ +template <typename Input_vertex_range> +inline void Lazy_toplex_map::erase_max(const Input_vertex_range &vertex_range) { + Simplex sigma(vertex_range.begin(), vertex_range.end()); + empty_toplex = false; + Simplex_ptr sptr = std::make_shared<Simplex>(sigma); + bool erased = false; + for (const Vertex &v : sigma) { + erased = t0.at(v).erase(sptr) > 0; + if (t0.at(v).size() == 0) t0.erase(v); + } + if (erased) size--; +} + +template <typename Input_vertex_range> +Lazy_toplex_map::Vertex Lazy_toplex_map::best_index(const Input_vertex_range &vertex_range) { + Simplex tau(vertex_range.begin(), vertex_range.end()); + std::size_t min = std::numeric_limits<size_t>::max(); + Vertex arg_min = -1; + for (const Vertex &v : tau) + if (!t0.count(v)) + return v; + else if (t0.at(v).size() < min) + min = t0.at(v).size(), arg_min = v; + if (min > ALPHA * get_gamma0_lbound(arg_min)) clean(arg_min); + return arg_min; +} + +std::size_t Lazy_toplex_map::get_gamma0_lbound(const Vertex v) const { + return gamma0_lbounds.count(v) ? gamma0_lbounds.at(v) : 0; +} + +void Lazy_toplex_map::clean(const Vertex v) { + Toplex_map toplices; + std::unordered_map<int, std::vector<Simplex>> dsorted_simplices; + std::size_t max_dim = 0; + for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(v))) { + if (sptr->size() > max_dim) { + for (std::size_t d = max_dim + 1; d <= sptr->size(); d++) dsorted_simplices.emplace(d, std::vector<Simplex>()); + max_dim = sptr->size(); + } + dsorted_simplices[sptr->size()].emplace_back(*sptr); + erase_max(*sptr); + } + for (std::size_t d = max_dim; d >= 1; d--) + for (const Simplex &s : dsorted_simplices.at(d)) + if (!toplices.membership(s)) toplices.insert_independent_simplex(s); + Simplex sv; + sv.insert(v); + auto clean_cofaces = toplices.maximal_cofaces(sv); + size_lbound = size_lbound - get_gamma0_lbound(v) + clean_cofaces.size(); + gamma0_lbounds[v] = clean_cofaces.size(); + for (const Simplex_ptr &sptr : clean_cofaces) insert_simplex(*sptr); +} + +} // namespace Gudhi + +#endif /* LAZY_TOPLEX_MAP_H */ diff --git a/src/Toplex_map/include/gudhi/Toplex_map.h b/src/Toplex_map/include/gudhi/Toplex_map.h new file mode 100644 index 00000000..4dc2331f --- /dev/null +++ b/src/Toplex_map/include/gudhi/Toplex_map.h @@ -0,0 +1,338 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2018 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef TOPLEX_MAP_H +#define TOPLEX_MAP_H + +#include <vector> +#include <set> +#include <unordered_set> +#include <unordered_map> +#include <memory> +#include <limits> + +namespace Gudhi { + +/** + * \brief Toplex map data structure for representing unfiltered simplicial complexes. + * + * \details A Toplex_map is an unordered map from vertices to maximal simplices (aka. toplices). + * + * \ingroup toplex_map */ +class Toplex_map { + public: + /** Vertex is the type of vertices. */ + using Vertex = std::size_t; + + /** Simplex is the type of simplices. */ + using Simplex = std::set<Vertex>; + + /** The type of the pointers to maximal simplices. */ + using Simplex_ptr = std::shared_ptr<Toplex_map::Simplex>; + + struct Sptr_hash { + std::size_t operator()(const Toplex_map::Simplex_ptr& s) const; + }; + + struct Sptr_equal { + std::size_t operator()(const Toplex_map::Simplex_ptr& a, const Toplex_map::Simplex_ptr& b) const; + }; + + /** The type of the sets of Toplex_map::Simplex_ptr. */ + using Simplex_ptr_set = std::unordered_set<Toplex_map::Simplex_ptr, Sptr_hash, Sptr_equal>; + + /** \brief Adds the given simplex to the complex. + * Nothing happens if the simplex is already in the complex (i.e. it is a face of one of the toplices). */ + template <typename Input_vertex_range> + void insert_simplex(const Input_vertex_range& vertex_range); + + /** \brief Removes the given simplex and its cofaces from the complex. + * Its faces are kept inside. */ + template <typename Input_vertex_range> + void remove_simplex(const Input_vertex_range& vertex_range); + + /** Does a simplex belong to the complex ? */ + template <typename Input_vertex_range> + bool membership(const Input_vertex_range& vertex_range) const; + + /** Does a simplex is a toplex ? */ + template <typename Input_vertex_range> + bool maximality(const Input_vertex_range& vertex_range) const; + + /** Gives a set of pointers to the maximal cofaces of a simplex. + * Gives all the toplices if given the empty simplex. + * Gives not more than max_number maximal cofaces if max_number is strictly positive. */ + template <typename Input_vertex_range> + Toplex_map::Simplex_ptr_set maximal_cofaces(const Input_vertex_range& vertex_range, + const std::size_t max_number = 0) const; + + /** Gives a set of pointers to the maximal simplices. + * Gives not more than max_number maximal cofaces if max_number is strictly positive. */ + Toplex_map::Simplex_ptr_set maximal_simplices(const std::size_t max_number = 0) const { + return maximal_cofaces(Simplex(), max_number); + } + + /** Contracts one edge in the complex. + * The edge has to verify the link condition if you want to preserve topology. + * Returns the remaining vertex. */ + Vertex contraction(const Vertex x, const Vertex y); + + /** Remove the vertex and all its cofaces from the complex. */ + void remove_vertex(const Vertex x); + + /** \brief Number of maximal simplices. */ + std::size_t num_maximal_simplices() const { return maximal_simplices().size(); } + + /** \brief Number of vertices. */ + std::size_t num_vertices() const { return t0.size(); } + + std::set<Vertex> unitary_collapse(const Vertex k, const Vertex d); + + /** Adds the given simplex to the complex. + * The simplex must not be in the complex already, and it must not contain one of the current toplices. */ + template <typename Input_vertex_range> + void insert_independent_simplex(const Input_vertex_range& vertex_range); + + protected: + /** \internal Gives an index in order to look for a simplex quickly. */ + template <typename Input_vertex_range> + Vertex best_index(const Input_vertex_range& vertex_range) const; + + /** \internal The map from vertices to toplices */ + std::unordered_map<Vertex, Toplex_map::Simplex_ptr_set> t0; + + const Vertex VERTEX_UPPER_BOUND = std::numeric_limits<Vertex>::max(); + + /** \internal Removes a toplex without adding facets after. */ + void erase_maximal(const Toplex_map::Simplex_ptr& sptr); +}; + +// Pointers are also used as key in the hash sets. +template <typename Input_vertex_range> +Toplex_map::Simplex_ptr get_key(const Input_vertex_range& vertex_range); + +// Is the first simplex a face of the second ? +template <typename Input_vertex_range1, typename Input_vertex_range2> +bool included(const Input_vertex_range1& vertex_range1, const Input_vertex_range2& vertex_range2); + +// All the facets of the given simplex. +template <typename Input_vertex_range> +std::vector<Toplex_map::Simplex> facets(const Input_vertex_range& vertex_range); + +template <typename Input_vertex_range> +void Toplex_map::insert_simplex(const Input_vertex_range& vertex_range) { + if (membership(vertex_range)) return; + bool replace_facets = true; + for (const Toplex_map::Simplex& facet : facets(vertex_range)) + if (!maximality(facet)) { + replace_facets = false; + break; + } + if (replace_facets) + for (const Toplex_map::Simplex& facet : facets(vertex_range)) erase_maximal(get_key(facet)); + else + for (const Vertex& v : vertex_range) + if (t0.count(v)) + for (const Toplex_map::Simplex_ptr& fptr : Simplex_ptr_set(t0.at(v))) + // Copy constructor needed because the set is modified + if (included(*fptr, vertex_range)) erase_maximal(fptr); + // We erase all the maximal faces of the simplex + insert_independent_simplex(vertex_range); +} + +template <typename Input_vertex_range> +void Toplex_map::remove_simplex(const Input_vertex_range& vertex_range) { + if (vertex_range.begin() == vertex_range.end()) t0.clear(); + // Removal of the empty simplex means cleaning everything + else { + const Vertex& v = best_index(vertex_range); + if (t0.count(v)) + for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(v))) + // Copy constructor needed because the set is modified + if (included(vertex_range, *sptr)) { + erase_maximal(sptr); + for (const Toplex_map::Simplex& f : facets(vertex_range)) + if (!membership(f)) insert_independent_simplex(f); + // We add the facets which are new maximal simplices + } + } +} + +template <typename Input_vertex_range> +bool Toplex_map::membership(const Input_vertex_range& vertex_range) const { + if (t0.size() == 0) return false; + const Vertex& v = best_index(vertex_range); + if (!t0.count(v)) return false; + if (maximality(vertex_range)) return true; + for (const Toplex_map::Simplex_ptr& sptr : t0.at(v)) + if (included(vertex_range, *sptr)) return true; + return false; +} + +template <typename Input_vertex_range> +bool Toplex_map::maximality(const Input_vertex_range& vertex_range) const { + const Vertex& v = best_index(vertex_range); + if (!t0.count(v)) return false; + return t0.at(v).count(get_key(vertex_range)); +} + +template <typename Input_vertex_range> +Toplex_map::Simplex_ptr_set Toplex_map::maximal_cofaces(const Input_vertex_range& vertex_range, + const std::size_t max_number) const { + Simplex_ptr_set cofaces; + if (maximality(vertex_range)) + cofaces.emplace(get_key(vertex_range)); + else if (vertex_range.begin() == vertex_range.end()) + for (const auto& kv : t0) + for (const Toplex_map::Simplex_ptr& sptr : kv.second) { + // kv.second is a Simplex_ptr_set + cofaces.emplace(sptr); + if (cofaces.size() == max_number) return cofaces; + } + else { + const Vertex& v = best_index(vertex_range); + if (t0.count(v)) + for (const Toplex_map::Simplex_ptr& sptr : t0.at(v)) + if (included(vertex_range, *sptr)) { + cofaces.emplace(sptr); + if (cofaces.size() == max_number) return cofaces; + } + } + return cofaces; +} + +Toplex_map::Vertex Toplex_map::contraction(const Toplex_map::Vertex x, const Toplex_map::Vertex y) { + if (!t0.count(x)) return y; + if (!t0.count(y)) return x; + int k, d; + if (t0.at(x).size() > t0.at(y).size()) + k = x, d = y; + else + k = y, d = x; + for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))) { + // Copy constructor needed because the set is modified + Simplex sigma(*sptr); + erase_maximal(sptr); + sigma.erase(d); + sigma.insert(k); + insert_simplex(sigma); + } + return k; +} + +std::set<Toplex_map::Vertex> Toplex_map::unitary_collapse(const Toplex_map::Vertex k, const Toplex_map::Vertex d) { + Toplex_map::Simplex r; + for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))) { + // Copy constructor needed because the set is modified + Simplex sigma(*sptr); + erase_maximal(sptr); + sigma.erase(d); + for (const Vertex v : sigma) r.insert(v); + sigma.insert(k); + insert_simplex(sigma); + } + return r; +} + +template <typename Input_vertex_range> +void Toplex_map::insert_independent_simplex(const Input_vertex_range& vertex_range) { + auto key = get_key(vertex_range); + for (const Vertex& v : vertex_range) { + if (!t0.count(v)) t0.emplace(v, Simplex_ptr_set()); + t0.at(v).emplace(key); + } +} + +void Toplex_map::remove_vertex(const Toplex_map::Vertex x) { + for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(x))) { + Simplex sigma(*sptr); + erase_maximal(sptr); + sigma.erase(x); + insert_simplex(sigma); + } +} + +inline void Toplex_map::erase_maximal(const Toplex_map::Simplex_ptr& sptr) { + Simplex sigma(*sptr); + if (sptr->size() == 0) sigma.insert(VERTEX_UPPER_BOUND); + for (const Vertex& v : sigma) { + t0.at(v).erase(sptr); + if (t0.at(v).size() == 0) t0.erase(v); + } +} + +template <typename Input_vertex_range> +Toplex_map::Vertex Toplex_map::best_index(const Input_vertex_range& vertex_range) const { + std::size_t min = std::numeric_limits<size_t>::max(); + Vertex arg_min = VERTEX_UPPER_BOUND; + for (const Vertex& v : vertex_range) + if (!t0.count(v)) + return v; + else if (t0.at(v).size() < min) + min = t0.at(v).size(), arg_min = v; + return arg_min; +} + +std::size_t Toplex_map::Sptr_equal::operator()(const Toplex_map::Simplex_ptr& s1, + const Toplex_map::Simplex_ptr& s2) const { + if (s1->size() != s2->size()) return false; + return included(*s1, *s2); + // inclusion tests equality for same size simplices +} + +std::size_t Toplex_map::Sptr_hash::operator()(const Toplex_map::Simplex_ptr& s) const { + std::hash<double> h_f; + // double hash works better than int hash + size_t h = 0; + for (const Vertex& v : *s) h += h_f(static_cast<double>(v)); + return h; +} + +template <typename Input_vertex_range> +Toplex_map::Simplex_ptr get_key(const Input_vertex_range& vertex_range) { + Toplex_map::Simplex s(vertex_range.begin(), vertex_range.end()); + return std::make_shared<Toplex_map::Simplex>(s); +} + +template <typename Input_vertex_range1, typename Input_vertex_range2> +bool included(const Input_vertex_range1& vertex_range1, const Input_vertex_range2& vertex_range2) { + Toplex_map::Simplex s2(vertex_range2.begin(), vertex_range2.end()); + for (const Toplex_map::Vertex& v : vertex_range1) + if (!s2.count(v)) return false; + return true; +} + +template <typename Input_vertex_range> +std::vector<Toplex_map::Simplex> facets(const Input_vertex_range& vertex_range) { + std::vector<Toplex_map::Simplex> facets; + Toplex_map::Simplex f(vertex_range.begin(), vertex_range.end()); + for (const Toplex_map::Vertex& v : vertex_range) { + f.erase(v); + facets.emplace_back(f); + f.insert(v); + } + return facets; +} + +} // namespace Gudhi + +#endif /* TOPLEX_MAP_H */ diff --git a/src/Toplex_map/test/CMakeLists.txt b/src/Toplex_map/test/CMakeLists.txt new file mode 100644 index 00000000..59517db5 --- /dev/null +++ b/src/Toplex_map/test/CMakeLists.txt @@ -0,0 +1,11 @@ +project(Toplex_map_tests) + +include(GUDHI_test_coverage) + +add_executable( Toplex_map_unit_test toplex_map_unit_test.cpp ) +target_link_libraries(Toplex_map_unit_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +gudhi_add_coverage_test(Toplex_map_unit_test) + +add_executable( Lazy_toplex_map_unit_test lazy_toplex_map_unit_test.cpp ) +target_link_libraries(Lazy_toplex_map_unit_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +gudhi_add_coverage_test(Lazy_toplex_map_unit_test) diff --git a/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp b/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp new file mode 100644 index 00000000..a050cc92 --- /dev/null +++ b/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp @@ -0,0 +1,170 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2018 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <vector> +#include <gudhi/Lazy_toplex_map.h> + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "lazy toplex map" +#include <boost/test/unit_test.hpp> + +BOOST_AUTO_TEST_CASE(toplex_map) { + using Vertex = Gudhi::Lazy_toplex_map::Vertex; + + Gudhi::Lazy_toplex_map tm; + std::cout << "insert_simplex {1, 2, 3, 4}" << std::endl; + std::vector<Vertex> sigma1 = {1, 2, 3, 4}; + tm.insert_simplex(sigma1); + std::cout << "insert_simplex {5, 2, 3, 6}" << std::endl; + std::vector<Vertex> sigma2 = {5, 2, 3, 6}; + tm.insert_simplex(sigma2); + std::cout << "insert_simplex {5}" << std::endl; + std::vector<Vertex> sigma3 = {5}; + tm.insert_simplex(sigma3); + std::cout << "insert_simplex {4, 5, 3}" << std::endl; + std::vector<Vertex> sigma6 = {4, 5, 3}; + tm.insert_simplex(sigma6); + std::cout << "insert_simplex {4, 5, 9}" << std::endl; + std::vector<Vertex> sigma7 = {4, 5, 9}; + tm.insert_simplex(sigma7); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 5); + + std::vector<Vertex> sigma4 = {5, 2, 3}; + std::vector<Vertex> sigma5 = {5, 2, 7}; + BOOST_CHECK(tm.membership(sigma4)); + BOOST_CHECK(!tm.membership(sigma5)); + std::cout << "insert_simplex {5, 2, 7}" << std::endl; + tm.insert_simplex(sigma5); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 6); + + BOOST_CHECK(tm.membership(sigma5)); + + std::cout << "contraction(4,5)" << std::endl; + auto r = tm.contraction(4, 5); + std::cout << "r=" << r << std::endl; + BOOST_CHECK(r == 5); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 6); + + std::vector<Vertex> sigma8 = {1, 2, 3}; + std::vector<Vertex> sigma9 = {2, 7}; + + sigma8.emplace_back(r); + sigma9.emplace_back(r); + BOOST_CHECK(!tm.membership(sigma6)); + BOOST_CHECK(tm.membership(sigma8)); + BOOST_CHECK(tm.membership(sigma9)); + + std::cout << "remove_simplex({2, 7, r = 5})" << std::endl; + tm.remove_simplex(sigma9); + BOOST_CHECK(!tm.membership(sigma9)); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 8); + + // {2, 7, 5} is removed, but verify its edges are still there + std::vector<Vertex> edge = {2, 7}; + BOOST_CHECK(tm.membership(edge)); + edge = {2, 5}; + BOOST_CHECK(tm.membership(edge)); + edge = {7, 5}; + BOOST_CHECK(tm.membership(edge)); +} + +BOOST_AUTO_TEST_CASE(toplex_map_empty_toplex) { + using Vertex = Gudhi::Lazy_toplex_map::Vertex; + + Gudhi::Lazy_toplex_map tm; + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 0); + std::cout << "num_vertices = " << tm.num_vertices() << std::endl; + BOOST_CHECK(tm.num_vertices() == 0); + + std::cout << "Check an empty simplex is a member." << std::endl; + std::vector<Vertex> empty_sigma = {}; + BOOST_CHECK(tm.membership(empty_sigma)); + + std::cout << "Check the edge 2,7 is not a member." << std::endl; + std::vector<Vertex> edge = {2, 7}; + BOOST_CHECK(!tm.membership(edge)); + + std::cout << "Insert an empty simplex." << std::endl; + tm.insert_simplex(empty_sigma); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 0); + std::cout << "num_vertices = " << tm.num_vertices() << std::endl; + BOOST_CHECK(tm.num_vertices() == 0); + + std::cout << "Check an empty simplex is a member." << std::endl; + BOOST_CHECK(tm.membership(empty_sigma)); + std::cout << "Check the edge 2,7 is not a member." << std::endl; + BOOST_CHECK(!tm.membership(edge)); + + std::cout << "Insert edge 2,7." << std::endl; + tm.insert_simplex(edge); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 1); + std::cout << "num_vertices = " << tm.num_vertices() << std::endl; + BOOST_CHECK(tm.num_vertices() == 2); + + std::cout << "Check an empty simplex is a member." << std::endl; + BOOST_CHECK(tm.membership(empty_sigma)); + std::cout << "Check the edge 2,7 is a member." << std::endl; + BOOST_CHECK(tm.membership(edge)); + + std::cout << "contraction(2,7)" << std::endl; + auto r = tm.contraction(2, 7); + std::cout << "r=" << r << std::endl; + BOOST_CHECK(r == 7); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 1); + std::cout << "num_vertices = " << tm.num_vertices() << std::endl; + BOOST_CHECK(tm.num_vertices() == 1); + + std::cout << "Check an empty simplex is a member." << std::endl; + BOOST_CHECK(tm.membership(empty_sigma)); + std::cout << "Check the edge 2,7 is not a member." << std::endl; + BOOST_CHECK(!tm.membership(edge)); + + std::cout << "Remove the vertex 7." << std::endl; + std::vector<Vertex> vertex = {7}; + tm.remove_simplex(vertex); + + std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 0); + std::cout << "num_vertices = " << tm.num_vertices() << std::endl; + BOOST_CHECK(tm.num_vertices() == 0); + + std::cout << "Check an empty simplex is a member." << std::endl; + BOOST_CHECK(tm.membership(empty_sigma)); + std::cout << "Check the edge 2,7 is not a member." << std::endl; + BOOST_CHECK(!tm.membership(edge)); +} diff --git a/src/Toplex_map/test/toplex_map_unit_test.cpp b/src/Toplex_map/test/toplex_map_unit_test.cpp new file mode 100644 index 00000000..2bd27936 --- /dev/null +++ b/src/Toplex_map/test/toplex_map_unit_test.cpp @@ -0,0 +1,138 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: François Godi, Vincent Rouvreau + * + * Copyright (C) 2018 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <iostream> +#include <vector> +#include <gudhi/Toplex_map.h> + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "toplex map" +#include <boost/test/unit_test.hpp> + +BOOST_AUTO_TEST_CASE(toplex_map) { + using Vertex = Gudhi::Toplex_map::Vertex; + + Gudhi::Toplex_map tm; + std::cout << "insert_simplex {1, 2, 3, 4}" << std::endl; + std::vector<Vertex> sigma1 = {1, 2, 3, 4}; + tm.insert_simplex(sigma1); + std::cout << "insert_simplex {5, 2, 3, 6}" << std::endl; + std::vector<Vertex> sigma2 = {5, 2, 3, 6}; + tm.insert_simplex(sigma2); + std::cout << "insert_simplex {5}" << std::endl; + std::vector<Vertex> sigma3 = {5}; + tm.insert_simplex(sigma3); + std::cout << "insert_simplex {4, 5, 3}" << std::endl; + std::vector<Vertex> sigma6 = {4, 5, 3}; + tm.insert_simplex(sigma6); + std::cout << "insert_simplex {4, 5, 9}" << std::endl; + std::vector<Vertex> sigma7 = {4, 5, 9}; + tm.insert_simplex(sigma7); + + std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 4); + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + BOOST_CHECK(tm.maximality(*simplex_ptr)); + } + + BOOST_CHECK(tm.maximality(sigma1)); + BOOST_CHECK(tm.maximality(sigma2)); + BOOST_CHECK(!tm.maximality(sigma3)); + BOOST_CHECK(tm.maximality(sigma6)); + BOOST_CHECK(tm.maximality(sigma7)); + + std::vector<Vertex> sigma4 = {5, 2, 3}; + std::vector<Vertex> sigma5 = {5, 2, 7}; + BOOST_CHECK(tm.membership(sigma4)); + BOOST_CHECK(!tm.membership(sigma5)); + std::cout << "insert_simplex {5, 2, 7}" << std::endl; + tm.insert_simplex(sigma5); + + std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 5); + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + BOOST_CHECK(tm.maximality(*simplex_ptr)); + } + + BOOST_CHECK(tm.membership(sigma5)); + + std::cout << "contraction(4,5)" << std::endl; + auto r = tm.contraction(4, 5); + std::cout << "r=" << r << std::endl; + BOOST_CHECK(r == 5); + + std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 4); + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + BOOST_CHECK(tm.maximality(*simplex_ptr)); + } + + std::vector<Vertex> sigma8 = {1, 2, 3}; + std::vector<Vertex> sigma9 = {2, 7}; + + sigma8.emplace_back(r); + sigma9.emplace_back(r); + BOOST_CHECK(!tm.membership(sigma6)); + BOOST_CHECK(tm.membership(sigma8)); + BOOST_CHECK(tm.membership(sigma9)); + + std::cout << "remove_simplex({2, 7, r = 5})" << std::endl; + tm.remove_simplex(sigma9); + BOOST_CHECK(!tm.membership(sigma9)); + + std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl; + BOOST_CHECK(tm.num_maximal_simplices() == 5); + // Browse maximal simplices + std::cout << "Maximal simplices are :" << std::endl; + for (auto simplex_ptr : tm.maximal_simplices()) { + for (auto v : *simplex_ptr) { + std::cout << v << ", "; + } + std::cout << std::endl; + BOOST_CHECK(tm.maximality(*simplex_ptr)); + } + // {2, 7, 5} is removed, but verify its edges are still there + std::vector<Vertex> edge = {2, 7}; + BOOST_CHECK(tm.membership(edge)); + edge = {2, 5}; + BOOST_CHECK(tm.membership(edge)); + edge = {7, 5}; + BOOST_CHECK(tm.membership(edge)); +} diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index 7b0d350d..57ea7d14 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -148,6 +148,7 @@ if( PYTHONINTERP_FOUND ) find_python_module("matplotlib") find_python_module("numpy") find_python_module("scipy") + find_python_module("sphinx") endif() if(NOT GUDHI_CYTHON_PATH) @@ -157,20 +158,16 @@ endif(NOT GUDHI_CYTHON_PATH) option(WITH_GUDHI_CYTHON_RUNTIME_LIBRARY_DIRS "Build with setting runtime_library_dirs. Usefull when setting rpath is not allowed" ON) if(PYTHONINTERP_FOUND AND CYTHON_FOUND) - # Default found version 2 - if(PYTHON_VERSION_MAJOR EQUAL 2) + if(SPHINX_FOUND) # Documentation generation is available through sphinx find_program( SPHINX_PATH sphinx-build ) - elseif(PYTHON_VERSION_MAJOR EQUAL 3) - # No sphinx-build in Pyton3, just hack it - set(SPHINX_PATH "${PYTHON_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build.py") - else() - message(FATAL_ERROR "ERROR: Try to compile the Cython interface. Python version ${PYTHON_VERSION_STRING} is not valid.") - endif(PYTHON_VERSION_MAJOR EQUAL 2) - # get PYTHON_SITE_PACKAGES relative path from a python command line - execute_process( - COMMAND "${PYTHON_EXECUTABLE}" -c "from distutils.sysconfig import get_python_lib; print (get_python_lib(prefix='', plat_specific=True))" - OUTPUT_VARIABLE PYTHON_SITE_PACKAGES - OUTPUT_STRIP_TRAILING_WHITESPACE) + + if(NOT SPHINX_PATH) + if(PYTHON_VERSION_MAJOR EQUAL 3) + # In Python3, just hack sphinx-build if it does not exist + set(SPHINX_PATH "${PYTHON_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build.py") + endif(PYTHON_VERSION_MAJOR EQUAL 3) + endif(NOT SPHINX_PATH) + endif(SPHINX_FOUND) endif(PYTHONINTERP_FOUND AND CYTHON_FOUND) diff --git a/src/cmake/modules/GUDHI_user_version_target.cmake b/src/cmake/modules/GUDHI_user_version_target.cmake index d43a6fa6..2ed48c48 100644 --- a/src/cmake/modules/GUDHI_user_version_target.cmake +++ b/src/cmake/modules/GUDHI_user_version_target.cmake @@ -27,7 +27,7 @@ add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_SOURCE_DIR}/Conventions.txt ${GUDHI_USER_VERSION_DIR}/Conventions.txt) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E - copy ${CMAKE_SOURCE_DIR}/README ${GUDHI_USER_VERSION_DIR}/README) + copy ${CMAKE_SOURCE_DIR}/README.md ${GUDHI_USER_VERSION_DIR}/README.md) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_SOURCE_DIR}/COPYING ${GUDHI_USER_VERSION_DIR}/COPYING) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index b33a95a1..afe9b68c 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -171,6 +171,23 @@ </td> </tr> </table> + \subsection ToplexMapDataStructure Toplex Map + \image html "map.png" "Toplex map representation" +<table border="0"> + <tr> + <td width="25%"> + <b>Author:</b> François Godi<br> + <b>Introduced in:</b> GUDHI 2.1.0<br> + <b>Copyright:</b> GPL v3<br> + </td> + <td width="75%"> + The Toplex map data structure is composed firstly of a raw storage of toplices (the maximal simplices) + and secondly of a map which associate any vertex to a set of pointers toward all toplices + containing this vertex. + <b>User manual:</b> \ref toplex_map - <b>Reference manual:</b> Gudhi::Toplex_map + </td> + </tr> + \subsection WitnessComplexDataStructure Witness complex \image html "Witness_complex_representation.png" "Witness complex representation" <table border="0"> diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt index 97859c98..480332d7 100644 --- a/src/cython/CMakeLists.txt +++ b/src/cython/CMakeLists.txt @@ -215,12 +215,7 @@ if(PYTHONINTERP_FOUND) add_custom_target(cython ALL DEPENDS gudhi.so COMMENT "Do not forget to add ${CMAKE_CURRENT_BINARY_DIR}/ to your PYTHONPATH before using examples or tests") - # For installation purpose - # TODO(VR) : files matching pattern mechanism is copying all cython directory - install(DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" DESTINATION "${PYTHON_SITE_PACKAGES}/" FILES_MATCHING - PATTERN "*.so" - PATTERN "*.dylib" - PATTERN "*.pyd") + install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/setup.py install)") # Test examples if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/src/cython/cython/tangential_complex.pyx b/src/cython/cython/tangential_complex.pyx index 4bb07076..293ef8cb 100644 --- a/src/cython/cython/tangential_complex.pyx +++ b/src/cython/cython/tangential_complex.pyx @@ -36,6 +36,7 @@ cdef extern from "Tangential_complex_interface.h" namespace "Gudhi": Tangential_complex_interface(int intrisic_dim, vector[vector[double]] points) # bool from_file is a workaround for cython to find the correct signature Tangential_complex_interface(int intrisic_dim, string off_file, bool from_file) + void compute_tangential_complex() except + vector[double] get_point(unsigned vertex) unsigned number_of_vertices() unsigned number_of_simplices() @@ -43,6 +44,7 @@ cdef extern from "Tangential_complex_interface.h" namespace "Gudhi": unsigned number_of_inconsistent_stars() void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree) void fix_inconsistencies_using_perturbation(double max_perturb, double time_limit) + void set_max_squared_edge_length(double max_squared_edge_length) # TangentialComplex python interface cdef class TangentialComplex: @@ -92,6 +94,17 @@ cdef class TangentialComplex: """ return self.thisptr != NULL + def compute_tangential_complex(self): + """This function computes the tangential complex. + + Raises: + ValueError: In debug mode, if the computed star dimension is too + low. Try to set a bigger maximal edge length value with + :func:`~gudhi.Tangential_complex.set_max_squared_edge_length` + if this happens. + """ + self.thisptr.compute_tangential_complex() + def get_point(self, vertex): """This function returns the point corresponding to a given vertex. @@ -152,3 +165,16 @@ cdef class TangentialComplex: """ self.thisptr.fix_inconsistencies_using_perturbation(max_perturb, time_limit) + + def set_max_squared_edge_length(self, max_squared_edge_length): + """Sets the maximal possible squared edge length for the edges in the + triangulations. + + :param max_squared_edge_length: Maximal possible squared edge length. + :type max_squared_edge_length: double + + If the maximal edge length value is too low + :func:`~gudhi.Tangential_complex.compute_tangential_complex` + will throw an exception in debug mode. + """ + self.thisptr.set_max_squared_edge_length(max_squared_edge_length) diff --git a/src/cython/doc/installation.rst b/src/cython/doc/installation.rst index 040f6b4a..855dea44 100644 --- a/src/cython/doc/installation.rst +++ b/src/cython/doc/installation.rst @@ -7,24 +7,23 @@ Installation Compiling ********* -The library uses c++11 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0 -and `CMake <https://www.cmake.org/>`_ ≥ 3.1. +The library uses c++11 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0, +`CMake <https://www.cmake.org/>`_ ≥ 3.1 to generate makefiles, and +`Cython <https://www.cython.org/>`_ to compile the GUDHI Python module. It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2015. -It also requires cmake to generate makefiles, and cython to compile the -library. On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python 3.5 and 3.6 are available because of the required Visual Studio version. -On other systems, if you have several Python/cython installed, the version 2.X +On other systems, if you have several Python/Cython installed, the version 2.X will be used by default, but you can force it by adding :code:`-DPython_ADDITIONAL_VERSIONS=3` to the cmake command. -GUDHI Cythonization -=================== +GUDHI Python module compilation +=============================== -To build the GUDHI cython module, run the following commands in a terminal: +To build the GUDHI Python module, run the following commands in a terminal: .. code-block:: bash @@ -32,7 +31,28 @@ To build the GUDHI cython module, run the following commands in a terminal: mkdir build cd build/ cmake .. - make cython + cd cython + make + +GUDHI Python module installation +================================ + +Once the compilation succeeds, one can add the GUDHI Python module path to the +PYTHONPATH: + +.. code-block:: bash + + # For windows, you have to set PYTHONPATH environment variable + export PYTHONPATH='$PYTHONPATH:/path-to-gudhi/build/cython' + +Or install it definitely in your Python packages folder: + +.. code-block:: bash + + cd /path-to-gudhi/build/cython + # May require sudo or administrator privileges + make install + Test suites =========== @@ -45,7 +65,7 @@ following command in a terminal: cd /path-to-gudhi/build/cython # For windows, you have to set PYTHONPATH environment variable export PYTHONPATH='$PYTHONPATH:/path-to-gudhi/build/cython' - ctest -R py_test + make test Debugging issues ================ @@ -54,7 +74,7 @@ If tests fail, please check your PYTHONPATH and try to :code:`import gudhi` and check the errors. The problem can come from a third-party library bad link or installation. -If :code:`import gudhi` succeeds, please have a look to debug informations: +If :code:`import gudhi` succeeds, please have a look to debug information: .. code-block:: python @@ -105,13 +125,17 @@ A complete configuration would be : Documentation ============= -To build the documentation, `sphinx-doc <http://http://www.sphinx-doc.org>`_ is -required. Please refer to *conf.py* file to see which -`sphinx-doc <http://http://www.sphinx-doc.org>`_ modules are required to -generate the documentation. Run the following commands in a terminal: +To build the documentation, `sphinx-doc <http://www.sphinx-doc.org>`_ and +`sphinxcontrib-bibtex <https://sphinxcontrib-bibtex.readthedocs.io>`_ are +required. As the documentation is auto-tested, `CGAL`_, `Eigen3`_, +`Matplotlib`_, `NumPy`_ and `SciPy`_ are also mandatory to build the +documentation. + +Run the following commands in a terminal: .. code-block:: bash + cd /path-to-gudhi/build/cython make sphinx Optional third-party library @@ -127,7 +151,7 @@ The :doc:`Alpha complex </alpha_complex_user>`, C++ library which provides easy access to efficient and reliable geometric algorithms. -Having CGAL, the Computational Geometry Algorithms Library, version 4.7.0 or +Having CGAL, the Computational Geometry Algorithms Library, version 4.7.0 or higher installed is recommended. The procedure to install this library according to your operating system is detailed `here <http://doc.cgal.org/latest/Manual/installation.html>`_. diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index e814b4c3..1d340dbe 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -50,8 +50,8 @@ by more than the length used to define "too close". A more general technique is to use a sparse approximation of the Rips introduced by Don Sheehy :cite:`sheehy13linear`. We are using the version described in :cite:`buchet16efficient` (except that we multiply all filtration -values by 2, to match the usual Rips complex), which proves a -:math:`\frac{1+\varepsilon}{1-\varepsilon}`-interleaving, although in practice the +values by 2, to match the usual Rips complex). :cite:`cavanna15geometric` proves +a :math:`\frac{1}{1-\varepsilon}`-interleaving, although in practice the error is usually smaller. A more intuitive presentation of the idea is available in :cite:`cavanna15geometric`, and in a video :cite:`cavanna15visualizing`. Passing an extra argument `sparse=0.3` at the diff --git a/src/cython/doc/tangential_complex_user.rst b/src/cython/doc/tangential_complex_user.rst index 5ce69e86..ebfe1e29 100644 --- a/src/cython/doc/tangential_complex_user.rst +++ b/src/cython/doc/tangential_complex_user.rst @@ -23,8 +23,10 @@ What is a Tangential Complex? ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Let us start with the description of the Tangential complex of a simple -example, with :math:`k = 1` and :math:`d = 2`. The input data is 4 points -:math:`P` located on a curve embedded in 2D. +example, with :math:`k = 1` and :math:`d = 2`. The point set +:math:`\mathscr P` is located on a closed curve embedded in 2D. +Only 4 points will be displayed (more are required for PCA) to simplify the +figures. .. figure:: ../../doc/Tangential_complex/tc_example_01.png :alt: The input @@ -32,8 +34,7 @@ example, with :math:`k = 1` and :math:`d = 2`. The input data is 4 points The input -For each point :math:`p`, estimate its tangent subspace :math:`T_p` (e.g. -using PCA). +For each point :math:`P`, estimate its tangent subspace :math:`T_P` using PCA. .. figure:: ../../doc/Tangential_complex/tc_example_02.png :alt: The estimated normals @@ -43,8 +44,8 @@ using PCA). Let us add the Voronoi diagram of the points in orange. For each point -:math:`p`, construct its star in the Delaunay triangulation of :math:`P` -restricted to :math:`T_p`. +:math:`P`, construct its star in the Delaunay triangulation of +:math:`\mathscr P` restricted to :math:`T_P`. .. figure:: ../../doc/Tangential_complex/tc_example_03.png :alt: The Voronoi diagram @@ -72,7 +73,7 @@ Let us take the same example. Before -Let us slightly move the tangent subspace :math:`T_q` +Let us slightly move the tangent subspace :math:`T_Q` .. figure:: ../../doc/Tangential_complex/tc_example_07_after.png :alt: After @@ -128,6 +129,7 @@ This example builds the Tangential complex of point set read in an OFF file. import gudhi tc = gudhi.TangentialComplex(intrisic_dim = 1, off_file=gudhi.__root_source_dir__ + '/data/points/alphacomplexdoc.off') + tc.compute_tangential_complex() result_str = 'Tangential contains ' + repr(tc.num_simplices()) + \ ' simplices - ' + repr(tc.num_vertices()) + ' vertices.' print(result_str) @@ -175,6 +177,7 @@ simplices. import gudhi tc = gudhi.TangentialComplex(intrisic_dim = 1, points=[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + tc.compute_tangential_complex() result_str = 'Tangential contains ' + repr(tc.num_vertices()) + ' vertices.' print(result_str) diff --git a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py index 0f8f5e80..536517d1 100755 --- a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py +++ b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py @@ -50,6 +50,7 @@ with open(args.file, 'r') as f: print("TangentialComplex creation from points read in a OFF file") tc = gudhi.TangentialComplex(intrisic_dim = args.intrisic_dim, off_file=args.file) + tc.compute_tangential_complex() st = tc.create_simplex_tree() message = "Number of simplices=" + repr(st.num_simplices()) diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 1a6e2477..40aff299 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -54,23 +54,17 @@ class Rips_complex_interface { } void init_points_sparse(const std::vector<std::vector<double>>& points, double threshold, double epsilon) { - sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon); - threshold_ = threshold; + sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon, -std::numeric_limits<double>::infinity(), threshold); } void init_matrix_sparse(const std::vector<std::vector<double>>& matrix, double threshold, double epsilon) { - sparse_rips_complex_.emplace(matrix, epsilon); - threshold_ = threshold; + sparse_rips_complex_.emplace(matrix, epsilon, -std::numeric_limits<double>::infinity(), threshold); } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) { if (rips_complex_) rips_complex_->create_complex(*simplex_tree, dim_max); - else { + else sparse_rips_complex_->create_complex(*simplex_tree, dim_max); - // This pruning should be done much earlier! It isn't that useful for sparse Rips, - // but it would be inconsistent not to do it. - simplex_tree->prune_above_filtration(threshold_); - } simplex_tree->initialize_filtration(); } @@ -79,7 +73,6 @@ class Rips_complex_interface { // Anyway, storing a graph would make more sense. Or changing the interface completely so there is no such storage. boost::optional<Rips_complex<Simplex_tree_interface<>::Filtration_value>> rips_complex_; boost::optional<Sparse_rips_complex<Simplex_tree_interface<>::Filtration_value>> sparse_rips_complex_; - double threshold_; }; } // namespace rips_complex diff --git a/src/cython/include/Tangential_complex_interface.h b/src/cython/include/Tangential_complex_interface.h index 71418886..c4ddbdbe 100644 --- a/src/cython/include/Tangential_complex_interface.h +++ b/src/cython/include/Tangential_complex_interface.h @@ -49,8 +49,6 @@ class Tangential_complex_interface { Dynamic_kernel k; tangential_complex_ = new TC(points, intrisic_dim, k); - tangential_complex_->compute_tangential_complex(); - num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); } Tangential_complex_interface(int intrisic_dim, const std::string& off_file_name, bool from_file = true) { @@ -60,14 +58,17 @@ class Tangential_complex_interface { std::vector<Point_d> points = off_reader.get_point_cloud(); tangential_complex_ = new TC(points, intrisic_dim, k); - tangential_complex_->compute_tangential_complex(); - num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); } ~Tangential_complex_interface() { delete tangential_complex_; } + void compute_tangential_complex() { + tangential_complex_->compute_tangential_complex(); + num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); + } + std::vector<double> get_point(unsigned vh) { std::vector<double> vd; if (vh < tangential_complex_->number_of_vertices()) { @@ -104,7 +105,11 @@ class Tangential_complex_interface { simplex_tree->initialize_filtration(); } - private: + void set_max_squared_edge_length(double max_squared_edge_length) { + tangential_complex_->set_max_squared_edge_length(max_squared_edge_length); + } + +private: TC* tangential_complex_; TC::Num_inconsistencies num_inconsistencies_; }; diff --git a/src/cython/test/test_tangential_complex.py b/src/cython/test/test_tangential_complex.py index 5385a0d3..5c62f278 100755 --- a/src/cython/test/test_tangential_complex.py +++ b/src/cython/test/test_tangential_complex.py @@ -32,6 +32,15 @@ def test_tangential(): tc = TangentialComplex(intrisic_dim = 1, points=point_list) assert tc.__is_defined() == True assert tc.num_vertices() == 4 + assert tc.num_simplices() == 0 + assert tc.num_inconsistent_simplices() == 0 + assert tc.num_inconsistent_stars() == 0 + + tc.compute_tangential_complex() + assert tc.num_vertices() == 4 + assert tc.num_simplices() == 4 + assert tc.num_inconsistent_simplices() == 0 + assert tc.num_inconsistent_stars() == 0 st = tc.create_simplex_tree() assert st.__is_defined() == True |