diff options
37 files changed, 335 insertions, 234 deletions
diff --git a/.github/workflows/pip-build-linux.yml b/.github/workflows/pip-build-linux.yml index a2b4f085..11b6271d 100644 --- a/.github/workflows/pip-build-linux.yml +++ b/.github/workflows/pip-build-linux.yml @@ -9,7 +9,7 @@ jobs: # cf. https://github.com/GUDHI/gudhi-deploy/blob/main/Dockerfile_for_pip container: gudhi/pip_for_gudhi steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - name: Build wheel for Python 3.10 diff --git a/.github/workflows/pip-build-osx.yml b/.github/workflows/pip-build-osx.yml index 99d515ff..59e94ca5 100644 --- a/.github/workflows/pip-build-osx.yml +++ b/.github/workflows/pip-build-osx.yml @@ -11,10 +11,10 @@ jobs: python-version: ['3.10'] name: Build wheels for Python ${{ matrix.python-version }} steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - - uses: actions/setup-python@v1 + - uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} architecture: x64 diff --git a/.github/workflows/pip-build-windows.yml b/.github/workflows/pip-build-windows.yml index 30b0bd94..b3d75706 100644 --- a/.github/workflows/pip-build-windows.yml +++ b/.github/workflows/pip-build-windows.yml @@ -11,10 +11,10 @@ jobs: python-version: ['3.10'] name: Build wheels for Python ${{ matrix.python-version }} steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - - uses: actions/setup-python@v1 + - uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} architecture: x64 diff --git a/.github/workflows/pip-packaging-linux.yml b/.github/workflows/pip-packaging-linux.yml index 98173ed3..285cfa00 100644 --- a/.github/workflows/pip-packaging-linux.yml +++ b/.github/workflows/pip-packaging-linux.yml @@ -11,7 +11,7 @@ jobs: # cf. https://github.com/GUDHI/gudhi-deploy/blob/main/Dockerfile_for_pip container: gudhi/pip_for_gudhi steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - name: Build wheel for Python 3.6 diff --git a/.github/workflows/pip-packaging-osx.yml b/.github/workflows/pip-packaging-osx.yml index 7417300a..3ae840c6 100644 --- a/.github/workflows/pip-packaging-osx.yml +++ b/.github/workflows/pip-packaging-osx.yml @@ -13,10 +13,10 @@ jobs: python-version: ['3.7', '3.8', '3.9', '3.10'] name: Build wheels for Python ${{ matrix.python-version }} steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - - uses: actions/setup-python@v1 + - uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} architecture: x64 diff --git a/.github/workflows/pip-packaging-windows.yml b/.github/workflows/pip-packaging-windows.yml index 48a98036..6f544499 100644 --- a/.github/workflows/pip-packaging-windows.yml +++ b/.github/workflows/pip-packaging-windows.yml @@ -13,10 +13,10 @@ jobs: python-version: ['3.7', '3.8', '3.9', '3.10'] name: Build wheels for Python ${{ matrix.python-version }} steps: - - uses: actions/checkout@v1 + - uses: actions/checkout@v3 with: submodules: true - - uses: actions/setup-python@v1 + - uses: actions/setup-python@v4 with: python-version: ${{ matrix.python-version }} architecture: x64 diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index aec8c1b1..a7372f19 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -17,8 +17,7 @@ // to construct Alpha_complex from a OFF file of points #include <gudhi/Points_off_io.h> -#include <stdlib.h> -#include <math.h> // isnan, fmax +#include <cmath> // isnan, fmax #include <memory> // for std::unique_ptr #include <cstddef> // for std::size_t @@ -45,6 +44,7 @@ #include <utility> // std::pair #include <stdexcept> #include <numeric> // for std::iota +#include <algorithm> // for std::sort // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 #if CGAL_VERSION_NR < 1041101000 @@ -101,13 +101,17 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val */ template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false> class Alpha_complex { + private: + // Vertex_handle internal type (required by triangulation_ and vertices_). + using Internal_vertex_handle = std::ptrdiff_t; + public: /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/ using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>; // Add an int in TDS to save point index in the structure using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension, - CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>, + CGAL::Triangulation_vertex<Geom_traits, Internal_vertex_handle>, CGAL::Triangulation_full_cell<Geom_traits> >; /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ @@ -132,9 +136,6 @@ class Alpha_complex { // Vertex_iterator type from CGAL. using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; - // size_type type from CGAL. - using size_type = typename Triangulation::size_type; - // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >; @@ -146,6 +147,10 @@ class Alpha_complex { std::unique_ptr<Triangulation> triangulation_; /** \brief Kernel for triangulation_ functions access.*/ A_kernel_d kernel_; + /** \brief Vertices to be inserted first by the create_complex method to avoid quadratic complexity. + * It isn't just [0, n) if some points have multiplicity (only one copy appears in the complex). + */ + std::vector<Internal_vertex_handle> vertices_; /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ std::vector<Sphere> cache_, old_cache_; @@ -257,11 +262,11 @@ class Alpha_complex { std::vector<Point_d> point_cloud(first, last); // Creates a vector {0, 1, ..., N-1} - std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0), - boost::counting_iterator<std::ptrdiff_t>(point_cloud.size())); + std::vector<Internal_vertex_handle> indices(boost::counting_iterator<Internal_vertex_handle>(0), + boost::counting_iterator<Internal_vertex_handle>(point_cloud.size())); using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator, - CGAL::Identity_property_map<std::ptrdiff_t>>; + CGAL::Identity_property_map<Internal_vertex_handle>>; using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>; CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); @@ -279,6 +284,9 @@ class Alpha_complex { // structure to retrieve CGAL points from vertex handle - one vertex handle per point. // Needs to be constructed before as vertex handles arrives in no particular order. vertex_handle_to_iterator_.resize(point_cloud.size()); + // List of sorted unique vertices in the triangulation. We take advantage of the existing loop to construct it + // Vertices list avoids quadratic complexity with the Simplex_tree. We should not fill it up with Toplex_map e.g. + vertices_.reserve(triangulation_->number_of_vertices()); // Loop on triangulation vertices list for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { if (!triangulation_->is_infinite(*vit)) { @@ -286,8 +294,10 @@ class Alpha_complex { std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; #endif // DEBUG_TRACES vertex_handle_to_iterator_[vit->data()] = vit; + vertices_.push_back(vit->data()); } } + std::sort(vertices_.begin(), vertices_.end()); // -------------------------------------------------------------------------------------------- } } @@ -384,12 +394,21 @@ class Alpha_complex { // -------------------------------------------------------------------------------------------- // Simplex_tree construction from loop on triangulation finite full cells list if (num_vertices() > 0) { + std::vector<Vertex_handle> one_vertex(1); + for (auto vertex : vertices_) { +#ifdef DEBUG_TRACES + std::clog << "SimplicialComplex insertion " << vertex << std::endl; +#endif // DEBUG_TRACES + one_vertex[0] = vertex; + complex.insert_simplex_and_subfaces(one_vertex, std::numeric_limits<double>::quiet_NaN()); + } + for (auto cit = triangulation_->finite_full_cells_begin(); cit != triangulation_->finite_full_cells_end(); ++cit) { Vector_vertex vertexVector; #ifdef DEBUG_TRACES - std::clog << "Simplex_tree insertion "; + std::clog << "SimplicialComplex insertion "; #endif // DEBUG_TRACES for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { if (*vit != nullptr) { diff --git a/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp new file mode 100644 index 00000000..e7c261f1 --- /dev/null +++ b/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp @@ -0,0 +1,117 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_dim3" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> + +#include <stdexcept> // std::out_of_range +#include <string> +#include <vector> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Simplex_tree.h> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > Exact_kernel_s; +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Inexact_kernel_s; +// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter + +typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double max_alpha_square_value = 60.0; + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); + + Gudhi::Simplex_tree<> simplex_tree_60; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); + + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + + std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; + BOOST_CHECK(simplex_tree_60.dimension() == 2); + + std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_vertices() == 7); + + std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_simplices() == 25); + + max_alpha_square_value = 59.0; + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::Simplex_tree<> simplex_tree_59; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); + + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + + std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; + BOOST_CHECK(simplex_tree_59.dimension() == 2); + + std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_vertices() == 7); + + std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_simplices() == 23); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { + std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl; + + // ---------------------------------------------------------------------------- + // Init of an empty list of points + // ---------------------------------------------------------------------------- + std::vector<typename TestedKernel::Point_d> points; + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); + + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); + + // Test to the limit + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); + + Gudhi::Simplex_tree<> simplex_tree; + BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); + + std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree.num_simplices() == 0); + + std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + BOOST_CHECK(simplex_tree.dimension() == -1); + + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree.num_vertices() == points.size()); +} diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index f74ad217..b474917f 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -13,81 +13,17 @@ #include <boost/test/unit_test.hpp> #include <boost/mpl/list.hpp> -#include <CGAL/Delaunay_triangulation.h> #include <CGAL/Epick_d.h> #include <CGAL/Epeck_d.h> -#include <cmath> // float comparison -#include <limits> +#include <stdexcept> // std::out_of_range #include <string> #include <vector> #include <gudhi/Alpha_complex.h> -// to construct a simplex_tree from Delaunay_triangulation -#include <gudhi/graph_simplicial_complex.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Unitary_tests_utils.h> -// Use dynamic_dimension_tag for the user to be able to set dimension -typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; -// Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > Exact_kernel_s; -// Use dynamic_dimension_tag for the user to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; -// Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Inexact_kernel_s; -// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter - -typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; - -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { - // ---------------------------------------------------------------------------- - // - // Init of an alpha-complex from a OFF file - // - // ---------------------------------------------------------------------------- - std::string off_file_name("alphacomplexdoc.off"); - double max_alpha_square_value = 60.0; - std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << - max_alpha_square_value << "==========" << std::endl; - - Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); - - Gudhi::Simplex_tree<> simplex_tree_60; - BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); - - std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; - BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); - - std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; - BOOST_CHECK(simplex_tree_60.dimension() == 2); - - std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree_60.num_vertices() == 7); - - std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree_60.num_simplices() == 25); - - max_alpha_square_value = 59.0; - std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << - max_alpha_square_value << "==========" << std::endl; - - Gudhi::Simplex_tree<> simplex_tree_59; - BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); - - std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; - BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); - - std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; - BOOST_CHECK(simplex_tree_59.dimension() == 2); - - std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree_59.num_vertices() == 7); - - std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree_59.num_simplices() == 23); -} - // Use static dimension_tag for the user not to be able to set dimension typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Kernel_4; typedef Kernel_4::Point_d Point_4; @@ -236,37 +172,6 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { } -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { - std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl; - - // ---------------------------------------------------------------------------- - // Init of an empty list of points - // ---------------------------------------------------------------------------- - std::vector<typename TestedKernel::Point_d> points; - - // ---------------------------------------------------------------------------- - // Init of an alpha complex from the list of points - // ---------------------------------------------------------------------------- - Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); - - std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; - BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); - - // Test to the limit - BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); - - Gudhi::Simplex_tree<> simplex_tree; - BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); - - std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree.num_simplices() == 0); - - std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; - BOOST_CHECK(simplex_tree.dimension() == -1); - - std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree.num_vertices() == points.size()); -} using Inexact_kernel_2 = CGAL::Epick_d< CGAL::Dimension_tag<2> >; using Exact_kernel_2 = CGAL::Epeck_d< CGAL::Dimension_tag<2> >; diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index 0595ca92..dd2c235f 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -8,14 +8,18 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp ) target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY}) + add_executable ( Alpha_complex_dim3_test_unit Alpha_complex_dim3_unit_test.cpp ) + target_link_libraries(Alpha_complex_dim3_test_unit ${CGAL_LIBRARY}) add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp ) target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Alpha_complex_dim3_test_unit ${TBB_LIBRARIES}) target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) + gudhi_add_boost_test(Alpha_complex_dim3_test_unit) gudhi_add_boost_test(Delaunay_complex_test_unit) endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) @@ -73,4 +77,4 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) endif() gudhi_add_boost_test(Zero_weighted_alpha_complex_test_unit) -endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
\ No newline at end of file +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index 91899040..e65d8c6f 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -263,7 +263,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "cuboid-file,c", po::value<std::string>(&cuboid_file), "Name of file describing the periodic domain. Format is:\n min_hx min_hy min_hz\n max_hx max_hy max_hz")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index e86b34e2..29edbd8e 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -163,7 +163,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "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 ")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal alpha square value for the Alpha complex construction.")( diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index 4a6af3a4..29fabc6c 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -241,10 +241,16 @@ class Bitmap_cubical_complex : public T { **/ class Filtration_simplex_range; - class Filtration_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> { + class Filtration_simplex_iterator { // Iterator over all simplices of the complex in the order of the indexing scheme. // 'value_type' must be 'Simplex_handle'. public: + typedef std::input_iterator_tag iterator_category; + typedef Simplex_handle value_type; + typedef std::ptrdiff_t difference_type; + typedef value_type* pointer; + typedef value_type reference; + Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) {} Filtration_simplex_iterator() : b(NULL), position(0) {} @@ -386,10 +392,16 @@ class Bitmap_cubical_complex : public T { **/ class Skeleton_simplex_range; - class Skeleton_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> { + class Skeleton_simplex_iterator { // Iterator over all simplices of the complex in the order of the indexing scheme. // 'value_type' must be 'Simplex_handle'. public: + typedef std::input_iterator_tag iterator_category; + typedef Simplex_handle value_type; + typedef std::ptrdiff_t difference_type; + typedef value_type* pointer; + typedef value_type reference; + Skeleton_simplex_iterator(Bitmap_cubical_complex* b, std::size_t d) : b(b), dimension(d) { if (globalDbg) { std::clog << "Skeleton_simplex_iterator ( Bitmap_cubical_complex* b , std::size_t d )\n"; diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index bafe7981..2bf62f9b 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -251,8 +251,14 @@ class Bitmap_cubical_complex_base { * @brief Iterator through all cells in the complex (in order they appear in the structure -- i.e. * in lexicographical order). **/ - class All_cells_iterator : std::iterator<std::input_iterator_tag, T> { + class All_cells_iterator { public: + typedef std::input_iterator_tag iterator_category; + typedef std::size_t value_type; + typedef std::ptrdiff_t difference_type; + typedef value_type* pointer; + typedef value_type reference; + All_cells_iterator() { this->counter = 0; } All_cells_iterator operator++() { @@ -355,8 +361,14 @@ class Bitmap_cubical_complex_base { * @brief Iterator through top dimensional cells of the complex. The cells appear in order they are stored * in the structure (i.e. in lexicographical order) **/ - class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> { + class Top_dimensional_cells_iterator { public: + typedef std::input_iterator_tag iterator_category; + typedef std::size_t value_type; + typedef std::ptrdiff_t difference_type; + typedef value_type* pointer; + typedef value_type reference; + Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) { this->counter = std::vector<std::size_t>(b.dimension()); // std::fill( this->counter.begin() , this->counter.end() , 0 ); diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index a07ba212..e6419f3d 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -118,7 +118,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, bool& "fast,f", po::bool_switch(&fast), "To activate fast version of Cech complex (default is false, not available if exact is set)")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-radius,r", po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Cech complex construction.")( diff --git a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp index 38efb9e6..70b489b5 100644 --- a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp +++ b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp @@ -111,7 +111,7 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.")( diff --git a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp index d8f42ab6..a8fd6f14 100644 --- a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp +++ b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp @@ -140,7 +140,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.")( diff --git a/src/Persistent_cohomology/example/persistence_from_file.cpp b/src/Persistent_cohomology/example/persistence_from_file.cpp index 38c44514..7f89c001 100644 --- a/src/Persistent_cohomology/example/persistence_from_file.cpp +++ b/src/Persistent_cohomology/example/persistence_from_file.cpp @@ -93,7 +93,7 @@ void program_options(int argc, char * argv[] visible.add_options() ("help,h", "produce help message") ("output-file,o", po::value<std::string>(&output_file)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog") + "Name of file in which the persistence diagram is written. Default print in standard output") ("field-charac,p", po::value<int>(&p)->default_value(11), "Characteristic p of the coefficient field Z/pZ for computing homology.") ("min-persistence,m", po::value<Filtration_value>(&min_persistence), diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index ca26a5b9..84453898 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -96,7 +96,7 @@ void program_options(int argc, char * argv[] visible.add_options() ("help,h", "produce help message") ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog") + "Name of file in which the persistence diagram is written. Default print in standard output") ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), "Maximal length of an edge for the Rips complex construction.") ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp index a503d983..6f37cf5c 100644 --- a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -112,7 +112,7 @@ void program_options(int argc, char * argv[] visible.add_options() ("help,h", "produce help message") ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog") + "Name of file in which the persistence diagram is written. Default print in standard output") ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.") diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp index 8c5742aa..6b60f603 100644 --- a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -109,7 +109,7 @@ void program_options(int argc, char * argv[] visible.add_options() ("help,h", "produce help message") ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog") + "Name of file in which the persistence diagram is written. Default print in standard output") ("max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(0), "Maximal length of an edge for the Rips complex construction.") ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), diff --git a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp index b473738e..72ddc797 100644 --- a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -118,7 +118,7 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "min-edge-corelation,c", po::value<Filtration_value>(&correlation_min)->default_value(0), "Minimal corelation of an edge for the Rips complex construction.")( "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), diff --git a/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp index 6306755d..77ad841a 100644 --- a/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp +++ b/src/Rips_complex/utilities/rips_distance_matrix_persistence.cpp @@ -79,7 +79,7 @@ void program_options(int argc, char* argv[], std::string& csv_matrix_file, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.")( diff --git a/src/Rips_complex/utilities/rips_persistence.cpp b/src/Rips_complex/utilities/rips_persistence.cpp index 9d7490b3..43194821 100644 --- a/src/Rips_complex/utilities/rips_persistence.cpp +++ b/src/Rips_complex/utilities/rips_persistence.cpp @@ -81,7 +81,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.")( diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp index ac935b41..829c85e6 100644 --- a/src/Rips_complex/utilities/sparse_rips_persistence.cpp +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -84,7 +84,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-edge-length,r", po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal length of an edge for the Rips complex construction.")( diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index cc424810..56a24af0 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -36,7 +36,6 @@ #include <Eigen/Eigen> #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST -#include <boost/optional.hpp> #include <boost/iterator/transform_iterator.hpp> #include <boost/range/adaptor/transformed.hpp> #include <boost/range/counting_range.hpp> @@ -56,6 +55,7 @@ #include <cmath> // for std::sqrt #include <string> #include <cstddef> // for std::size_t +#include <optional> #ifdef GUDHI_USE_TBB #include <tbb/parallel_for.h> @@ -994,7 +994,7 @@ class Tangential_complex { // circumspheres of the star of "center_vertex" // If th the m_max_squared_edge_length is set the maximal radius of the "star sphere" // is at most square root of m_max_squared_edge_length - boost::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length; + std::optional<FT> squared_star_sphere_radius_plus_margin = m_max_squared_edge_length; // Insert points until we find a point which is outside "star sphere" for (auto nn_it = ins_range.begin(); nn_it != ins_range.end(); ++nn_it) { @@ -1036,7 +1036,7 @@ class Tangential_complex { // Let's recompute squared_star_sphere_radius_plus_margin if (triangulation.current_dimension() >= tangent_space_dim) { - squared_star_sphere_radius_plus_margin = boost::none; + squared_star_sphere_radius_plus_margin = std::nullopt; // Get the incident cells and look for the biggest circumsphere std::vector<Tr_full_cell_handle> incident_cells; triangulation.incident_full_cells(center_vertex, std::back_inserter(incident_cells)); @@ -1044,7 +1044,7 @@ class Tangential_complex { cit != incident_cells.end(); ++cit) { Tr_full_cell_handle cell = *cit; if (triangulation.is_infinite(cell)) { - squared_star_sphere_radius_plus_margin = boost::none; + squared_star_sphere_radius_plus_margin = std::nullopt; break; } else { // Note that this uses the perturbed point since it uses @@ -2030,7 +2030,7 @@ class Tangential_complex { // and their center vertex Stars_container m_stars; std::vector<FT> m_squared_star_spheres_radii_incl_margin; - boost::optional<FT> m_max_squared_edge_length; + std::optional<FT> m_max_squared_edge_length; #ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM Points m_points_for_tse; diff --git a/src/Witness_complex/utilities/strong_witness_persistence.cpp b/src/Witness_complex/utilities/strong_witness_persistence.cpp index 614de0d4..b2ecad82 100644 --- a/src/Witness_complex/utilities/strong_witness_persistence.cpp +++ b/src/Witness_complex/utilities/strong_witness_persistence.cpp @@ -108,7 +108,7 @@ void program_options(int argc, char* argv[], int& nbL, std::string& file_name, s visible.add_options()("help,h", "produce help message")("landmarks,l", po::value<int>(&nbL), "Number of landmarks to choose from the point cloud.")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha), "Maximal squared relaxation parameter.")( "field-charac,p", po::value<int>(&p)->default_value(11), diff --git a/src/Witness_complex/utilities/weak_witness_persistence.cpp b/src/Witness_complex/utilities/weak_witness_persistence.cpp index 5ea31d6b..c7ead7de 100644 --- a/src/Witness_complex/utilities/weak_witness_persistence.cpp +++ b/src/Witness_complex/utilities/weak_witness_persistence.cpp @@ -108,7 +108,7 @@ void program_options(int argc, char* argv[], int& nbL, std::string& file_name, s visible.add_options()("help,h", "produce help message")("landmarks,l", po::value<int>(&nbL), "Number of landmarks to choose from the point cloud.")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha), "Maximal squared relaxation parameter.")( "field-charac,p", po::value<int>(&p)->default_value(11), diff --git a/src/Witness_complex/utilities/witnesscomplex.md b/src/Witness_complex/utilities/witnesscomplex.md index 3a3a7d83..e994e0b8 100644 --- a/src/Witness_complex/utilities/witnesscomplex.md +++ b/src/Witness_complex/utilities/witnesscomplex.md @@ -29,7 +29,7 @@ and `p` is the characteristic of the field *Z/pZ* used for homology coefficients * `-h [ --help ]` Produce help message
* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud.
-* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::clog.
+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in standard output.
* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter.
* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology.
* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
@@ -60,7 +60,7 @@ and `p` is the characteristic of the field *Z/pZ* used for homology coefficients * `-h [ --help ]` Produce help message
* `-l [ --landmarks ]` Number of landmarks to choose from the point cloud.
-* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::clog.
+* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in standard output.
* `-a [ --max-sq-alpha ]` (default = inf) Maximal squared relaxation parameter.
* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology.
* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
diff --git a/src/cmake/modules/GUDHI_compilation_flags.cmake b/src/cmake/modules/GUDHI_compilation_flags.cmake index e2d3d872..b43ccf73 100644 --- a/src/cmake/modules/GUDHI_compilation_flags.cmake +++ b/src/cmake/modules/GUDHI_compilation_flags.cmake @@ -12,6 +12,7 @@ macro(add_cxx_compiler_flag _flag) endmacro() set (CMAKE_CXX_STANDARD 17) +# This number needs to be changed in python/CMakeLists.txt at the same time enable_testing() diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 5f323935..8f8df138 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -129,9 +129,10 @@ if(PYTHONINTERP_FOUND) # Gudhi and CGAL compilation option if(MSVC) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/std:c++17', ") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/fp:strict', ") else(MSVC) - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++14', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++17', ") endif(MSVC) if(CMAKE_COMPILER_IS_GNUCXX) set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-frounding-math', ") @@ -246,8 +247,8 @@ if(PYTHONINTERP_FOUND) # Specific for Mac if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.12', ") - set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.12', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.14', ") + set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.14', ") endif(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") # Loop on INCLUDE_DIRECTORIES PROPERTY diff --git a/src/python/doc/clustering.rst b/src/python/doc/clustering.rst index c5a57d3c..62422682 100644 --- a/src/python/doc/clustering.rst +++ b/src/python/doc/clustering.rst @@ -17,9 +17,8 @@ As a by-product, we produce the persistence diagram of the merge tree of the ini :include-source: import gudhi - f = open(gudhi.__root_source_dir__ + '/data/points/spiral_2d.csv', 'r') - import numpy as np - data = np.loadtxt(f) + from gudhi.datasets.remote import fetch_spiral_2d + data = fetch_spiral_2d() import matplotlib.pyplot as plt plt.scatter(data[:,0],data[:,1],marker='.',s=1) plt.show() diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index b704f778..276ac4e2 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -391,7 +391,7 @@ The :doc:`persistence graphical tools </persistence_graphical_tools_user>` and mathematics, science, and engineering. :class:`~gudhi.point_cloud.knn.KNearestNeighbors` can use the Python package -`SciPy <http://scipy.org>`_ as a backend if explicitly requested. +`SciPy <http://scipy.org>`_ :math:`\geq` 1.6.0 as a backend if explicitly requested. TensorFlow ---------- diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py index 21275cdd..e438aa66 100644 --- a/src/python/gudhi/persistence_graphical_tools.py +++ b/src/python/gudhi/persistence_graphical_tools.py @@ -194,19 +194,21 @@ def plot_persistence_barcode( y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence] c=[colormap[dim] for (dim,(birth,death)) in persistence] - axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0) + axes.barh(range(len(x)), y, left=x, alpha=alpha, color=c, linewidth=0) if legend: - dimensions = list(set(item[0] for item in persistence)) + dimensions = set(item[0] for item in persistence) axes.legend( handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right", ) axes.set_title("Persistence barcode", fontsize=fontsize) + axes.set_yticks([]) + axes.invert_yaxis() # Ends plot on infinity value and starts a little bit before min_birth if len(x) != 0: - axes.axis([axis_start, infinity, 0, len(x)]) + axes.set_xlim((axis_start, infinity)) return axes except ImportError as import_error: diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py index de5844f9..7dc83817 100644 --- a/src/python/gudhi/point_cloud/knn.py +++ b/src/python/gudhi/point_cloud/knn.py @@ -314,7 +314,9 @@ class KNearestNeighbors: return None if self.params["implementation"] == "ckdtree": - qargs = {key: val for key, val in self.params.items() if key in {"p", "eps", "n_jobs"}} + qargs = {key: val for key, val in self.params.items() if key in {"p", "eps"}} + # SciPy renamed n_jobs to workers + qargs["workers"] = self.params.get("workers") or self.params.get("n_jobs") or 1 distances, neighbors = self.kdtree.query(X, k=self.k, **qargs) if k == 1: # SciPy decided to squeeze the last dimension for k=1 diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py index 69ff5e1e..a169aee8 100644 --- a/src/python/gudhi/representations/vector_methods.py +++ b/src/python/gudhi/representations/vector_methods.py @@ -85,7 +85,7 @@ class PersistenceImage(BaseEstimator, TransformerMixin): Xfit.append(image.flatten()[np.newaxis,:]) - Xfit = np.concatenate(Xfit,0) + Xfit = np.concatenate(Xfit, 0) return Xfit @@ -123,6 +123,15 @@ def _automatic_sample_range(sample_range, X, y): pass return sample_range + +def _trim_on_edges(x, are_endpoints_nan): + if are_endpoints_nan[0]: + x = x[1:] + if are_endpoints_nan[1]: + x = x[:-1] + return x + + class Landscape(BaseEstimator, TransformerMixin): """ This is a class for computing persistence landscapes from a list of persistence diagrams. A persistence landscape is a collection of 1D piecewise-linear functions computed from the rank function associated to the persistence diagram. These piecewise-linear functions are then sampled evenly on a given range and the corresponding vectors of samples are concatenated and returned. See http://jmlr.org/papers/v16/bubenik15a.html for more details. @@ -149,6 +158,8 @@ class Landscape(BaseEstimator, TransformerMixin): y (n x 1 array): persistence diagram labels (unused). """ self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + self.im_range = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) + self.im_range = _trim_on_edges(self.im_range, self.nan_in_range) return self def transform(self, X): @@ -161,53 +172,26 @@ class Landscape(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of diagrams) x (number of samples = **num_landscapes** x **resolution**): output persistence landscapes. """ - num_diag, Xfit = len(X), [] - x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) - step_x = x_values[1] - x_values[0] - - for i in range(num_diag): - - diagram, num_pts_in_diag = X[i], X[i].shape[0] - ls = np.zeros([self.num_landscapes, self.new_resolution]) + Xfit = [] + x_values = self.im_range + for diag in X: + midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2. + tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0) + n_points = diag.shape[0] + # Complete the array with zeros to get the right number of landscapes + if self.num_landscapes > n_points: + tent_functions = np.concatenate( + [tent_functions, np.zeros((tent_functions.shape[0], self.num_landscapes-n_points))], + axis=1 + ) + tent_functions.partition(tent_functions.shape[1]-self.num_landscapes, axis=1) + landscapes = np.sort(tent_functions[:, -self.num_landscapes:], axis=1)[:, ::-1].T - events = [] - for j in range(self.new_resolution): - events.append([]) + landscapes = np.sqrt(2) * np.ravel(landscapes) + Xfit.append(landscapes) - for j in range(num_pts_in_diag): - [px,py] = diagram[j,:2] - min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - - if min_idx < self.new_resolution and max_idx > 0: - - landscape_value = self.sample_range[0] + min_idx * step_x - px - for k in range(min_idx, mid_idx): - events[k].append(landscape_value) - landscape_value += step_x - - landscape_value = py - self.sample_range[0] - mid_idx * step_x - for k in range(mid_idx, max_idx): - events[k].append(landscape_value) - landscape_value -= step_x - - for j in range(self.new_resolution): - events[j].sort(reverse=True) - for k in range( min(self.num_landscapes, len(events[j])) ): - ls[k,j] = events[j][k] - - if self.nan_in_range[0]: - ls = ls[:,1:] - if self.nan_in_range[1]: - ls = ls[:,:-1] - ls = np.sqrt(2)*np.reshape(ls,[1,-1]) - Xfit.append(ls) - - Xfit = np.concatenate(Xfit,0) - - return Xfit + return np.stack(Xfit, axis=0) def __call__(self, diag): """ @@ -219,7 +203,7 @@ class Landscape(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of samples = **num_landscapes** x **resolution**): output persistence landscape. """ - return self.fit_transform([diag])[0,:] + return self.fit_transform([diag])[0, :] class Silhouette(BaseEstimator, TransformerMixin): """ @@ -235,6 +219,8 @@ class Silhouette(BaseEstimator, TransformerMixin): sample_range ([double, double]): minimum and maximum for the weighted average domain, of the form [x_min, x_max] (default [numpy.nan, numpy.nan]). It is the interval on which samples will be drawn evenly. If one of the values is numpy.nan, it can be computed from the persistence diagrams with the fit() method. """ self.weight, self.resolution, self.sample_range = weight, resolution, sample_range + self.nan_in_range = np.isnan(np.array(self.sample_range)) + self.new_resolution = self.resolution + self.nan_in_range.sum() def fit(self, X, y=None): """ @@ -245,6 +231,8 @@ class Silhouette(BaseEstimator, TransformerMixin): y (n x 1 array): persistence diagram labels (unused). """ self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + self.im_range = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) + self.im_range = _trim_on_edges(self.im_range, self.nan_in_range) return self def transform(self, X): @@ -257,44 +245,19 @@ class Silhouette(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of diagrams) x (**resolution**): output persistence silhouettes. """ - num_diag, Xfit = len(X), [] - x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution) - step_x = x_values[1] - x_values[0] - - for i in range(num_diag): - - diagram, num_pts_in_diag = X[i], X[i].shape[0] + Xfit = [] + x_values = self.im_range - sh, weights = np.zeros(self.resolution), np.zeros(num_pts_in_diag) - for j in range(num_pts_in_diag): - weights[j] = self.weight(diagram[j,:]) + for diag in X: + midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2. + weights = np.array([self.weight(pt) for pt in diag]) total_weight = np.sum(weights) - for j in range(num_pts_in_diag): - - [px,py] = diagram[j,:2] - weight = weights[j] / total_weight - min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) - mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) - max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) - - if min_idx < self.resolution and max_idx > 0: - - silhouette_value = self.sample_range[0] + min_idx * step_x - px - for k in range(min_idx, mid_idx): - sh[k] += weight * silhouette_value - silhouette_value += step_x - - silhouette_value = py - self.sample_range[0] - mid_idx * step_x - for k in range(mid_idx, max_idx): - sh[k] += weight * silhouette_value - silhouette_value -= step_x - - Xfit.append(np.reshape(np.sqrt(2) * sh, [1,-1])) - - Xfit = np.concatenate(Xfit, 0) + tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0) + silhouette = np.sum(weights[None, :] / total_weight * tent_functions, axis=1) + Xfit.append(silhouette * np.sqrt(2)) - return Xfit + return np.stack(Xfit, axis=0) def __call__(self, diag): """ diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py index 4a455bb6..58caab21 100755 --- a/src/python/test/test_representations.py +++ b/src/python/test/test_representations.py @@ -187,3 +187,67 @@ def test_kernel_empty_diagrams(): # PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.)(empty_diag, empty_diag) # PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1., kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag) + +def test_silhouette_permutation_invariance(): + dgm = _n_diags(1)[0] + dgm_permuted = dgm[np.random.permutation(dgm.shape[0]).astype(int)] + random_resolution = random.randint(50, 100) * 10 + slt = Silhouette(resolution=random_resolution, weight=pow(2)) + + assert np.all(np.isclose(slt(dgm), slt(dgm_permuted))) + + +def test_silhouette_multiplication_invariance(): + dgm = _n_diags(1)[0] + n_repetitions = np.random.randint(2, high=10) + dgm_augmented = np.repeat(dgm, repeats=n_repetitions, axis=0) + + random_resolution = random.randint(50, 100) * 10 + slt = Silhouette(resolution=random_resolution, weight=pow(2)) + assert np.all(np.isclose(slt(dgm), slt(dgm_augmented))) + + +def test_silhouette_numeric(): + dgm = np.array([[2., 3.], [5., 6.]]) + slt = Silhouette(resolution=9, weight=pow(1), sample_range=[2., 6.]) + #slt.fit([dgm]) + # x_values = array([2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.]) + + expected_silhouette = np.array([0., 0.5, 0., 0., 0., 0., 0., 0.5, 0.])/np.sqrt(2) + output_silhouette = slt(dgm) + assert np.all(np.isclose(output_silhouette, expected_silhouette)) + + +def test_landscape_small_persistence_invariance(): + dgm = np.array([[2., 6.], [2., 5.], [3., 7.]]) + small_persistence_pts = np.random.rand(10, 2) + small_persistence_pts[:, 1] += small_persistence_pts[:, 0] + small_persistence_pts += np.min(dgm) + dgm_augmented = np.concatenate([dgm, small_persistence_pts], axis=0) + + lds = Landscape(num_landscapes=2, resolution=5) + lds_dgm, lds_dgm_augmented = lds(dgm), lds(dgm_augmented) + + assert np.all(np.isclose(lds_dgm, lds_dgm_augmented)) + + +def test_landscape_numeric(): + dgm = np.array([[2., 6.], [3., 5.]]) + lds_ref = np.array([ + 0., 0.5, 1., 1.5, 2., 1.5, 1., 0.5, 0., # tent of [2, 6] + 0., 0., 0., 0.5, 1., 0.5, 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., + ]) + lds_ref *= np.sqrt(2) + lds = Landscape(num_landscapes=4, resolution=9, sample_range=[2., 6.]) + lds_dgm = lds(dgm) + assert np.all(np.isclose(lds_dgm, lds_ref)) + + +def test_landscape_nan_range(): + dgm = np.array([[2., 6.], [3., 5.]]) + lds = Landscape(num_landscapes=2, resolution=9, sample_range=[np.nan, 6.]) + lds_dgm = lds(dgm) + assert (lds.sample_range[0] == 2) & (lds.sample_range[1] == 6) + assert lds.new_resolution == 10 |