diff options
author | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-03-25 12:38:56 +0000 |
---|---|---|
committer | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2015-03-25 12:38:56 +0000 |
commit | 9280807119434e1ca4dd90267233ba01dfb0eeb6 (patch) | |
tree | 41bd747e1ca55b52a0a4d21e088ec7297aae813f | |
parent | 1ec3e1ae82dab1109ce871f690df716fb05c6a16 (diff) | |
parent | 87032aea5bf14683d964ecd30bd9d744da975e1b (diff) |
Merged latest trunk to witness
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/witness@505 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: a90ff18c8991da76cf66daec3dc8f68c049cfe57
36 files changed, 4268 insertions, 2593 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index f3b29994..fabba412 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -69,6 +69,8 @@ else() message(STATUS "boost include dirs:" ${Boost_INCLUDE_DIRS}) include_directories(src/common/include/) + include_directories(src/Alpha_shapes/include/) + include_directories(src/Bottleneck/include/) include_directories(src/Contraction/include/) include_directories(src/Hasse_complex/include/) include_directories(src/Persistent_cohomology/include/) @@ -86,6 +88,13 @@ else() add_subdirectory(src/Hasse_complex/example) add_subdirectory(src/Witness_complex/test) add_subdirectory(src/Witness_complex/example) + add_subdirectory(src/Alpha_shapes/example) + add_subdirectory(src/Alpha_shapes/test) + add_subdirectory(src/Bottleneck/example) + add_subdirectory(src/Bottleneck/test) + + # data points generator + add_subdirectory(data/points/generator) endif() diff --git a/data/points/generator/CMakeLists.txt b/data/points/generator/CMakeLists.txt new file mode 100644 index 00000000..0f2674c4 --- /dev/null +++ b/data/points/generator/CMakeLists.txt @@ -0,0 +1,26 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIPointsGenerator) + +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + add_executable ( hypergenerator hypergenerator.cpp ) + add_test(hypergenerator_on_sphere_3000_10_5.0 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 3000 10 5.0) + add_test(hypergenerator_on_sphere_10000_3 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 10000 3) + add_test(hypergenerator_in_sphere_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 7000 12 10.8) + add_test(hypergenerator_in_sphere_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 50000 2) + # on cube is not available in CGAL + add_test(hypergenerator_in_cube_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 7000 12 10.8) + add_test(hypergenerator_in_cube_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 50000 3) + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/data/points/generator/README b/data/points/generator/README new file mode 100644 index 00000000..e6b28eb0 --- /dev/null +++ b/data/points/generator/README @@ -0,0 +1,26 @@ +To build the example, run in a Terminal: + +cd /path-to-gudhi/ +cmake . +cd /path-to-data-generator/ +make + + +Example of use : + +*** Hyper sphere|cube generator + +./hypergenerator on sphere onSphere.off 1000 3 15.2 + + => generates a onSphere.off file with 1000 points randomized on a sphere of dimension 3 and radius 15.2 + +./hypergenerator in sphere inSphere.off 100 2 + + => generates a inSphere.off file with 100 points randomized in a sphere of dimension 2 (circle) and radius 1.0 (default) + +./hypergenerator in cube inCube.off 10000 3 5.8 + + => generates a inCube.off file with 10000 points randomized in a cube of dimension 3 and side 5.8 + +!! Warning: hypegenerator on cube is not available !! + diff --git a/data/points/generator/hypergenerator.cpp b/data/points/generator/hypergenerator.cpp new file mode 100644 index 00000000..f4ea6b07 --- /dev/null +++ b/data/points/generator/hypergenerator.cpp @@ -0,0 +1,129 @@ +/* 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) 2014 INRIA Saclay (France) + * + * 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 <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <fstream> // for std::ofstream + + +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K; +typedef K::Point_d Point; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " in|on sphere|cube off_file_name points_number[integer > 0] dimension[integer > 1] radius[double > 0.0 | default = 1.0]" << std::endl; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + // program args management + if ((argc != 6) && (argc != 7)) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + int points_number = 0; + int returnedScanValue = sscanf(argv[4], "%d", &points_number); + if ((returnedScanValue == EOF) || (points_number <= 0)) { + std::cerr << "Error: " << argv[4] << " is not correct" << std::endl; + usage(argv[0]); + } + + int dimension = 0; + returnedScanValue = sscanf(argv[5], "%d", &dimension); + if ((returnedScanValue == EOF) || (dimension <= 0)) { + std::cerr << "Error: " << argv[5] << " is not correct" << std::endl; + usage(argv[0]); + } + + double radius = 1.0; + if (argc == 7) { + returnedScanValue = sscanf(argv[6], "%lf", &radius); + if ((returnedScanValue == EOF) || (radius <= 0.0)) { + std::cerr << "Error: " << argv[6] << " is not correct" << std::endl; + usage(argv[0]); + } + } + + bool in = false; + if (strcmp(argv[1], "in") == 0) { + in = true; + } else if (strcmp(argv[1], "on") != 0) { + std::cerr << "Error: " << argv[1] << " is not correct" << std::endl; + usage(argv[0]); + } + + bool sphere = false; + if (memcmp(argv[2], "sphere", sizeof("sphere")) == 0) { + sphere = true; + } else if (memcmp(argv[2], "cube", sizeof("cube")) != 0) { + std::cerr << "Error: " << argv[2] << " is not correct" << std::endl; + usage(argv[0]); + } + + std::ofstream diagram_out(argv[3]); + diagram_out << "OFF" << std::endl; + diagram_out << points_number << " 0 0" << std::endl; + + if (diagram_out.is_open()) { + // Instanciate a random point generator + CGAL::Random rng(0); + // Generate "points_number" random points in a vector + std::vector<Point> points; + if (in) { + if (sphere) { + CGAL::Random_points_in_ball_d<Point> rand_it(dimension, radius, rng); + CGAL::cpp11::copy_n(rand_it, points_number, std::back_inserter(points)); + } else { + CGAL::Random_points_in_cube_d<Point> rand_it(dimension, radius, rng); + CGAL::cpp11::copy_n(rand_it, points_number, std::back_inserter(points)); + } + } else { // means "on" + if (sphere) { + CGAL::Random_points_on_sphere_d<Point> rand_it(dimension, radius, rng); + CGAL::cpp11::copy_n(rand_it, points_number, std::back_inserter(points)); + } else { + std::cerr << "Sorry: on cube is not available" << std::endl; + usage(argv[0]); + } + } + + for (auto thePoint : points) { + int i = 0; + for (;i < dimension - 1; i++) { + diagram_out << thePoint[i] << " "; + } + diagram_out << thePoint[i] << std::endl; // last point + Carriage Return + } + } else { + std::cerr << "Error: " << argv[3] << " cannot be opened" << std::endl; + usage(argv[0]); + } + + return 0; +} + diff --git a/src/Alpha_shapes/example/CMakeLists.txt b/src/Alpha_shapes/example/CMakeLists.txt new file mode 100644 index 00000000..fb94ca05 --- /dev/null +++ b/src/Alpha_shapes/example/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIAlphaShapesExample) + +# need CGAL 4.6 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + add_executable ( dtoffrw Delaunay_triangulation_off_rw.cpp ) + target_link_libraries(dtoffrw ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3) + + #add_definitions(-DDEBUG_TRACES) + add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp ) + target_link_libraries(stfromdt ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp new file mode 100644 index 00000000..03f6440d --- /dev/null +++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp @@ -0,0 +1,105 @@ +/* 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) 2014 INRIA Saclay (France) + * + * 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/>. + */ + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> +#include <vector> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K; +typedef CGAL::Delaunay_triangulation<K> T; +// The triangulation uses the default instanciation of the +// TriangulationDataStructure template parameter + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if (argc != 3) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + int dimension = 0; + int returnedScanValue = sscanf(argv[2], "%d", &dimension); + if ((returnedScanValue == EOF) || (dimension <= 0)) { + std::cerr << "Error: " << argv[2] << " is not correct" << std::endl; + usage(argv[0]); + } + + T dt(dimension); + std::string offFileName(argv[1]); + Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(offFileName, dt, true, true); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offFileName << std::endl; + exit(-1); // ----- >> + } + + std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl; + std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl; + std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl; + + // Points list + /*for (T::Vertex_iterator vit = dt.vertices_begin(); vit != dt.vertices_end(); ++vit) { + for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { + std::cout << *Coord << " "; + } + std::cout << std::endl; + } + std::cout << std::endl;*/ + + int i = 0, j = 0; + typedef T::Full_cell_iterator Full_cell_iterator; + typedef T::Facet Facet; + + for (Full_cell_iterator cit = dt.full_cells_begin(); cit != dt.full_cells_end(); ++cit) { + if (!dt.is_infinite(cit)) { + j++; + continue; + } + Facet fct(cit, cit->index(dt.infinite_vertex())); + i++; + } + std::cout << "There are " << i << " facets on the convex hull." << std::endl; + std::cout << "There are " << j << " facets not on the convex hull." << std::endl; + + + std::string offOutputFile("out.off"); + Gudhi::alphashapes::Delaunay_triangulation_off_writer<T> off_writer(offOutputFile, dt); + + return 0; +}
\ No newline at end of file diff --git a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp new file mode 100644 index 00000000..3a17b519 --- /dev/null +++ b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp @@ -0,0 +1,103 @@ +/* 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) 2014 INRIA Saclay (France) + * + * 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/>. + */ + +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_shapes.h" + +// to construct a simplex_tree from Delaunay_triangulation +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Simplex_tree.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K; +typedef CGAL::Delaunay_triangulation<K> T; +// The triangulation uses the default instanciation of the +// TriangulationDataStructure template parameter + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if (argc != 3) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl; + usage(argv[0]); + } + + int dimension = 0; + int returnedScanValue = sscanf(argv[2], "%d", &dimension); + if ((returnedScanValue == EOF) || (dimension <= 0)) { + std::cerr << "Error: " << argv[2] << " is not correct" << std::endl; + usage(argv[0]); + } + + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a Delaunay triangulation + // + // ---------------------------------------------------------------------------- + T dt(dimension); + std::string off_file_name(argv[1]); + + Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, false, false); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_name << std::endl; + exit(-1); // ----- >> + } + + std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl; + std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl; + std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl; + + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt); + //std::cout << alpha_shapes_from_dt << std::endl; + + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a OFF file + // + // ---------------------------------------------------------------------------- + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, dimension); + //std::cout << alpha_shapes_from_file << std::endl; + + std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl; + std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl; + std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl; + std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl; + + return 0; +}
\ No newline at end of file diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h new file mode 100644 index 00000000..b8efdb4d --- /dev/null +++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h @@ -0,0 +1,200 @@ +/* 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) 2015 INRIA Saclay (France) + * + * 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 SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ +#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ + +// to construct a Delaunay_triangulation from a OFF file +#include <gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h> + +// to construct a simplex_tree from Delaunay_triangulation +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> + +#include <stdio.h> +#include <stdlib.h> + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <string> + +namespace Gudhi { + +namespace alphashapes { + +/** \defgroup alpha_shapes Alpha shapes in dimension N + * + <DT>Implementations:</DT> + Alpha shapes in dimension N are a subset of Delaunay Triangulation in dimension N. + + + * \author Vincent Rouvreau + * \version 1.0 + * \date 2015 + * \copyright GNU General Public License v3. + * @{ + */ + +/** + * \brief Alpha shapes data structure. + * + * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation + * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$. + * + * Details may be found in \cite boissonnatmariasimplextreealgorithmica. + * + * \implements FilteredComplex + * + */ +class Alpha_shapes { + private: + // From Simplex_tree + /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/ + typedef std::vector<Vertex_handle> typeVectorVertex; + + // From CGAL + /** \brief Kernel for the Delaunay_triangulation. + * Dimension can be set dynamically. + */ + typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel; + /** \brief Delaunay_triangulation type required to create an alpha-shape. + */ + typedef CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation; + + private: + /** \brief Upper bound on the simplex tree of the simplicial complex.*/ + Gudhi::Simplex_tree<> _st; + + public: + + Alpha_shapes(std::string off_file_name, int dimension) { + Delaunay_triangulation dt(dimension); + Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation> + off_reader(off_file_name, dt, true, true); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_name << std::endl; + exit(-1); // ----- >> + } +#ifdef DEBUG_TRACES + std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl; + std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl; + std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl; +#endif // DEBUG_TRACES + init<Delaunay_triangulation>(dt); + } + + template<typename T> + Alpha_shapes(T triangulation) { + init<T>(triangulation); + } + + ~Alpha_shapes() { } + + private: + + template<typename T> + void init(T triangulation) { + _st.set_dimension(triangulation.maximal_dimension()); + _st.set_filtration(0.0); + // triangulation points list + for (auto vit = triangulation.finite_vertices_begin(); + vit != triangulation.finite_vertices_end(); ++vit) { + typeVectorVertex vertexVector; + Vertex_handle vertexHdl = std::distance(triangulation.finite_vertices_begin(), vit); + vertexVector.push_back(vertexHdl); + + // Insert each point in the simplex tree + _st.insert_simplex(vertexVector, 0.0); + +#ifdef DEBUG_TRACES + std::cout << "P" << vertexHdl << ":"; + for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { + std::cout << *Coord << " "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + } + // triangulation finite full cells list + for (auto cit = triangulation.finite_full_cells_begin(); + cit != triangulation.finite_full_cells_end(); ++cit) { + typeVectorVertex vertexVector; + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + // Vertex handle is distance - 1 + Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1; + vertexVector.push_back(vertexHdl); + } + // Insert each point in the simplex tree + _st.insert_simplex_and_subfaces(vertexVector, 0.0); + +#ifdef DEBUG_TRACES + std::cout << "C" << std::distance(triangulation.finite_full_cells_begin(), cit) << ":"; + for (auto value : vertexVector) { + std::cout << value << ' '; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + } + } + + public: + + /** \brief Returns the number of vertices in the complex. */ + size_t num_vertices() { + return _st.num_vertices(); + } + + /** \brief Returns the number of simplices in the complex. + * + * Does not count the empty simplex. */ + const unsigned int& num_simplices() const { + return _st.num_simplices(); + } + + /** \brief Returns an upper bound on the dimension of the simplicial complex. */ + int dimension() { + return _st.dimension(); + } + + /** \brief Returns an upper bound of the filtration values of the simplices. */ + Filtration_value filtration() { + return _st.filtration(); + } + + friend std::ostream& operator<<(std::ostream& os, const Alpha_shapes& alpha_shape) { + // TODO: Program terminated with signal SIGABRT, Aborted - Maybe because of copy constructor + Gudhi::Simplex_tree<> st = alpha_shape._st; + os << st << std::endl; + return os; + } +}; + +} // namespace alphashapes + +} // namespace Gudhi + +#endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_ diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h new file mode 100644 index 00000000..693b393e --- /dev/null +++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h @@ -0,0 +1,213 @@ +/* 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) 2015 INRIA Saclay (France) + * + * 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 SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ +#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ + +#include <string> +#include <vector> +#include <fstream> + +#include "gudhi/Off_reader.h" + +namespace Gudhi { + +namespace alphashapes { + +/** + *@brief Off reader visitor with flag that can be passed to Off_reader to read a Delaunay_triangulation_complex. + */ +template<typename Complex> +class Delaunay_triangulation_off_flag_visitor_reader { + Complex& complex_; + typedef typename Complex::Point Point; + + const bool load_only_points_; + + public: + explicit Delaunay_triangulation_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { +#ifdef DEBUG_TRACES + std::cout << "init" << std::endl; +#endif // DEBUG_TRACES + } + + void point(const std::vector<double>& point) { +#ifdef DEBUG_TRACES + std::cout << "p "; + for (auto coordinate: point) { + std::cout << coordinate << " | "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + complex_.insert(Point(point.size(), point.begin(), point.end())); + } + + void maximal_face(const std::vector<int>& face) { + // For alpha shapes, only points are read + } + + void done() { +#ifdef DEBUG_TRACES + std::cout << "done" << std::endl; +#endif // DEBUG_TRACES + } +}; + +/** + *@brief Off reader visitor that can be passed to Off_reader to read a Delaunay_triangulation_complex. + */ +template<typename Complex> +class Delaunay_triangulation_off_visitor_reader { + Complex& complex_; + // typedef typename Complex::Vertex_handle Vertex_handle; + // typedef typename Complex::Simplex_handle Simplex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + std::vector<Point> points_; + // std::vector<Simplex_handle> maximal_faces_; + + public: + explicit Delaunay_triangulation_off_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { +#ifdef DEBUG_TRACES + std::cout << "init - " << num_vertices << std::endl; +#endif // DEBUG_TRACES + // maximal_faces_.reserve(num_faces); + points_.reserve(num_vertices); + } + + void point(const std::vector<double>& point) { +#ifdef DEBUG_TRACES + std::cout << "p "; + for (auto coordinate: point) { + std::cout << coordinate << " | "; + } + std::cout << std::endl; +#endif // DEBUG_TRACES + points_.emplace_back(Point(point.size(), point.begin(), point.end())); + } + + void maximal_face(const std::vector<int>& face) { + // For alpha shapes, only points are read + } + + void done() { + complex_.insert(points_.begin(), points_.end()); +#ifdef DEBUG_TRACES + std::cout << "done" << std::endl; +#endif // DEBUG_TRACES + } +}; + +/** + *@brief Class that allows to load a Delaunay_triangulation_complex from an off file. + */ +template<typename Complex> +class Delaunay_triangulation_off_reader { + public: + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false, + bool is_flag = false) : valid_(false) { + std::ifstream stream(name_file); + if (stream.is_open()) { + if (is_flag) { + // For alpha shapes, only points are read + Delaunay_triangulation_off_flag_visitor_reader<Complex> off_visitor(read_complex, true); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } else { + // For alpha shapes, only points are read + Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex, true); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + } + } + + /** + * return true if reading did not meet problems. + */ + bool is_valid() const { + return valid_; + } + + private: + bool valid_; +}; + +template<typename Complex> +class Delaunay_triangulation_off_writer { + public: + /** + * name_file : file where the off will be written + * save_complex : complex that be outputted in the file + * for now only save triangles. + */ + Delaunay_triangulation_off_writer(const std::string & name_file, const Complex& save_complex) { + std::ofstream stream(name_file); + if (stream.is_open()) { + // OFF header + stream << "OFF" << std::endl; + // no endl on next line - don't know why... + stream << save_complex.number_of_vertices() << " " << save_complex.number_of_finite_full_cells() << " 0"; + + // Points list + for (auto vit = save_complex.vertices_begin(); vit != save_complex.vertices_end(); ++vit) { + for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) { + stream << *Coord << " "; + } + stream << std::endl; + } + + // Finite cells list + for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) { + stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; // Dimension + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + auto vertexHdl = *vit; + // auto vertexHdl = std::distance(save_complex.vertices_begin(), *vit) - 1; + // stream << std::distance(save_complex.vertices_begin(), *(vit)) - 1 << " "; + } + stream << std::endl; + } + stream.close(); + } else { + std::cerr << "could not open file " << name_file << std::endl; + } + } +}; + +} // namespace alphashapes + +} // namespace Gudhi + +#endif // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_ diff --git a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp new file mode 100644 index 00000000..a90704b6 --- /dev/null +++ b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp @@ -0,0 +1,126 @@ +/* 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) 2015 INRIA Saclay (France) + * + * 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/>. + */ + +#define BOOST_TEST_MODULE alpha_shapes +#include <boost/test/included/unit_test.hpp> +#include <boost/system/error_code.hpp> +#include <boost/chrono/thread_clock.hpp> +// to construct a Delaunay_triangulation from a OFF file +#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h" +#include "gudhi/Alpha_shapes.h" + +// to construct a simplex_tree from Delaunay_triangulation +#include "gudhi/graph_simplicial_complex.h" +#include "gudhi/Simplex_tree.h" + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <iostream> +#include <iterator> + +#include <stdio.h> +#include <stdlib.h> +#include <string> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K; +typedef CGAL::Delaunay_triangulation<K> T; +// The triangulation uses the default instanciation of the +// TriangulationDataStructure template parameter + +BOOST_AUTO_TEST_CASE( OFF_file ) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("S4_100.off"); + std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl; + + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, 4); + + const int DIMENSION = 4; + std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.dimension() == DIMENSION); + + const double FILTRATION = 0.0; + std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.filtration() == FILTRATION); + + const int NUMBER_OF_VERTICES = 100; + std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.num_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_SIMPLICES = 6779; + std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl; + BOOST_CHECK(alpha_shapes_from_file.num_simplices() == NUMBER_OF_SIMPLICES); + +} + +BOOST_AUTO_TEST_CASE( Delaunay_triangulation ) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-shape from a Delauny triangulation + // + // ---------------------------------------------------------------------------- + T dt(8); + std::string off_file_name("S8_10.off"); + std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl; + + Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, true, true); + std::cout << "off_reader.is_valid()=" << off_reader.is_valid() << std::endl; + BOOST_CHECK(off_reader.is_valid()); + + const int NUMBER_OF_VERTICES = 10; + std::cout << "dt.number_of_vertices()=" << dt.number_of_vertices() << std::endl; + BOOST_CHECK(dt.number_of_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_FULL_CELLS = 30; + std::cout << "dt.number_of_full_cells()=" << dt.number_of_full_cells() << std::endl; + BOOST_CHECK(dt.number_of_full_cells() == NUMBER_OF_FULL_CELLS); + + const int NUMBER_OF_FINITE_FULL_CELLS = 6; + std::cout << "dt.number_of_finite_full_cells()=" << dt.number_of_finite_full_cells() << std::endl; + BOOST_CHECK(dt.number_of_finite_full_cells() == NUMBER_OF_FINITE_FULL_CELLS); + + Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt); + + const int DIMENSION = 8; + std::cout << "alpha_shapes_from_dt.dimension()=" << alpha_shapes_from_dt.dimension() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.dimension() == DIMENSION); + + const double FILTRATION = 0.0; + std::cout << "alpha_shapes_from_dt.filtration()=" << alpha_shapes_from_dt.filtration() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.filtration() == FILTRATION); + + std::cout << "alpha_shapes_from_dt.num_vertices()=" << alpha_shapes_from_dt.num_vertices() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.num_vertices() == NUMBER_OF_VERTICES); + + const int NUMBER_OF_SIMPLICES = 997; + std::cout << "alpha_shapes_from_dt.num_simplices()=" << alpha_shapes_from_dt.num_simplices() << std::endl; + BOOST_CHECK(alpha_shapes_from_dt.num_simplices() == NUMBER_OF_SIMPLICES); +} + diff --git a/src/Alpha_shapes/test/CMakeLists.txt b/src/Alpha_shapes/test/CMakeLists.txt new file mode 100644 index 00000000..a48c1a8f --- /dev/null +++ b/src/Alpha_shapes/test/CMakeLists.txt @@ -0,0 +1,31 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIAlphaShapesTest) + +# need CGAL 4.6 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + include_directories (BEFORE "../../include") + + add_definitions(-DDEBUG_TRACES) + add_executable ( AlphaShapesUnitTest Alpha_shapes_unit_test.cpp ) + target_link_libraries(AlphaShapesUnitTest ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_test(AlphaShapesUnitTest ${CMAKE_CURRENT_BINARY_DIR}/AlphaShapesUnitTest) + + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") + endif () +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Alpha_shapes/include/gudhi") diff --git a/src/Alpha_shapes/test/README b/src/Alpha_shapes/test/README new file mode 100644 index 00000000..244a2b84 --- /dev/null +++ b/src/Alpha_shapes/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./AlphaShapesUnitTest --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Alpha_shapes/test/S4_100.off b/src/Alpha_shapes/test/S4_100.off new file mode 100644 index 00000000..0a5dc58c --- /dev/null +++ b/src/Alpha_shapes/test/S4_100.off @@ -0,0 +1,102 @@ +OFF +100 0 0 +0.562921 -0.735261 -0.256472 0.277007 +-0.803733 -0.0527915 -0.315125 0.501918 +-0.24946 -0.354982 -0.410773 -0.801887 +0.916381 -0.0512295 0.371049 0.141223 +0.182222 0.940836 -0.171362 -0.228599 +-0.787145 -0.129213 0.568102 -0.202402 +0.0187866 -0.22093 -0.832882 -0.507095 +-0.4702 -0.533814 0.353776 0.607286 +-0.159798 -0.504771 0.263586 0.806346 +0.295546 0.162541 0.452931 0.825279 +0.242043 -0.107437 -0.913612 -0.308521 +0.875759 -0.113035 -0.469189 -0.0114505 +0.547877 -0.762247 -0.256972 -0.229729 +-0.172302 0.521057 0.412013 -0.727363 +-0.724729 -0.0574074 -0.0290602 0.686023 +0.700434 -0.102636 0.687285 0.162779 +-0.681386 0.0946893 0.610047 0.393178 +-0.847553 -0.357132 0.383743 0.0827718 +0.72297 -0.161631 -0.608517 0.284424 +0.757394 0.141549 0.196065 -0.606528 +0.78094 0.00901997 0.434536 0.448586 +0.14166 -0.619339 0.614589 -0.467582 +0.473105 -0.537832 -0.0103746 0.697711 +-0.208004 0.536218 0.818027 0.00605288 +0.743694 -0.628926 0.188072 0.126488 +-0.462228 -0.278147 0.35514 -0.76345 +-0.17361 0.249211 0.758567 -0.57648 +0.416958 -0.254924 -0.576373 -0.654946 +-0.590751 -0.286089 -0.424896 0.623402 +-0.639538 -0.739693 -0.203745 0.0482932 +0.0731787 0.132121 0.864022 0.480266 +-0.149644 -0.164724 -0.249746 -0.94239 +-0.348592 0.0120379 -0.928656 -0.126239 +0.395328 -0.54513 0.149976 -0.723917 +-0.974164 0.14707 0.157191 0.068302 +-0.166425 0.119943 -0.627304 -0.751269 +0.031947 -0.358518 -0.708301 -0.607251 +0.93781 -0.155368 -0.30951 0.0240237 +0.276094 -0.753155 -0.597088 -0.00387459 +-0.642876 -0.200551 -0.263517 -0.690687 +0.178711 0.604987 -0.262989 0.729993 +-0.520347 0.497922 -0.676144 0.155374 +-0.703999 0.500219 -0.484381 0.139789 +-0.131013 0.835735 0.506779 0.166004 +-0.536116 -0.566557 0.229226 0.582279 +-0.334105 0.158252 0.926091 0.0754059 +-0.0362677 0.296076 0.897108 0.325915 +-0.57486 0.798575 0.15324 -0.0912754 +0.498602 0.0186805 0.72824 0.469801 +-0.960329 0.0473356 0.261005 -0.0860505 +0.899134 -0.381392 -0.214508 0.00921711 +0.570576 0.567224 0.393019 -0.445237 +-0.761763 -0.614589 -0.0546476 -0.197513 +0.188584 0.289531 0.174031 0.922129 +-0.458506 -0.583876 0.639297 -0.2004 +0.785343 -0.21571 0.0794082 -0.574804 +0.0819036 0.65961 -0.247426 0.704973 +0.573125 0.49706 0.373026 0.534145 +-0.513286 -0.626226 0.208535 -0.548536 +0.460558 0.468686 0.507832 -0.55707 +0.716158 -0.488201 0.388209 -0.313164 +0.881074 0.152441 0.380128 -0.236589 +0.885793 0.0386389 0.161009 -0.433537 +-0.365162 0.298384 0.292846 0.831784 +0.364934 0.632269 -0.197205 -0.654346 +-0.31469 -0.429991 0.665304 -0.522923 +-0.734198 0.462914 -0.135691 -0.477756 +-0.422885 0.674444 -0.364143 -0.483419 +0.829218 -0.154622 -0.381147 0.378439 +0.887881 0.310479 -0.109528 0.321363 +-0.354398 -0.693974 0.456019 -0.429941 +-0.492045 -0.160008 0.044387 0.854587 +0.0595532 0.158421 0.412577 -0.895062 +-0.211441 0.491794 -0.153521 0.83058 +-0.33558 -0.504711 0.353831 -0.71236 +-0.735211 -0.197714 0.525626 0.379593 +0.465818 -0.424245 0.769469 -0.104627 +-0.641071 -0.286339 -0.704442 -0.103923 +-0.00446569 0.0249849 -0.194417 -0.980591 +-0.610081 -0.252448 0.176698 -0.729966 +-0.0859217 -0.154471 0.715027 0.676382 +0.091315 0.0723382 -0.855023 -0.505337 +0.165362 0.200983 -0.428242 -0.865373 +-0.587465 0.303019 -0.152442 0.734729 +0.454946 -0.319828 0.437063 -0.706902 +-0.384368 0.277509 0.879225 -0.0470385 +0.523335 -0.330233 -0.208592 0.757335 +0.895086 0.0448492 0.268089 -0.353466 +-0.0272491 -0.567336 -0.72254 -0.39411 +-0.0745014 -0.121818 -0.882466 0.448179 +0.382304 -0.240135 0.851109 -0.267941 +-0.418057 -0.852847 -0.3128 0.00606452 +-0.554046 0.304237 0.272381 -0.725453 +0.155115 -0.0894732 -0.245017 -0.952838 +0.114459 -0.130722 0.953669 0.245614 +0.0913002 -0.462466 0.244433 0.847374 +-0.198849 0.0785111 0.131441 -0.967997 +-0.303154 -0.686484 0.639333 0.167604 +0.521455 0.256835 -0.0584503 -0.811606 +-0.109787 0.870544 0.161523 0.451676 diff --git a/src/Alpha_shapes/test/S8_10.off b/src/Alpha_shapes/test/S8_10.off new file mode 100644 index 00000000..1d67e10f --- /dev/null +++ b/src/Alpha_shapes/test/S8_10.off @@ -0,0 +1,12 @@ +OFF +10 0 0 +0.440036 -0.574754 -0.200485 0.216537 -0.501251 -0.0329236 -0.196529 0.313023 +-0.129367 -0.184089 -0.213021 -0.415848 0.783529 -0.0438025 0.317256 0.120749 +0.132429 0.683748 -0.124536 -0.166133 -0.540695 -0.0887576 0.390234 -0.139031 +0.0137399 -0.161581 -0.609142 -0.370872 -0.320669 -0.364053 0.24127 0.41416 +-0.115313 -0.36425 0.190208 0.581871 0.204605 0.112527 0.313562 0.571337 +0.168272 -0.0746917 -0.635156 -0.214488 0.629498 -0.0812499 -0.337255 -0.00823068 +0.369896 -0.514626 -0.173493 -0.1551 -0.127105 0.384377 0.303936 -0.536566 +-0.49013 -0.0388242 -0.0196532 0.463953 0.515962 -0.0756047 0.506276 0.119908 +-0.434258 0.060347 0.388793 0.250579 -0.653127 -0.275207 0.295714 0.0637842 +0.596172 -0.133284 -0.501793 0.234541 0.428452 0.0800735 0.110912 -0.343109 diff --git a/src/Bottleneck/concept/Persistence_diagram.h b/src/Bottleneck/concept/Persistence_diagram.h new file mode 100644 index 00000000..eaaf8bc5 --- /dev/null +++ b/src/Bottleneck/concept/Persistence_diagram.h @@ -0,0 +1,7 @@ +typedef typename std::pair<double,double> Diagram_point; + +struct Persistence_Diagram +{ + const_iterator<Diagram_point> cbegin() const; + const_iterator<Diagram_point> cend() const; +}; diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt new file mode 100644 index 00000000..2ff009c4 --- /dev/null +++ b/src/Bottleneck/example/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBottleneckExample) + +add_executable ( RandomDiagrams random_diagrams.cpp ) +add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams) diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck/example/random_diagrams.cpp new file mode 100644 index 00000000..71f152a6 --- /dev/null +++ b/src/Bottleneck/example/random_diagrams.cpp @@ -0,0 +1,40 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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/Graph_matching.h" +#include <iostream> + +using namespace Gudhi::bottleneck; + +int main() { + int n = 100; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + int a = rand() % n; + v1.emplace_back(a, a + rand() % (n - a)); + int b = rand() % n; + v2.emplace_back(b, b + rand() % (n - b)); + } + // v1 and v2 are persistence diagrams containing each 100 randoms points. + double b = bottleneck_distance(v1, v2, 0); + std::cout << b << std::endl; +} diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..ea47e1d5 --- /dev/null +++ b/src/Bottleneck/include/gudhi/Graph_matching.h @@ -0,0 +1,197 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ + +#include <deque> +#include <list> +#include <vector> + +#include "gudhi/Layered_neighbors_finder.h" + +namespace Gudhi { + +namespace bottleneck { + +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); + +class Graph_matching { + public: + Graph_matching(const Persistence_diagrams_graph& g); + Graph_matching& operator=(const Graph_matching& m); + bool perfect() const; + bool multi_augment(); + void set_r(double r); + + private: + const Persistence_diagrams_graph& g; + double r; + std::vector<int> v_to_u; + std::list<int> unmatched_in_u; + + Layered_neighbors_finder* layering() const; + bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth); + void update(std::deque<int>& path); +}; + +Graph_matching::Graph_matching(const Persistence_diagrams_graph& g) + : g(g), r(0), v_to_u(g.size()), unmatched_in_u() { + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.emplace_back(u_point_index); +} + +Graph_matching& Graph_matching::operator=(const Graph_matching& m) { + r = m.r; + v_to_u = m.v_to_u; + unmatched_in_u = m.unmatched_in_u; + return *this; +} + +inline bool Graph_matching::perfect() const { + return unmatched_in_u.empty(); +} + +inline bool Graph_matching::multi_augment() { + if (perfect()) + return false; + Layered_neighbors_finder* layered_nf = layering(); + double rn = sqrt(g.size()); + int nblmax = layered_nf->vlayers_number()*2 + 1; + // verification of a necessary criterion + if ((unmatched_in_u.size() > rn && nblmax > rn) || nblmax == 0) + return false; + bool successful = false; + std::list<int>* tries = new std::list<int>(unmatched_in_u); + for (auto it = tries->cbegin(); it != tries->cend(); it++) + successful = successful || augment(layered_nf, *it, nblmax); + delete tries; + delete layered_nf; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +Layered_neighbors_finder* Graph_matching::layering() const { + bool end = false; + int layer = 0; + std::list<int> u_vertices(unmatched_in_u); + std::list<int> v_vertices; + Neighbors_finder nf(g, r); + Layered_neighbors_finder* layered_nf = new Layered_neighbors_finder(g, r); + for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + nf.add(v_point_index); + while (!u_vertices.empty()) { + for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) { + std::list<int>* u_succ = nf.pull_all_near(*it); + for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) { + layered_nf->add(*it, layer); + v_vertices.emplace_back(*it); + } + delete u_succ; + } + u_vertices.clear(); + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) { + if (v_to_u.at(*it) == null_point_index()) + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + } + if (end) + return layered_nf; + v_vertices.clear(); + layer++; + } + return layered_nf; +} + +bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) { + std::deque<int> path; + path.emplace_back(u_start_index); + // start is a point from U + do { + if (static_cast<int>(path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + int w = path.back(); + path.emplace_back(layered_nf->pull_near(w, path.size() / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + path.pop_back(); + update(path); + return true; +} + +void Graph_matching::update(std::deque<int>& path) { + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + int tmp = *it; + ++it; + v_to_u[*it] = tmp; + } +} + +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) { + Persistence_diagrams_graph g(diag1, diag2, e); + std::vector<double>* sd = g.sorted_distances(); + int idmin = 0; + int idmax = sd->size() - 1; + double alpha = pow(sd->size(), 0.25); + Graph_matching m(g); + Graph_matching biggest_unperfect = m; + while (idmin != idmax) { + int pas = static_cast<int>((idmax - idmin) / alpha); + m.set_r(sd->at(idmin + pas)); + while (m.multi_augment()) {} + if (m.perfect()) { + idmax = idmin + pas; + m = biggest_unperfect; + } else { + biggest_unperfect = m; + idmin = idmin + pas + 1; + } + } + double b = sd->at(idmin); + delete sd; + return b; +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_ diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h new file mode 100644 index 00000000..de36e00b --- /dev/null +++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h @@ -0,0 +1,74 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ + +#include <vector> + +#include "Neighbors_finder.h" + +// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given +// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index. +// A neighbor returned is automatically removed. + +namespace Gudhi { + +namespace bottleneck { + +class Layered_neighbors_finder { + public: + Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r); + void add(int v_point_index, int vlayer); + int pull_near(int u_point_index, int vlayer); + int vlayers_number() const; + + private: + const Persistence_diagrams_graph& g; + const double r; + std::vector<Neighbors_finder> neighbors_finder; +}; + +Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) : + g(g), r(r), neighbors_finder() { } + +inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(Neighbors_finder(g, r)); + neighbors_finder.at(vlayer).add(v_point_index); +} + +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer).pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return neighbors_finder.size(); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..98256571 --- /dev/null +++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h @@ -0,0 +1,96 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ + +#include <unordered_set> +#include <list> + +#include "gudhi/Planar_neighbors_finder.h" + +namespace Gudhi { + +namespace bottleneck { + +// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V +// in the persistence diagrams graph. +// V's points have to be added manually using their index. A neighbor returned is automatically removed. + +class Neighbors_finder { + public: + Neighbors_finder(const Persistence_diagrams_graph& g, double r); + void add(int v_point_index); + int pull_near(int u_point_index); + std::list<int>* pull_all_near(int u_point_index); + + private: + const Persistence_diagrams_graph& g; + const double r; + Planar_neighbors_finder planar_neighbors_f; + std::unordered_set<int> projections_f; +}; + +Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) : + g(g), r(r), planar_neighbors_f(g, r), projections_f() { } + +inline void Neighbors_finder::add(int v_point_index) { + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + planar_neighbors_f.add(v_point_index); +} + +inline int Neighbors_finder::pull_near(int u_point_index) { + int v_challenger = g.corresponding_point_in_v(u_point_index); + if (planar_neighbors_f.contains(v_challenger) && g.distance(u_point_index, v_challenger) < r) { + planar_neighbors_f.remove(v_challenger); + return v_challenger; + } + if (g.on_the_u_diagonal(u_point_index)) { + auto it = projections_f.cbegin(); + if (it != projections_f.cend()) { + int tmp = *it; + projections_f.erase(it); + return tmp; + } + } else { + return planar_neighbors_f.pull_near(u_point_index); + } + return null_point_index(); +} + +inline std::list<int>* Neighbors_finder::pull_all_near(int u_point_index) { + std::list<int>* all_pull = planar_neighbors_f.pull_all_near(u_point_index); + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull->emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h new file mode 100644 index 00000000..7e278209 --- /dev/null +++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h @@ -0,0 +1,147 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ + +#include <vector> +#include <set> +#include <cmath> +#include <utility> // for pair<> +#include <algorithm> // for max + +namespace Gudhi { + +namespace bottleneck { + +// Diagram_point is the type of the persistence diagram's points +typedef typename std::pair<double, double> Diagram_point; + +// Return the used index for encoding none of the points +int null_point_index(); + +// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and +// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity. + +class Persistence_diagrams_graph { + public: + // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams. + // They have to have an iterator over points, which have to have fields first (for birth) and second (for death). + template<typename Persistence_diagram1, typename Persistence_diagram2> + Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.); + Persistence_diagrams_graph(); + bool on_the_u_diagonal(int u_point_index) const; + bool on_the_v_diagonal(int v_point_index) const; + int corresponding_point_in_u(int v_point_index) const; + int corresponding_point_in_v(int u_point_index) const; + double distance(int u_point_index, int v_point_index) const; + int size() const; + std::vector<double>* sorted_distances(); + + private: + std::vector<Diagram_point> u; + std::vector<Diagram_point> v; + Diagram_point get_u_point(int u_point_index) const; + Diagram_point get_v_point(int v_point_index) const; +}; + +inline int null_point_index() { + return -1; +} + +template<typename Persistence_diagram1, typename Persistence_diagram2> +Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) + : u(), v() { + for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) + if (it->second - it->first > e) + u.emplace_back(*it); + for (auto it = diag2.cbegin(); it != diag2.cend(); ++it) + if (it->second - it->first > e) + v.emplace_back(*it); + if (u.size() < v.size()) + swap(u, v); +} + +Persistence_diagrams_graph::Persistence_diagrams_graph::Persistence_diagrams_graph() + : u(), v() { } + +inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const { + return u_point_index >= static_cast<int> (u.size()); +} + +inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const { + return v_point_index >= static_cast<int> (v.size()); +} + +inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); +} + +inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); +} + +inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const { + // could be optimized for the case where one point is the projection of the other + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0; + Diagram_point p_u = get_u_point(u_point_index); + Diagram_point p_v = get_v_point(v_point_index); + return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second)); +} + +inline int Persistence_diagrams_graph::size() const { + return static_cast<int> (u.size() + v.size()); +} + +inline std::vector<double>* Persistence_diagrams_graph::sorted_distances() { + // could be optimized + std::set<double> sorted_distances; + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + sorted_distances.emplace(distance(u_point_index, v_point_index)); + return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()); +} + +inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const { + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Diagram_point projector = v.at(corresponding_point_in_v(u_point_index)); + double x = (projector.first + projector.second) / 2; + return Diagram_point(x, x); +} + +inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const { + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Diagram_point projector = u.at(corresponding_point_in_u(v_point_index)); + double x = (projector.first + projector.second) / 2; + return Diagram_point(x, x); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h new file mode 100644 index 00000000..4af672e4 --- /dev/null +++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h @@ -0,0 +1,119 @@ +/* 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): Francois Godi + * + * Copyright (C) 2015 INRIA Saclay (France) + * + * 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 SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ +#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ + +#include <list> +#include <iostream> +#include <set> + +#include "Persistence_diagrams_graph.h" + +namespace Gudhi { + +namespace bottleneck { + +// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the +// planar distance. +// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also +// remove points manually using their index. + +class Abstract_planar_neighbors_finder { + public: + Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r); + virtual ~Abstract_planar_neighbors_finder() = 0; + virtual void add(int v_point_index) = 0; + virtual void remove(int v_point_index) = 0; + virtual bool contains(int v_point_index) const = 0; + virtual int pull_near(int u_point_index) = 0; + virtual std::list<int>* pull_all_near(int u_point_index); + + protected: + const Persistence_diagrams_graph& g; + const double r; +}; + + +// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder + +class Naive_pnf : public Abstract_planar_neighbors_finder { + public: + Naive_pnf(const Persistence_diagrams_graph& g, double r); + void add(int v_point_index); + void remove(int v_point_index); + bool contains(int v_point_index) const; + int pull_near(int u_point_index); + + private: + std::set<int> candidates; +}; + + +// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation +typedef Naive_pnf Planar_neighbors_finder; + +Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) : + g(g), r(r) { } + +inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { } + +inline std::list<int>* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) { + std::list<int>* all_pull = new std::list<int>(); + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull->emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) : + Abstract_planar_neighbors_finder(g, r), candidates() { } + +inline void Naive_pnf::add(int v_point_index) { + candidates.emplace(v_point_index); +} + +inline void Naive_pnf::remove(int v_point_index) { + candidates.erase(v_point_index); +} + +inline bool Naive_pnf::contains(int v_point_index) const { + return (candidates.count(v_point_index) > 0); +} + +inline int Naive_pnf::pull_near(int u_point_index) { + for (auto it = candidates.begin(); it != candidates.end(); ++it) + if (g.distance(u_point_index, *it) <= r) { + int tmp = *it; + candidates.erase(it); + return tmp; + } + return null_point_index(); +} + +} // namespace bottleneck + +} // namespace Gudhi + +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt new file mode 100644 index 00000000..7044372e --- /dev/null +++ b/src/Bottleneck/test/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIBottleneckUnitTest) + +if(NOT MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage") +endif() + +add_executable ( BottleneckUnitTest bottleneck_unit_test.cpp ) +target_link_libraries(BottleneckUnitTest ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests +add_test(BottleneckUnitTest ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUnitTest) + +if (LCOV_PATH) + # Lcov code coverage of unitary test + add_test(src/Bottleneck/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Bottleneck) +endif() + +cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Bottleneck/include/gudhi") diff --git a/src/Bottleneck/test/README b/src/Bottleneck/test/README new file mode 100644 index 00000000..0e7b8673 --- /dev/null +++ b/src/Bottleneck/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./BottleneckUnitTest --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..068b8690 --- /dev/null +++ b/src/Bottleneck/test/bottleneck_unit_test.cpp @@ -0,0 +1,26 @@ +#define BOOST_TEST_MODULE bottleneck test + +#include <boost/test/included/unit_test.hpp> + +#include "gudhi/Graph_matching.h" +#include <iostream> + +using namespace Gudhi::bottleneck; + +BOOST_AUTO_TEST_CASE(random_diagrams) { + int n = 100; + // Random construction + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + int a = rand() % n; + v1.emplace_back(a, a + rand() % (n - a)); + int b = rand() % n; + v2.emplace_back(b, b + rand() % (n - b)); + } + // v1 and v2 are persistence diagrams containing each 100 randoms points. + double b = bottleneck_distance(v1, v2, 0); + // + std::cout << b << std::endl; + const double EXPECTED_DISTANCE = 98.5; + BOOST_CHECK(b == EXPECTED_DISTANCE); +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 04c4369c..362ff7a8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -38,5 +38,7 @@ else() add_subdirectory(example/Contraction) add_subdirectory(example/Hasse_complex) add_subdirectort(example/Witness_complex) + add_subdirectory(example/Alpha_shapes) + add_subdirectory(example/Bottleneck) endif() diff --git a/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp b/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp index c7cff843..1de43ab7 100644 --- a/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp +++ b/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA Saclay (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp index 4ecad752..bde224f1 100644 --- a/src/Simplex_tree/example/simple_simplex_tree.cpp +++ b/src/Simplex_tree/example/simple_simplex_tree.cpp @@ -31,273 +31,285 @@ typedef std::vector< Vertex_handle > typeVectorVertex; typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex; typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool; -int main (int argc, char * const argv[]) -{ - const Filtration_value FIRST_FILTRATION_VALUE = 0.1; - const Filtration_value SECOND_FILTRATION_VALUE = 0.2; - const Filtration_value THIRD_FILTRATION_VALUE = 0.3; - const Filtration_value FOURTH_FILTRATION_VALUE = 0.4; - Vertex_handle FIRST_VERTEX_HANDLE = (Vertex_handle)0; - Vertex_handle SECOND_VERTEX_HANDLE = (Vertex_handle)1; - Vertex_handle THIRD_VERTEX_HANDLE = (Vertex_handle)2; - Vertex_handle FOURTH_VERTEX_HANDLE = (Vertex_handle)3; - - // TEST OF INSERTION - std::cout << "********************************************************************" << std::endl; - std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl; - //Construct the Simplex Tree - Simplex_tree<> simplexTree; - - /* Simplex to be inserted: */ - /* 1 */ - /* o */ - /* /X\ */ - /* o---o---o */ - /* 2 0 3 */ - - // ++ FIRST - std::cout << " * INSERT 0" << std::endl; - typeVectorVertex firstSimplexVector; - firstSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); - typePairSimplexBool returnValue = - simplexTree.insert_simplex ( firstSimplex.first, firstSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + 0 INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - 0 NOT INSERTED" << std::endl; - } - - // ++ SECOND - std::cout << " * INSERT 1" << std::endl; - typeVectorVertex secondSimplexVector; - secondSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( secondSimplex.first, secondSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + 1 INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - 1 NOT INSERTED" << std::endl; - } - - // ++ THIRD - std::cout << " * INSERT (0,1)" << std::endl; - typeVectorVertex thirdSimplexVector; - thirdSimplexVector.push_back(FIRST_VERTEX_HANDLE); - thirdSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( thirdSimplex.first, thirdSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (0,1) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (0,1) NOT INSERTED" << std::endl; - } - - // ++ FOURTH - std::cout << " * INSERT 2" << std::endl; - typeVectorVertex fourthSimplexVector; - fourthSimplexVector.push_back(THIRD_VERTEX_HANDLE); - typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( fourthSimplex.first, fourthSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + 2 INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - 2 NOT INSERTED" << std::endl; - } - - // ++ FIFTH - std::cout << " * INSERT (2,0)" << std::endl; - typeVectorVertex fifthSimplexVector; - fifthSimplexVector.push_back(THIRD_VERTEX_HANDLE); - fifthSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( fifthSimplex.first, fifthSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (2,0) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (2,0) NOT INSERTED" << std::endl; - } - - // ++ SIXTH - std::cout << " * INSERT (2,1)" << std::endl; - typeVectorVertex sixthSimplexVector; - sixthSimplexVector.push_back(THIRD_VERTEX_HANDLE); - sixthSimplexVector.push_back(SECOND_VERTEX_HANDLE); - typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( sixthSimplex.first, sixthSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (2,1) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (2,1) NOT INSERTED" << std::endl; - } - - // ++ SEVENTH - std::cout << " * INSERT (2,1,0)" << std::endl; - typeVectorVertex seventhSimplexVector; - seventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); - seventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); - seventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( seventhSimplex.first, seventhSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (2,1,0) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (2,1,0) NOT INSERTED" << std::endl; - } - - // ++ EIGHTH - std::cout << " * INSERT 3" << std::endl; - typeVectorVertex eighthSimplexVector; - eighthSimplexVector.push_back(FOURTH_VERTEX_HANDLE); - typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( eighthSimplex.first, eighthSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + 3 INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - 3 NOT INSERTED" << std::endl; - } - - // ++ NINETH - std::cout << " * INSERT (3,0)" << std::endl; - typeVectorVertex ninethSimplexVector; - ninethSimplexVector.push_back(FOURTH_VERTEX_HANDLE); - ninethSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( ninethSimplex.first, ninethSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (3,0) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (3,0) NOT INSERTED" << std::endl; - } - - // ++ TENTH - std::cout << " * INSERT 0 (already inserted)" << std::endl; - typeVectorVertex tenthSimplexVector; - tenthSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value - returnValue = - simplexTree.insert_simplex ( tenthSimplex.first, tenthSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + 0 INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - 0 NOT INSERTED" << std::endl; - } - - // ++ ELEVENTH - std::cout << " * INSERT (2,1,0) (already inserted)" << std::endl; - typeVectorVertex eleventhSimplexVector; - eleventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); - eleventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); - eleventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); - typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); - returnValue = - simplexTree.insert_simplex ( eleventhSimplex.first, eleventhSimplex.second ); - - if (returnValue.second == true) - { - std::cout << " + (2,1,0) INSERTED" << std::endl; - int nb_simplices = simplexTree.num_simplices() + 1; - simplexTree.set_num_simplices(nb_simplices); - } - else - { - std::cout << " - (2,1,0) NOT INSERTED" << std::endl; - } - - // ++ GENERAL VARIABLE SET - simplexTree.set_filtration(FOURTH_FILTRATION_VALUE); // Max filtration value - simplexTree.set_dimension(2); // Max dimension = 2 -> (2,1,0) - - std::cout << "********************************************************************" << std::endl; - // Display the Simplex_tree - Can not be done in the middle of 2 inserts - std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices" << std::endl; - std::cout << " - dimension " << simplexTree.dimension() << " - filtration " << simplexTree.filtration() << std::endl; - std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; - for( auto f_simplex : simplexTree.filtration_simplex_range() ) - { - std::cout << " " << "[" << simplexTree.filtration(f_simplex) << "] "; - for( auto vertex : simplexTree.simplex_vertex_range(f_simplex) ) - { - std::cout << (int)vertex << " "; - } - std::cout << std::endl; - } - // [0.1] 0 - // [0.1] 1 - // [0.1] 2 - // [0.1] 3 - // [0.2] 1 0 - // [0.2] 2 0 - // [0.2] 2 1 - // [0.2] 3 0 - // [0.3] 2 1 0 - return 0; +int main(int argc, char * const argv[]) { + const Filtration_value FIRST_FILTRATION_VALUE = 0.1; + const Filtration_value SECOND_FILTRATION_VALUE = 0.2; + const Filtration_value THIRD_FILTRATION_VALUE = 0.3; + const Filtration_value FOURTH_FILTRATION_VALUE = 0.4; + Vertex_handle FIRST_VERTEX_HANDLE = (Vertex_handle) 0; + Vertex_handle SECOND_VERTEX_HANDLE = (Vertex_handle) 1; + Vertex_handle THIRD_VERTEX_HANDLE = (Vertex_handle) 2; + Vertex_handle FOURTH_VERTEX_HANDLE = (Vertex_handle) 3; + + // TEST OF INSERTION + std::cout << "********************************************************************" << std::endl; + std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl; + //Construct the Simplex Tree + Simplex_tree<> simplexTree; + + /* Simplex to be inserted: */ + /* 1 */ + /* o */ + /* /X\ */ + /* o---o---o */ + /* 2 0 3 */ + + // ++ FIRST + std::cout << " * INSERT 0" << std::endl; + typeVectorVertex firstSimplexVector; + firstSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + typePairSimplexBool returnValue = + simplexTree.insert_simplex(firstSimplex.first, firstSimplex.second); + + if (returnValue.second == true) { + std::cout << " + 0 INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - 0 NOT INSERTED" << std::endl; + } + + // ++ SECOND + std::cout << " * INSERT 1" << std::endl; + typeVectorVertex secondSimplexVector; + secondSimplexVector.push_back(SECOND_VERTEX_HANDLE); + typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(secondSimplex.first, secondSimplex.second); + + if (returnValue.second == true) { + std::cout << " + 1 INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - 1 NOT INSERTED" << std::endl; + } + + // ++ THIRD + std::cout << " * INSERT (0,1)" << std::endl; + typeVectorVertex thirdSimplexVector; + thirdSimplexVector.push_back(FIRST_VERTEX_HANDLE); + thirdSimplexVector.push_back(SECOND_VERTEX_HANDLE); + typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(thirdSimplex.first, thirdSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (0,1) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (0,1) NOT INSERTED" << std::endl; + } + + // ++ FOURTH + std::cout << " * INSERT 2" << std::endl; + typeVectorVertex fourthSimplexVector; + fourthSimplexVector.push_back(THIRD_VERTEX_HANDLE); + typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(fourthSimplex.first, fourthSimplex.second); + + if (returnValue.second == true) { + std::cout << " + 2 INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - 2 NOT INSERTED" << std::endl; + } + + // ++ FIFTH + std::cout << " * INSERT (2,0)" << std::endl; + typeVectorVertex fifthSimplexVector; + fifthSimplexVector.push_back(THIRD_VERTEX_HANDLE); + fifthSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(fifthSimplex.first, fifthSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (2,0) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (2,0) NOT INSERTED" << std::endl; + } + + // ++ SIXTH + std::cout << " * INSERT (2,1)" << std::endl; + typeVectorVertex sixthSimplexVector; + sixthSimplexVector.push_back(THIRD_VERTEX_HANDLE); + sixthSimplexVector.push_back(SECOND_VERTEX_HANDLE); + typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(sixthSimplex.first, sixthSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (2,1) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (2,1) NOT INSERTED" << std::endl; + } + + // ++ SEVENTH + std::cout << " * INSERT (2,1,0)" << std::endl; + typeVectorVertex seventhSimplexVector; + seventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); + seventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); + seventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(seventhSimplex.first, seventhSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (2,1,0) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (2,1,0) NOT INSERTED" << std::endl; + } + + // ++ EIGHTH + std::cout << " * INSERT 3" << std::endl; + typeVectorVertex eighthSimplexVector; + eighthSimplexVector.push_back(FOURTH_VERTEX_HANDLE); + typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(eighthSimplex.first, eighthSimplex.second); + + if (returnValue.second == true) { + std::cout << " + 3 INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - 3 NOT INSERTED" << std::endl; + } + + // ++ NINETH + std::cout << " * INSERT (3,0)" << std::endl; + typeVectorVertex ninethSimplexVector; + ninethSimplexVector.push_back(FOURTH_VERTEX_HANDLE); + ninethSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(ninethSimplex.first, ninethSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (3,0) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (3,0) NOT INSERTED" << std::endl; + } + + // ++ TENTH + std::cout << " * INSERT 0 (already inserted)" << std::endl; + typeVectorVertex tenthSimplexVector; + tenthSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value + returnValue = + simplexTree.insert_simplex(tenthSimplex.first, tenthSimplex.second); + + if (returnValue.second == true) { + std::cout << " + 0 INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - 0 NOT INSERTED" << std::endl; + } + + // ++ ELEVENTH + std::cout << " * INSERT (2,1,0) (already inserted)" << std::endl; + typeVectorVertex eleventhSimplexVector; + eleventhSimplexVector.push_back(THIRD_VERTEX_HANDLE); + eleventhSimplexVector.push_back(SECOND_VERTEX_HANDLE); + eleventhSimplexVector.push_back(FIRST_VERTEX_HANDLE); + typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); + returnValue = + simplexTree.insert_simplex(eleventhSimplex.first, eleventhSimplex.second); + + if (returnValue.second == true) { + std::cout << " + (2,1,0) INSERTED" << std::endl; + int nb_simplices = simplexTree.num_simplices() + 1; + simplexTree.set_num_simplices(nb_simplices); + } else { + std::cout << " - (2,1,0) NOT INSERTED" << std::endl; + } + + // ++ GENERAL VARIABLE SET + simplexTree.set_filtration(FOURTH_FILTRATION_VALUE); // Max filtration value + simplexTree.set_dimension(2); // Max dimension = 2 -> (2,1,0) + + std::cout << "********************************************************************" << std::endl; + // Display the Simplex_tree - Can not be done in the middle of 2 inserts + std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices" << std::endl; + std::cout << " - dimension " << simplexTree.dimension() << " - filtration " << simplexTree.filtration() << std::endl; + std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:" << std::endl; + for (auto f_simplex : simplexTree.filtration_simplex_range()) { + std::cout << " " << "[" << simplexTree.filtration(f_simplex) << "] "; + for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) { + std::cout << (int) vertex << " "; + } + std::cout << std::endl; + } + // [0.1] 0 + // [0.1] 1 + // [0.1] 2 + // [0.1] 3 + // [0.2] 1 0 + // [0.2] 2 0 + // [0.2] 2 1 + // [0.2] 3 0 + // [0.3] 2 1 0 + + // ------------------------------------------------------------------------------------------------------------------ + // Find in the simplex_tree + // ------------------------------------------------------------------------------------------------------------------ + Simplex_tree<>::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector); + std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != simplexTree.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + + Vertex_handle UNKNOWN_VERTEX_HANDLE = (Vertex_handle) 15; + typeVectorVertex unknownSimplexVector; + unknownSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE); + simplexFound = simplexTree.find(unknownSimplexVector); + std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != simplexTree.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + + simplexFound = simplexTree.find(fifthSimplexVector); + std::cout << "**************IS THE SIMPLEX {2,0} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != simplexTree.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + + typeVectorVertex otherSimplexVector; + otherSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE); + otherSimplexVector.push_back(SECOND_VERTEX_HANDLE); + simplexFound = simplexTree.find(otherSimplexVector); + std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != simplexTree.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + + typeVectorVertex invSimplexVector; + invSimplexVector.push_back(SECOND_VERTEX_HANDLE); + invSimplexVector.push_back(THIRD_VERTEX_HANDLE); + invSimplexVector.push_back(FIRST_VERTEX_HANDLE); + simplexFound = simplexTree.find(invSimplexVector); + std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != simplexTree.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + return 0; } diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index e52fe1ae..f95679cb 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -397,7 +397,7 @@ class Simplex_tree { * <CODE>Vertex_handle</CODE>. */ template<class RandomAccessVertexRange> - Simplex_handle find(const RandomAccessVertexRange & s) { + Simplex_handle find(RandomAccessVertexRange & s) { if (s.begin() == s.end()) std::cerr << "Empty simplex \n"; diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp index eaec4881..c0cfced1 100644 --- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp +++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp @@ -512,6 +512,67 @@ BOOST_AUTO_TEST_CASE( NSimplexAndSubfaces_tree_insertion ) test_simplex_tree_contains(st,simplexPair4,2); // (1,0) is in position 2 test_simplex_tree_contains(st,simplexPair5,14); // (3,4,5) is in position 14 test_simplex_tree_contains(st,simplexPair6,26); // (7,6,1,0) is in position 26 + + // ------------------------------------------------------------------------------------------------------------------ + // Find in the simplex_tree + // ------------------------------------------------------------------------------------------------------------------ + typeVectorVertex simpleSimplexVector; + simpleSimplexVector.push_back(SECOND_VERTEX_HANDLE); + Simplex_tree<>::Simplex_handle simplexFound = st.find(simpleSimplexVector); + std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + + Vertex_handle UNKNOWN_VERTEX_HANDLE = (Vertex_handle) 15; + typeVectorVertex unknownSimplexVector; + unknownSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE); + simplexFound = st.find(unknownSimplexVector); + std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + // Check it is NOT found + BOOST_CHECK(simplexFound == st.null_simplex()); + + simplexFound = st.find(SimplexVector6); + std::cout << "**************IS THE SIMPLEX {0,1,6,7} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + + typeVectorVertex otherSimplexVector; + otherSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE); + otherSimplexVector.push_back(SECOND_VERTEX_HANDLE); + simplexFound = st.find(otherSimplexVector); + std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + // Check it is NOT found + BOOST_CHECK(simplexFound == st.null_simplex()); + + typeVectorVertex invSimplexVector; + invSimplexVector.push_back(SECOND_VERTEX_HANDLE); + invSimplexVector.push_back(THIRD_VERTEX_HANDLE); + invSimplexVector.push_back(FIRST_VERTEX_HANDLE); + simplexFound = st.find(invSimplexVector); + std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n"; + if (simplexFound != st.null_simplex()) + std::cout << "***+ YES IT IS!\n"; + else + std::cout << "***- NO IT ISN'T\n"; + // Check it is found + BOOST_CHECK(simplexFound != st.null_simplex()); + // Display the Simplex_tree - Can not be done in the middle of 2 inserts std::cout << "The complex contains " << st.num_simplices() << " simplices" << std::endl; std::cout << " - dimension " << st.dimension() << " - filtration " << st.filtration() << std::endl; diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h index c98b0b45..aaaab8b0 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h @@ -24,6 +24,7 @@ #include <string> #include <vector> +#include <map> #include "gudhi/Off_reader.h" @@ -36,38 +37,36 @@ namespace skbl { */ template<typename Complex> class Skeleton_blocker_off_flag_visitor_reader { - Complex& complex_; - typedef typename Complex::Vertex_handle Vertex_handle; - typedef typename Complex::Point Point; - - const bool load_only_points_; - -public: - explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false): - complex_(complex), - load_only_points_(load_only_points) {} - - - void init(int dim, int num_vertices, int num_faces, int num_edges) { - // todo do an assert to check that this number are correctly read - // todo reserve size for vector points - } - - - void point(const std::vector<double>& point) { - complex_.add_vertex(Point(point.begin(),point.end())); - } - - void maximal_face(const std::vector<int>& face) { - if (!load_only_points_) { - for (size_t i = 0; i < face.size(); ++i) - for (size_t j = i+1; j < face.size(); ++j) { - complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j])); - } - } - } - - void done() {} + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + + public: + explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { + // todo do an assert to check that this number are correctly read + // todo reserve size for vector points + } + + void point(const std::vector<double>& point) { + complex_.add_vertex(Point(point.begin(), point.end())); + } + + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + for (size_t i = 0; i < face.size(); ++i) + for (size_t j = i + 1; j < face.size(); ++j) { + complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j])); + } + } + } + + void done() { } }; /** @@ -75,137 +74,127 @@ public: */ template<typename Complex> class Skeleton_blocker_off_visitor_reader { - Complex& complex_; - typedef typename Complex::Vertex_handle Vertex_handle; - typedef typename Complex::Simplex_handle Simplex_handle; - typedef typename Complex::Point Point; - - const bool load_only_points_; - std::vector<Point> points_; - std::vector<Simplex_handle> maximal_faces_; - -public: - explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false): - complex_(complex), - load_only_points_(load_only_points) {} - - - void init(int dim, int num_vertices, int num_faces, int num_edges){ - maximal_faces_.reserve(num_faces); - points_.reserve(num_vertices); - } - - - void point(const std::vector<double>& point) { - points_.emplace_back(point.begin(),point.end()); - } - - void maximal_face(const std::vector<int>& face) { - if (!load_only_points_) { - Simplex_handle s; - for(auto x : face) - s.add_vertex(Vertex_handle(x)); - maximal_faces_.emplace_back(s); - } - } - - void done() { - complex_ = make_complex_from_top_faces( - maximal_faces_.begin(), - maximal_faces_.end(), - points_.begin(), - points_.end() - ); - } + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Simplex_handle Simplex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + std::vector<Point> points_; + std::vector<Simplex_handle> maximal_faces_; + + public: + explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { + maximal_faces_.reserve(num_faces); + points_.reserve(num_vertices); + } + + void point(const std::vector<double>& point) { + points_.emplace_back(point.begin(), point.end()); + } + + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + Simplex_handle s; + for (auto x : face) + s.add_vertex(Vertex_handle(x)); + maximal_faces_.emplace_back(s); + } + } + + void done() { + complex_ = make_complex_from_top_faces(maximal_faces_.begin(), maximal_faces_.end(), + points_.begin(), points_.end() ); + } }; - /** *@brief Class that allows to load a Skeleton_blocker_complex from an off file. */ template<typename Complex> class Skeleton_blocker_off_reader { -public: - /** - * name_file : file to read - * read_complex : complex that will receive the file content - * read_only_points : specify true if only the points must be read - */ - Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,bool is_flag = false):valid_(false) { - std::ifstream stream(name_file); - if (stream.is_open()) { - if(is_flag || read_only_points){ - Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points); - Off_reader off_reader(stream); - valid_ = off_reader.read(off_visitor); - } - else{ - Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points); - Off_reader off_reader(stream); - valid_ = off_reader.read(off_visitor); - } - } - } - - /** - * return true iff reading did not meet problems. - */ - bool is_valid() const { - return valid_; - } - -private: - bool valid_; + public: + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, + bool read_only_points = false, bool is_flag = false) : valid_(false) { + std::ifstream stream(name_file); + if (stream.is_open()) { + if (is_flag || read_only_points) { + Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } else { + Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + } + } + + /** + * return true iff reading did not meet problems. + */ + bool is_valid() const { + return valid_; + } + + private: + bool valid_; }; - template<typename Complex> -class Skeleton_blocker_off_writer{ -public: - /** - * name_file : file where the off will be written - * save_complex : complex that be outputted in the file - * for now only save triangles. - */ - Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex){ - typedef typename Complex::Vertex_handle Vertex_handle; - - std::ofstream stream(name_file); - if (stream.is_open()) { - stream<<"OFF\n"; - size_t num_triangles = std::distance(save_complex.triangle_range().begin(),save_complex.triangle_range().end()); - stream<<save_complex.num_vertices()<<" "<<num_triangles<< " 0 \n"; - - //in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance - //we compute a map from 0 2 5 7 to 0 1 2 3 - std::map<Vertex_handle,size_t> vertex_num; - size_t current_vertex=0; - - for(auto v : save_complex.vertex_range()){ - vertex_num[v]=current_vertex++; - const auto& pt(save_complex.point(v)); - for(auto x : pt) - stream<<x<<" "; - stream<<std::endl; - } - - for(const auto & t : save_complex.triangle_range()){ - stream<<"3 "; - for(auto x : t) - stream<<vertex_num[x]<<" "; - stream<<std::endl; - } - stream.close(); - } else std::cerr<<"could not open file "<<name_file<<std::endl; - } - +class Skeleton_blocker_off_writer { + public: + /** + * name_file : file where the off will be written + * save_complex : complex that be outputted in the file + * for now only save triangles. + */ + Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex) { + typedef typename Complex::Vertex_handle Vertex_handle; + + std::ofstream stream(name_file); + if (stream.is_open()) { + stream << "OFF\n"; + size_t num_triangles = std::distance(save_complex.triangle_range().begin(), save_complex.triangle_range().end()); + stream << save_complex.num_vertices() << " " << num_triangles << " 0 \n"; + + // in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance + // we compute a map from 0 2 5 7 to 0 1 2 3 + std::map<Vertex_handle, size_t> vertex_num; + size_t current_vertex = 0; + + for (auto v : save_complex.vertex_range()) { + vertex_num[v] = current_vertex++; + const auto& pt(save_complex.point(v)); + for (auto x : pt) + stream << x << " "; + stream << std::endl; + } + + for (const auto & t : save_complex.triangle_range()) { + stream << "3 "; + for (auto x : t) + stream << vertex_num[x] << " "; + stream << std::endl; + } + stream.close(); + } else { + std::cerr << "could not open file " << name_file << std::endl; + } + } }; - } // namespace skbl - } // namespace Gudhi - #endif // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h index 5b2c1f27..700830f2 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h @@ -63,1479 +63,1466 @@ namespace skbl { */ template<class SkeletonBlockerDS> class Skeleton_blocker_complex { - template<class ComplexType> friend class Vertex_iterator; - template<class ComplexType> friend class Neighbors_vertices_iterator; - template<class ComplexType> friend class Edge_iterator; - template<class ComplexType> friend class Edge_around_vertex_iterator; - - template<class ComplexType> friend class Skeleton_blocker_link_complex; - template<class ComplexType> friend class Skeleton_blocker_link_superior; - template<class ComplexType> friend class Skeleton_blocker_sub_complex; - -public: - /** - * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS - */ - typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; - - /** - * @brief The type of stored edge node, specified by the template SkeletonBlockerDS - */ - typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; - - typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; - - /** - * @brief The type of an handle to a vertex of the complex. - */ - typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; - typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; - - /** - * @brief A ordered set of integers that represents a simplex. - */ - typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; - typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; - - /** - * @brief Handle to a blocker of the complex. - */ - typedef Simplex_handle* Blocker_handle; - - typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; - typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; - -protected: - typedef typename boost::adjacency_list<boost::setS, // edges - boost::vecS, // vertices - boost::undirectedS, Graph_vertex, Graph_edge> Graph; - // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink - // (eg Link with greater vertices) - // that burdens simplex iteration / complex initialization via list of simplices. - // to avoid that, one should modify the graph by storing two lists of adjacency for every - // vertex, the one with superior and the one with lower vertices, that way there is - // no more extra cost for computation of SuperiorLink - typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; - typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; - -protected: - typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; - -public: - /** - * @brief Handle to an edge of the complex. - */ - typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; - -protected: - typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator; - typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator; - -protected: - int num_vertices_; - int num_blockers_; - - typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; - // typedef Visitor* Visitor_ptr; - Visitor* visitor; - - /** - * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. - * - * This quantity is updated when adding/removing edge. - * - * This is useful because the operation - * list.size() is done in linear time. - */ - std::vector<boost_vertex_handle> degree_; - Graph skeleton; /** 1-skeleton of the simplicial complex. */ - - /** Each vertex can access to the blockers passing through it. */ - BlockerMap blocker_map_; - -public: - ///////////////////////////////////////////////////////////////////////////// - /** @name Constructors, Destructors - */ - //@{ - - /** - *@brief constructs a simplicial complex with a given number of vertices and a visitor. - */ - explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) - : visitor(visitor_) { - clear(); - for (int i = 0; i < num_vertices_; ++i) { - add_vertex(); - } - } - -private: - // typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie; - typedef Trie<Simplex_handle> STrie; - -public: - /** - * @brief Constructor with a list of simplices. - * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). - */ - template<typename SimpleHandleOutputIterator> - Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL) - : num_vertices_(0), - num_blockers_(0), - visitor(visitor_) { - add_vertex_and_edges(simplex_begin,simplex_end); - - if (!is_flag_complex) - // need to compute blockers - add_blockers(simplex_begin,simplex_end); - } - -private: - template<typename SimpleHandleOutputIterator> - void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){ - std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; - // first pass to add vertices and edges - int num_vertex = -1; - for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { - if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex); - if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex()); - } - while (num_vertex-- >= 0) add_vertex(); - - for (const auto& e : edges) - add_edge(e.first, e.second); - } - - template<typename SimpleHandleOutputIterator> - void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){ - Tries<Simplex_handle> tries(num_vertices(),simplex_begin,simplex_end); - tries.init_next_dimension(); - auto simplices(tries.next_dimension_simplices()); - - while(!simplices.empty()){ - simplices = tries.next_dimension_simplices(); - for(auto& sigma : simplices){ - //common_positive_neighbors is the set of vertices u such that - //for all s in sigma, us is an edge and u>s - Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex())); - for(auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it/**/) - if(sigma_it != sigma.rbegin()) - common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it)); - - for(auto x : common_positive_neighbors){ - //first test that all edges sx are here for all s in sigma - bool all_edges_here = true; - for(auto s : sigma) - if(!contains_edge(x,s)){ - all_edges_here = false; - break; - } - if(!all_edges_here) continue; - - // all edges of sigma \cup x are here - //we have a blocker if all proper faces of sigma \cup x - // are in the complex and if sigma \cup x is not in the complex - //the first is equivalent at checking if blocks(sigma \cup x) is true - //as all blockers of lower dimension have already been computed - sigma.add_vertex(x); - if(!tries.contains(sigma)&&!blocks(sigma)) - add_blocker(sigma); - sigma.remove_vertex(x); - } - - } - } - } - - -public: - - - // We cannot use the default copy constructor since we need - // to make a copy of each of the blockers - - Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) { - visitor = NULL; - degree_ = copy.degree_; - skeleton = Graph(copy.skeleton); - num_vertices_ = copy.num_vertices_; - - num_blockers_ = 0; - // we copy the blockers - for (auto blocker : copy.const_blocker_range()) { - add_blocker(*blocker); - } - } - - /** - */ - Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) { - clear(); - visitor = NULL; - degree_ = copy.degree_; - skeleton = Graph(copy.skeleton); - num_vertices_ = copy.num_vertices_; - - num_blockers_ = 0; - // we copy the blockers - for (auto blocker : copy.const_blocker_range()) - add_blocker(*blocker); - return *this; - } - - - /** - * return true if both complexes have the same simplices. - */ - bool operator==(const Skeleton_blocker_complex& other) const{ - if(other.num_vertices() != num_vertices()) return false; - if(other.num_edges() != num_edges()) return false; - if(other.num_blockers() != num_blockers()) return false; - - for(auto v : vertex_range()) - if(!other.contains_vertex(v)) return false; - - for(auto e : edge_range()) - if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false; - - for(const auto b : const_blocker_range()) - if(!other.contains_blocker(*b)) return false; - - return true; - } - - bool operator!=(const Skeleton_blocker_complex& other) const{ - return !(*this == other); - } - - /** - * The destructor delete all blockers allocated. - */ - virtual ~Skeleton_blocker_complex() { - clear(); - } - - /** - * @details Clears the simplicial complex. After a call to this function, - * blockers are destroyed. The 1-skeleton and the set of blockers - * are both empty. - */ - virtual void clear() { - // xxx for now the responsabilty of freeing the visitor is for - // the user - visitor = NULL; - - degree_.clear(); - num_vertices_ = 0; - - remove_blockers(); - - skeleton.clear(); - } - - /** - *@brief allows to change the visitor. - */ - void set_visitor(Visitor* other_visitor) { - visitor = other_visitor; - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Vertices operations - */ - //@{ -public: - /** - * @brief Return a local Vertex_handle of a vertex given a global one. - * @remark Assume that the vertex is present in the complex. - */ - Vertex_handle operator[](Root_vertex_handle global) const { - auto local(get_address(global)); - assert(local); - return *local; - } - - /** - * @brief Return the vertex node associated to local Vertex_handle. - * @remark Assume that the vertex is present in the complex. - */ - Graph_vertex& operator[](Vertex_handle address) { - assert( - 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); - return skeleton[address.vertex]; - } - - /** - * @brief Return the vertex node associated to local Vertex_handle. - * @remark Assume that the vertex is present in the complex. - */ - const Graph_vertex& operator[](Vertex_handle address) const { - assert( - 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); - return skeleton[address.vertex]; - } - - /** - * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. - */ - Vertex_handle add_vertex() { - Vertex_handle address(boost::add_vertex(skeleton)); - num_vertices_++; - (*this)[address].activate(); - // safe since we now that we are in the root complex and the field 'address' and 'id' - // are identical for every vertices - (*this)[address].set_id(Root_vertex_handle(address.vertex)); - degree_.push_back(0); - if (visitor) - visitor->on_add_vertex(address); - return address; - } - - /** - * @brief Remove a vertex from the simplicial complex - * @remark It just deactivates the vertex with a boolean flag but does not - * remove it from vertices from complexity issues. - */ - void remove_vertex(Vertex_handle address) { - assert(contains_vertex(address)); - // We remove b - boost::clear_vertex(address.vertex, skeleton); - (*this)[address].deactivate(); - num_vertices_--; - degree_[address.vertex] = -1; - if (visitor) - visitor->on_remove_vertex(address); - } - - /** - */ - bool contains_vertex(Vertex_handle u) const { - if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton)) - return false; - return (*this)[u].is_active(); - } - - /** - */ - bool contains_vertex(Root_vertex_handle u) const { - boost::optional<Vertex_handle> address = get_address(u); - return address && (*this)[*address].is_active(); - } - - /** - * @return true iff the simplicial complex contains all vertices - * of simplex sigma - */ - bool contains_vertices(const Simplex_handle & sigma) const { - for (auto vertex : sigma) - if (!contains_vertex(vertex)) - return false; - return true; - } - - /** - * @brief Given an Id return the address of the vertex having this Id in the complex. - * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex. - */ - virtual boost::optional<Vertex_handle> get_address( - Root_vertex_handle id) const { - boost::optional<Vertex_handle> res; - if (id.vertex < boost::num_vertices(skeleton)) - res = Vertex_handle(id.vertex); // xxx - return res; - } - - /** - * return the id of a vertex of adress local present in the graph - */ - Root_vertex_handle get_id(Vertex_handle local) const { - assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); - return (*this)[local].get_id(); - } - - /** - * @brief Convert an address of a vertex of a complex to the address in - * the current complex. - * @details - * If the current complex is a sub (or sup) complex of 'other', it converts - * the address of a vertex v expressed in 'other' to the address of the vertex - * v in the current one. - * @remark this methods uses Root_vertex_handle to identify the vertex and - * assumes the vertex is present in the current complex. - */ - Vertex_handle convert_handle_from_another_complex( - const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const { - auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); - assert(vh_in_current_complex); - return *vh_in_current_complex; - } - - /** - * @brief return the graph degree of a vertex. - */ - int degree(Vertex_handle local) const { - assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); - return degree_[local.vertex]; - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Edges operations - */ - //@{ -public: - /** - * @brief return an edge handle if the two vertices forms - * an edge in the complex - */ - boost::optional<Edge_handle> operator[]( - const std::pair<Vertex_handle, Vertex_handle>& ab) const { - boost::optional<Edge_handle> res; - std::pair<Edge_handle, bool> edge_pair( - boost::edge(ab.first.vertex, ab.second.vertex, skeleton)); - if (edge_pair.second) - res = edge_pair.first; - return res; - } - - /** - * @brief returns the stored node associated to an edge - */ - Graph_edge& operator[](Edge_handle edge_handle) { - return skeleton[edge_handle]; - } - - /** - * @brief returns the stored node associated to an edge - */ - const Graph_edge& operator[](Edge_handle edge_handle) const { - return skeleton[edge_handle]; - } - - /** - * @brief returns the first vertex of an edge - * @details it assumes that the edge is present in the complex - */ - Vertex_handle first_vertex(Edge_handle edge_handle) const { - return static_cast<Vertex_handle> (source(edge_handle, skeleton)); - } - - /** - * @brief returns the first vertex of an edge - * @details it assumes that the edge is present in the complex - */ - Vertex_handle second_vertex(Edge_handle edge_handle) const { - return static_cast<Vertex_handle> (target(edge_handle, skeleton)); - } - - /** - * @brief returns the simplex made with the two vertices of the edge - * @details it assumes that the edge is present in the complex - - */ - Simplex_handle get_vertices(Edge_handle edge_handle) const { - auto edge((*this)[edge_handle]); - return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]); - } - - /** - * @brief Adds an edge between vertices a and b and all its cofaces. - */ - Edge_handle add_edge(Vertex_handle a, Vertex_handle b) { - assert(contains_vertex(a) && contains_vertex(b)); - assert(a != b); - - auto edge_handle((*this)[std::make_pair(a, b)]); - // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)]; - // Edge_handle edge_descr; - // bool edge_present = pair_descr_bool.second; - if (!edge_handle) { - edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first; - (*this)[*edge_handle].setId(get_id(a), get_id(b)); - degree_[a.vertex]++; - degree_[b.vertex]++; - if (visitor) - visitor->on_add_edge(a, b); - } - return *edge_handle; - } - - /** - * @brief Adds all edges and their cofaces of a simplex to the simplicial complex. - */ - void add_edges(const Simplex_handle & sigma) { - Simplex_handle_iterator i, j; - for (i = sigma.begin(); i != sigma.end(); ++i) - for (j = i, j++; j != sigma.end(); ++j) - add_edge(*i, *j); - } - - /** - * @brief Removes an edge from the simplicial complex and all its cofaces. - * @details returns the former Edge_handle representing the edge - */ - virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) { - bool found; - Edge_handle edge; - tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton); - if (found) { - if (visitor) - visitor->on_remove_edge(a, b); - boost::remove_edge(a.vertex, b.vertex, skeleton); - degree_[a.vertex]--; - degree_[b.vertex]--; - } - return edge; - } - - /** - * @brief Removes edge and its cofaces from the simplicial complex. - */ - void remove_edge(Edge_handle edge) { - assert(contains_vertex(first_vertex(edge))); - assert(contains_vertex(second_vertex(edge))); - remove_edge(first_vertex(edge), second_vertex(edge)); - } - - /** - * @brief The complex is reduced to its set of vertices. - * All the edges and blockers are removed. - */ - void keep_only_vertices() { - remove_blockers(); - - for (auto u : vertex_range()) { - while (this->degree(u) > 0) { - Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); - this->remove_edge(u, v); - } - } - } - - /** - * @return true iff the simplicial complex contains an edge between - * vertices a and b - */ - bool contains_edge(Vertex_handle a, Vertex_handle b) const { - // if (a.vertex<0 || b.vertex <0) return false; - return boost::edge(a.vertex, b.vertex, skeleton).second; - } - - /** - * @return true iff the simplicial complex contains all vertices - * and all edges of simplex sigma - */ - bool contains_edges(const Simplex_handle & sigma) const { - for (auto i = sigma.begin(); i != sigma.end(); ++i) { - if (!contains_vertex(*i)) - return false; - for (auto j = i; ++j != sigma.end();) { - if (!contains_edge(*i, *j)) - return false; - } - } - return true; - } - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Blockers operations - */ - //@{ - - /** - * @brief Adds the simplex to the set of blockers and - * returns a Blocker_handle toward it if was not present before and 0 otherwise. - */ - Blocker_handle add_blocker(const Simplex_handle& blocker) { - assert(blocker.dimension() > 1); - if (contains_blocker(blocker)) { - // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; - return 0; - } else { - if (visitor) - visitor->on_add_blocker(blocker); - Blocker_handle blocker_pt = new Simplex_handle(blocker); - num_blockers_++; - auto vertex = blocker_pt->begin(); - while (vertex != blocker_pt->end()) { - blocker_map_.insert(BlockerPair(*vertex, blocker_pt)); - ++vertex; - } - return blocker_pt; - } - } - -protected: - /** - * @brief Adds the simplex to the set of blockers - */ - void add_blocker(Blocker_handle blocker) { - if (contains_blocker(*blocker)) { - // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; - return; - } else { - if (visitor) - visitor->on_add_blocker(*blocker); - num_blockers_++; - auto vertex = blocker->begin(); - while (vertex != blocker->end()) { - blocker_map_.insert(BlockerPair(*vertex, blocker)); - ++vertex; - } - } - } - -protected: - /** - * Removes sigma from the blocker map of vertex v - */ - void remove_blocker(const Blocker_handle sigma, Vertex_handle v) { - Complex_blocker_around_vertex_iterator blocker; - for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end(); - ++blocker) { - if (*blocker == sigma) - break; - } - if (*blocker != sigma) { - std::cerr - << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; - assert(false); - } else { - blocker_map_.erase(blocker.current_position()); - } - } - -public: - /** - * @brief Removes the simplex from the set of blockers. - * @remark sigma has to belongs to the set of blockers - */ - void remove_blocker(const Blocker_handle sigma) { - for (auto vertex : *sigma) - remove_blocker(sigma, vertex); - num_blockers_--; - } - - /** - * @brief Remove all blockers, in other words, it expand the simplicial - * complex to the smallest flag complex that contains it. - */ - void remove_blockers() { - // Desallocate the blockers - while (!blocker_map_.empty()) { - delete_blocker(blocker_map_.begin()->second); - } - num_blockers_ = 0; - blocker_map_.clear(); - } - -protected: - /** - * Removes the simplex sigma from the set of blockers. - * sigma has to belongs to the set of blockers - * - * @remark contrarily to delete_blockers does not call the destructor - */ - void remove_blocker(const Simplex_handle& sigma) { - assert(contains_blocker(sigma)); - for (auto vertex : sigma) - remove_blocker(sigma, vertex); - num_blockers_--; - } - -public: - /** - * Removes the simplex s from the set of blockers - * and desallocate s. - */ - void delete_blocker(Blocker_handle sigma) { - if (visitor) - visitor->on_delete_blocker(sigma); - remove_blocker(sigma); - delete sigma; - } - - /** - * @return true iff s is a blocker of the simplicial complex - */ - bool contains_blocker(const Blocker_handle s) const { - if (s->dimension() < 2) - return false; - - Vertex_handle a = s->first_vertex(); - - for (const auto blocker : const_blocker_range(a)) { - if (s == *blocker) - return true; - } - return false; - } - - /** - * @return true iff s is a blocker of the simplicial complex - */ - bool contains_blocker(const Simplex_handle & s) const { - if (s.dimension() < 2) - return false; - - Vertex_handle a = s.first_vertex(); - - for (auto blocker : const_blocker_range(a)) { - if (s == *blocker) - return true; - } - return false; - } - -private: - /** - * @return true iff a blocker of the simplicial complex - * is a face of sigma. - */ - bool blocks(const Simplex_handle & sigma) const { - for(auto s : sigma) - for (auto blocker : const_blocker_range(s)) - if (sigma.contains(*blocker)) - return true; - return false; - - } - - //@} - -protected: - /** - * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$. - * If keep_only_superior is true then only vertices that are greater than v are added. - */ - virtual void add_neighbours(Vertex_handle v, Simplex_handle & n, - bool keep_only_superior = false) const { - boost_adjacency_iterator ai, ai_end; - for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; - ++ai) { - if (keep_only_superior) { - if (*ai > v.vertex) { - n.add_vertex(Vertex_handle(*ai)); - } - } else { - n.add_vertex(Vertex_handle(*ai)); - } - } - } - - /** - * @details Add to simplex res all vertices which are - * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. - * - * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. - * todo revoir - * - */ - virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res, - bool keep_only_superior = false) const { - res.clear(); - auto alpha_vertex = alpha.begin(); - add_neighbours(*alpha_vertex, res, keep_only_superior); - for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end(); - ++alpha_vertex) - keep_neighbours(*alpha_vertex, res, keep_only_superior); - } - - /** - * @details Remove from simplex n all vertices which are - * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$. - * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. - */ - virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res, - bool keep_only_superior = false) const { - Simplex_handle nv; - add_neighbours(v, nv, keep_only_superior); - res.intersection(nv); - } - - /** - * @details Remove from simplex all vertices which are - * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$. - * If 'keep_only_superior' is true then only vertices that are greater than v are added. - */ - virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res, - bool keep_only_superior = false) const { - Simplex_handle nv; - add_neighbours(v, nv, keep_only_superior); - res.difference(nv); - } - -public: - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex; - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(Vertex_handle v) const { - return Link_complex(*this, Simplex_handle(v)); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(Edge_handle edge) const { - return Link_complex(*this, edge); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Link_complex link(const Simplex_handle& simplex) const { - return Link_complex(*this, simplex); - } - - /** - * @brief Compute the local vertices of 's' in the current complex - * If one of them is not present in the complex then the return value is uninitialized. - * - * - */ - // xxx rename get_address et place un using dans sub_complex - - boost::optional<Simplex_handle> get_simplex_address( - const Root_simplex_handle& s) const { - boost::optional<Simplex_handle> res; - - Simplex_handle s_address; - // Root_simplex_const_iterator i; - for (auto i = s.begin(); i != s.end(); ++i) { - boost::optional<Vertex_handle> address = get_address(*i); - if (!address) - return res; - else - s_address.add_vertex(*address); - } - res = s_address; - return res; - } - - /** - * @brief returns a simplex with vertices which are the id of vertices of the - * argument. - */ - Root_simplex_handle get_id(const Simplex_handle& local_simplex) const { - Root_simplex_handle global_simplex; - for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) { - global_simplex.add_vertex(get_id(*x)); - } - return global_simplex; - } - - /** - * @brief returns true iff the simplex s belongs to the simplicial - * complex. - */ - virtual bool contains(const Simplex_handle & s) const { - if (s.dimension() == -1) { - return false; - } else if (s.dimension() == 0) { - return contains_vertex(s.first_vertex()); - } else { - return (contains_edges(s) && !blocks(s)); - } - } - - /* - * @brief returnrs true iff the complex is empty. - */ - bool empty() const { - return num_vertices() == 0; - } - - /* - * @brief returns the number of vertices in the complex. - */ - int num_vertices() const { - // remark boost::num_vertices(skeleton) counts deactivated vertices - return num_vertices_; - } - - /* - * @brief returns the number of edges in the complex. - * @details currently in O(n) - */ - // todo cache the value - - int num_edges() const { - return boost::num_edges(skeleton); - } - - int num_triangles() const{ - auto triangles = triangle_range(); - return std::distance(triangles.begin(),triangles.end()); - } - /* - * @brief returns the number of blockers in the complex. - */ - int num_blockers() const { - return num_blockers_; - } - - /* - * @brief returns true iff the graph of the 1-skeleton of the complex is complete. - */ - bool complete() const { - return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges(); - } - - /** - * @brief returns the number of connected components in the graph of the 1-skeleton. - */ - int num_connected_components() const { - int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); - std::vector<int> component(skeleton.vertex_set().size()); - return boost::connected_components(this->skeleton, &component[0]) - - num_vert_collapsed; - } - - /** - * @brief %Test if the complex is a cone. - * @details Runs in O(n) where n is the number of vertices. - */ - bool is_cone() const { - if (num_vertices() == 0) - return false; - if (num_vertices() == 1) - return true; - for (auto vi : vertex_range()) { - // xxx todo faire une methode bool is_in_blocker(Vertex_handle) - if (blocker_map_.find(vi) == blocker_map_.end()) { - // no blocker passes through the vertex, we just need to - // check if the current vertex is linked to all others vertices of the complex - if (degree_[vi.vertex] == num_vertices() - 1) - return true; - } - } - return false; - } - - //@} - /** @Simplification operations - */ - //@{ - - /** - * Returns true iff the blocker 'sigma' is popable. - * To define popable, let us call 'L' the complex that - * consists in the current complex without the blocker 'sigma'. - * A blocker 'sigma' is then "popable" if the link of 'sigma' - * in L is reducible. - * - */ - bool is_popable_blocker(Blocker_handle sigma) const; - - /** - * Removes all the popable blockers of the complex and delete them. - * @returns the number of popable blockers deleted - */ - void remove_popable_blockers(); - - /** - * Removes all the popable blockers of the complex passing through v and delete them. - */ - void remove_popable_blockers(Vertex_handle v); - - /** - * @brief Removes all the popable blockers of the complex passing through v and delete them. - * Also remove popable blockers in the neighborhood if they became popable. - * - */ - void remove_all_popable_blockers(Vertex_handle v); - - /** - * Remove the star of the vertex 'v' - */ - void remove_star(Vertex_handle v) ; - -private: - /** - * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. - * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added - * whenever the dimension of tau is at least 2. - */ - void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed); - -public: - /** - * Remove the star of the edge connecting vertices a and b. - * @returns the number of blocker that have been removed - */ - void remove_star(Vertex_handle a, Vertex_handle b) ; - - /** - * Remove the star of the edge 'e'. - */ - void remove_star(Edge_handle e) ; - - /** - * Remove the star of the simplex 'sigma' which needs to belong to the complex - */ - void remove_star(const Simplex_handle& sigma); - - /** - * @brief add a maximal simplex plus all its cofaces. - * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge). - */ - void add_simplex(const Simplex_handle& sigma); - -private: - /** - * remove all blockers that contains sigma - */ - void remove_blocker_containing_simplex(const Simplex_handle& sigma) ; - - /** - * remove all blockers that contains sigma - */ - void remove_blocker_include_in_simplex(const Simplex_handle& sigma); - -public: - enum simplifiable_status { - NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ - }; - - simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) { - // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex - // then call it there to build the link and return the value of link.is_contractible() - return MAYBE_HOMOTOPY_EQ; - } - - enum contractible_status { - NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE - }; - - /** - * @brief %Test if the complex is reducible using a strategy defined in the class - * (by default it tests if the complex is a cone) - * @details Note that NO could be returned if some invariant ensures that the complex - * is not a point (for instance if the euler characteristic is different from 1). - * This function will surely have to return MAYBE in some case because the - * associated problem is undecidable but it in practice, it can often - * be solved with the help of geometry. - */ - virtual contractible_status is_contractible() const { - if (this->is_cone()) { - return CONTRACTIBLE; - } else { - return MAYBE_CONTRACTIBLE; - } - } - //@} - - /** @Edge contraction operations - */ - //@{ - - /** - * @return If ignore_popable_blockers is true - * then the result is true iff the link condition at edge ab is satisfied - * or equivalently iff no blocker contains ab. - * If ignore_popable_blockers is false then the - * result is true iff all blocker containing ab are popable. - */ - bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{ - for (auto blocker : this->const_blocker_range(a)) - if (blocker->contains(b)) { - // false if ignore_popable_blockers is false - // otherwise the blocker has to be popable - return ignore_popable_blockers && is_popable_blocker(blocker); - } - return true; - - } - - /** - * @return If ignore_popable_blockers is true - * then the result is true iff the link condition at edge ab is satisfied - * or equivalently iff no blocker contains ab. - * If ignore_popable_blockers is false then the - * result is true iff all blocker containing ab are popable. - */ - bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const { - const Graph_edge& edge = (*this)[e]; - assert(this->get_address(edge.first())); - assert(this->get_address(edge.second())); - Vertex_handle a(*this->get_address(edge.first())); - Vertex_handle b(*this->get_address(edge.second())); - return link_condition(a, b, ignore_popable_blockers); - } - -protected: - /** - * Compute simplices beta such that a.beta is an order 0 blocker - * that may be used to construct a new blocker after contracting ab. - * It requires that the link condition is satisfied. - */ - void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const; - -private: - /** - * @brief "Replace" the edge 'bx' by the edge 'ax'. - * Assume that the edge 'bx' was present whereas 'ax' was not. - * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. - * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added - * and just before edge 'bx' had been removed. That way, it can - * eventually access to information of 'ax'. - */ - void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x); - -private: - /** - * @brief removes all blockers passing through the edge 'ab' - */ - void delete_blockers_around_vertex(Vertex_handle v); - - /** - * @brief removes all blockers passing through the edge 'ab' - */ - void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b); - -public: - /** - * Contracts the edge. - * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, - * it removes first all blockers passing through 'ab' - */ - void contract_edge(Edge_handle edge) { - contract_edge(this->first_vertex(edge), this->second_vertex(edge)); - } - - /** - * Contracts the edge connecting vertices a and b. - * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, - * it removes first all blockers passing through 'ab' - */ - void contract_edge(Vertex_handle a, Vertex_handle b); - -private: - void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add); - - /** - * delete all blockers that passes through a or b - */ - void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b); - void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ; - - void notify_changed_edges(Vertex_handle a) ; - //@} - -public: - ///////////////////////////////////////////////////////////////////////////// - /** @name Vertex iterators - */ - //@{ - typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator; - - // - // Range over the vertices of the simplicial complex. - // Methods .begin() and .end() return a Complex_vertex_iterator. - // - typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range; - - /** - * @brief Returns a Complex_vertex_range over all vertices of the complex - */ - Complex_vertex_range vertex_range() const { - auto begin = Complex_vertex_iterator(this); - auto end = Complex_vertex_iterator(this, 0); - return Complex_vertex_range(begin, end); - } - - typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator; - - - typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v - */ - Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const { - auto begin = Complex_neighbors_vertices_iterator(this, v); - auto end = Complex_neighbors_vertices_iterator(this, v, 0); - return Complex_neighbors_vertices_range(begin, end); - } - - //@} - - /** @name Edge iterators - */ - //@{ - - typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator; - - - typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex - */ - Complex_edge_range edge_range() const { - auto begin = Complex_edge_iterator(this); - auto end = Complex_edge_iterator(this, 0); - return Complex_edge_range(begin, end); - } - - - typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator; - - - typedef boost::iterator_range <Complex_edge_around_vertex_iterator> Complex_edge_around_vertex_range; - - /** - * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes - * through 'v' - */ - Complex_edge_around_vertex_range edge_range(Vertex_handle v) const { - auto begin = Complex_edge_around_vertex_iterator(this, v); - auto end = Complex_edge_around_vertex_iterator(this, v, 0); - return Complex_edge_around_vertex_range(begin, end); - } - - //@} - - /** @name Triangles iterators - */ - //@{ -private: - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; - typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; - -public: - typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator; - - typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range; - - /** - * @brief Range over triangles around a vertex of the simplicial complex. - * Methods .begin() and .end() return a Triangle_around_vertex_iterator. - * - */ - Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const { - auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v); - auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0); - return Complex_triangle_around_vertex_range(begin, end); - } - - typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; - - - typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator; - - - /** - * @brief Range over triangles of the simplicial complex. - * Methods .begin() and .end() return a Triangle_around_vertex_iterator. - * - */ - Complex_triangle_range triangle_range() const { - auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0); - if (empty()) { - return Complex_triangle_range(end, end); - } else { - auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); - return Complex_triangle_range(begin, end); - } - } - - //@} - - /** @name Simplices iterators - */ - //@{ - typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator; - - /** - * @brief Range over the simplices of the simplicial complex around a vertex. - * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. - */ - typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; - - /** - * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex - */ - Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const { - assert(contains_vertex(v)); - return Complex_simplex_around_vertex_range( - Complex_simplex_around_vertex_iterator(this, v), - Complex_simplex_around_vertex_iterator(this, v, true)); - } - - // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; - typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; - - typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; - - /** - * @brief Returns a Complex_simplex_range over all the simplices of the complex - */ - Complex_simplex_range simplex_range() const { - Complex_simplex_iterator end(this, true); - if (empty()) { - return Complex_simplex_range(end, end); - } else { - Complex_simplex_iterator begin(this); - return Complex_simplex_range(begin, end); - } - } - - //@} - - /** @name Blockers iterators - */ - //@{ -private: - /** - * @brief Iterator over the blockers adjacent to a vertex - */ - typedef Blocker_iterator_around_vertex_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, - Blocker_handle> - Complex_blocker_around_vertex_iterator; - - /** - * @brief Iterator over (constant) blockers adjacent to a vertex - */ - typedef Blocker_iterator_around_vertex_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, - const Blocker_handle> - Const_complex_blocker_around_vertex_iterator; - - typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; - typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range; - -public: - /** - * @brief Returns a range of the blockers of the complex passing through a vertex - */ - Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) { - auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); - auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); - return Complex_blocker_around_vertex_range(begin, end); - } - - /** - * @brief Returns a range of the blockers of the complex passing through a vertex - */ - Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const { - auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); - auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); - return Const_complex_blocker_around_vertex_range(begin, end); - } - -private: - /** - * @brief Iterator over the blockers. - */ - typedef Blocker_iterator_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, - Blocker_handle> - Complex_blocker_iterator; - - /** - * @brief Iterator over the (constant) blockers. - */ - typedef Blocker_iterator_internal< - typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, - const Blocker_handle> - Const_complex_blocker_iterator; - - typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; - typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; - -public: - /** - * @brief Returns a range of the blockers of the complex - */ - Complex_blocker_range blocker_range() { - auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); - auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); - return Complex_blocker_range(begin, end); - } - - /** - * @brief Returns a range of the blockers of the complex - */ - Const_complex_blocker_range const_blocker_range() const { - auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); - auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); - return Const_complex_blocker_range(begin, end); - } - - //@} - - ///////////////////////////////////////////////////////////////////////////// - /** @name Print and IO methods - */ - //@{ -public: - std::string to_string() const { - std::ostringstream stream; - stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl; - stream << num_edges() << " edges:\n" << edges_to_string() << std::endl; - stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl; - return stream.str(); - } - - std::string vertices_to_string() const { - std::ostringstream stream; - for (auto vertex : vertex_range()) { - stream << "{" << (*this)[vertex].get_id() << "} "; - } - stream << std::endl; - return stream.str(); - } - - std::string edges_to_string() const { - std::ostringstream stream; - for (auto edge : edge_range()) - stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} "; - stream << std::endl; - return stream.str(); - } - - std::string blockers_to_string() const { - std::ostringstream stream; - - for (auto b : const_blocker_range()) - stream << *b << std::endl; - return stream.str(); - } - //@} - - - - - + template<class ComplexType> friend class Vertex_iterator; + template<class ComplexType> friend class Neighbors_vertices_iterator; + template<class ComplexType> friend class Edge_iterator; + template<class ComplexType> friend class Edge_around_vertex_iterator; + + template<class ComplexType> friend class Skeleton_blocker_link_complex; + template<class ComplexType> friend class Skeleton_blocker_link_superior; + template<class ComplexType> friend class Skeleton_blocker_sub_complex; + + public: + /** + * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; + + /** + * @brief The type of stored edge node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; + + typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; + + /** + * @brief The type of an handle to a vertex of the complex. + */ + typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; + typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; + + /** + * @brief A ordered set of integers that represents a simplex. + */ + typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle; + typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; + + /** + * @brief Handle to a blocker of the complex. + */ + typedef Simplex_handle* Blocker_handle; + + typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; + typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator; + + protected: + typedef typename boost::adjacency_list<boost::setS, // edges + boost::vecS, // vertices + boost::undirectedS, Graph_vertex, Graph_edge> Graph; + // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink + // (eg Link with greater vertices) + // that burdens simplex iteration / complex initialization via list of simplices. + // to avoid that, one should modify the graph by storing two lists of adjacency for every + // vertex, the one with superior and the one with lower vertices, that way there is + // no more extra cost for computation of SuperiorLink + typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; + typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; + + protected: + typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; + + public: + /** + * @brief Handle to an edge of the complex. + */ + typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; + + protected: + typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator; + typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator; + + protected: + int num_vertices_; + int num_blockers_; + + typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; + // typedef Visitor* Visitor_ptr; + Visitor* visitor; + + /** + * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. + * + * This quantity is updated when adding/removing edge. + * + * This is useful because the operation + * list.size() is done in linear time. + */ + std::vector<boost_vertex_handle> degree_; + Graph skeleton; /** 1-skeleton of the simplicial complex. */ + + /** Each vertex can access to the blockers passing through it. */ + BlockerMap blocker_map_; + + public: + ///////////////////////////////////////////////////////////////////////////// + /** @name Constructors, Destructors + */ + //@{ + + /** + *@brief constructs a simplicial complex with a given number of vertices and a visitor. + */ + explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) + : visitor(visitor_) { + clear(); + for (int i = 0; i < num_vertices_; ++i) { + add_vertex(); + } + } + + private: + // typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie; + typedef Trie<Simplex_handle> STrie; + + public: + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + bool is_flag_complex = false, Visitor* visitor_ = NULL) + : num_vertices_(0), + num_blockers_(0), + visitor(visitor_) { + add_vertex_and_edges(simplex_begin, simplex_end); + + if (!is_flag_complex) + // need to compute blockers + add_blockers(simplex_begin, simplex_end); + } + + private: + template<typename SimpleHandleOutputIterator> + void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) { + std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; + // first pass to add vertices and edges + int num_vertex = -1; + for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { + if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex); + if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex()); + } + while (num_vertex-- >= 0) add_vertex(); + + for (const auto& e : edges) + add_edge(e.first, e.second); + } + + template<typename SimpleHandleOutputIterator> + void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) { + Tries<Simplex_handle> tries(num_vertices(), simplex_begin, simplex_end); + tries.init_next_dimension(); + auto simplices(tries.next_dimension_simplices()); + + while (!simplices.empty()) { + simplices = tries.next_dimension_simplices(); + for (auto& sigma : simplices) { + // common_positive_neighbors is the set of vertices u such that + // for all s in sigma, us is an edge and u>s + Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex())); + for (auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it) + if (sigma_it != sigma.rbegin()) + common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it)); + + for (auto x : common_positive_neighbors) { + // first test that all edges sx are here for all s in sigma + bool all_edges_here = true; + for (auto s : sigma) + if (!contains_edge(x, s)) { + all_edges_here = false; + break; + } + if (!all_edges_here) continue; + + // all edges of sigma \cup x are here + // we have a blocker if all proper faces of sigma \cup x + // are in the complex and if sigma \cup x is not in the complex + // the first is equivalent at checking if blocks(sigma \cup x) is true + // as all blockers of lower dimension have already been computed + sigma.add_vertex(x); + if (!tries.contains(sigma) && !blocks(sigma)) + add_blocker(sigma); + sigma.remove_vertex(x); + } + } + } + } + + public: + // We cannot use the default copy constructor since we need + // to make a copy of each of the blockers + + Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) { + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) { + add_blocker(*blocker); + } + } + + /** + */ + Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) { + clear(); + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) + add_blocker(*blocker); + return *this; + } + + /** + * return true if both complexes have the same simplices. + */ + bool operator==(const Skeleton_blocker_complex& other) const { + if (other.num_vertices() != num_vertices()) return false; + if (other.num_edges() != num_edges()) return false; + if (other.num_blockers() != num_blockers()) return false; + + for (auto v : vertex_range()) + if (!other.contains_vertex(v)) return false; + + for (auto e : edge_range()) + if (!other.contains_edge(first_vertex(e), second_vertex(e))) return false; + + for (const auto b : const_blocker_range()) + if (!other.contains_blocker(*b)) return false; + + return true; + } + + bool operator!=(const Skeleton_blocker_complex& other) const { + return !(*this == other); + } + + /** + * The destructor delete all blockers allocated. + */ + virtual ~Skeleton_blocker_complex() { + clear(); + } + + /** + * @details Clears the simplicial complex. After a call to this function, + * blockers are destroyed. The 1-skeleton and the set of blockers + * are both empty. + */ + virtual void clear() { + // xxx for now the responsabilty of freeing the visitor is for + // the user + visitor = NULL; + + degree_.clear(); + num_vertices_ = 0; + + remove_blockers(); + + skeleton.clear(); + } + + /** + *@brief allows to change the visitor. + */ + void set_visitor(Visitor* other_visitor) { + visitor = other_visitor; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertices operations + */ + //@{ + public: + /** + * @brief Return a local Vertex_handle of a vertex given a global one. + * @remark Assume that the vertex is present in the complex. + */ + Vertex_handle operator[](Root_vertex_handle global) const { + auto local(get_address(global)); + assert(local); + return *local; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + Graph_vertex& operator[](Vertex_handle address) { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + const Graph_vertex& operator[](Vertex_handle address) const { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. + */ + Vertex_handle add_vertex() { + Vertex_handle address(boost::add_vertex(skeleton)); + num_vertices_++; + (*this)[address].activate(); + // safe since we now that we are in the root complex and the field 'address' and 'id' + // are identical for every vertices + (*this)[address].set_id(Root_vertex_handle(address.vertex)); + degree_.push_back(0); + if (visitor) + visitor->on_add_vertex(address); + return address; + } + + /** + * @brief Remove a vertex from the simplicial complex + * @remark It just deactivates the vertex with a boolean flag but does not + * remove it from vertices from complexity issues. + */ + void remove_vertex(Vertex_handle address) { + assert(contains_vertex(address)); + // We remove b + boost::clear_vertex(address.vertex, skeleton); + (*this)[address].deactivate(); + num_vertices_--; + degree_[address.vertex] = -1; + if (visitor) + visitor->on_remove_vertex(address); + } + + /** + */ + bool contains_vertex(Vertex_handle u) const { + if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton)) + return false; + return (*this)[u].is_active(); + } + + /** + */ + bool contains_vertex(Root_vertex_handle u) const { + boost::optional<Vertex_handle> address = get_address(u); + return address && (*this)[*address].is_active(); + } + + /** + * @return true iff the simplicial complex contains all vertices + * of simplex sigma + */ + bool contains_vertices(const Simplex_handle & sigma) const { + for (auto vertex : sigma) + if (!contains_vertex(vertex)) + return false; + return true; + } + + /** + * @brief Given an Id return the address of the vertex having this Id in the complex. + * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex. + */ + virtual boost::optional<Vertex_handle> get_address( + Root_vertex_handle id) const { + boost::optional<Vertex_handle> res; + if (id.vertex < boost::num_vertices(skeleton)) + res = Vertex_handle(id.vertex); // xxx + return res; + } + + /** + * return the id of a vertex of adress local present in the graph + */ + Root_vertex_handle get_id(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return (*this)[local].get_id(); + } + + /** + * @brief Convert an address of a vertex of a complex to the address in + * the current complex. + * @details + * If the current complex is a sub (or sup) complex of 'other', it converts + * the address of a vertex v expressed in 'other' to the address of the vertex + * v in the current one. + * @remark this methods uses Root_vertex_handle to identify the vertex and + * assumes the vertex is present in the current complex. + */ + Vertex_handle convert_handle_from_another_complex(const Skeleton_blocker_complex& other, + Vertex_handle vh_in_other) const { + auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); + assert(vh_in_current_complex); + return *vh_in_current_complex; + } + + /** + * @brief return the graph degree of a vertex. + */ + int degree(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return degree_[local.vertex]; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Edges operations + */ + //@{ + public: + /** + * @brief return an edge handle if the two vertices forms + * an edge in the complex + */ + boost::optional<Edge_handle> operator[]( + const std::pair<Vertex_handle, Vertex_handle>& ab) const { + boost::optional<Edge_handle> res; + std::pair<Edge_handle, bool> edge_pair( + boost::edge(ab.first.vertex, ab.second.vertex, skeleton)); + if (edge_pair.second) + res = edge_pair.first; + return res; + } + + /** + * @brief returns the stored node associated to an edge + */ + Graph_edge& operator[](Edge_handle edge_handle) { + return skeleton[edge_handle]; + } + + /** + * @brief returns the stored node associated to an edge + */ + const Graph_edge& operator[](Edge_handle edge_handle) const { + return skeleton[edge_handle]; + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle first_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (source(edge_handle, skeleton)); + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle second_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (target(edge_handle, skeleton)); + } + + /** + * @brief returns the simplex made with the two vertices of the edge + * @details it assumes that the edge is present in the complex + + */ + Simplex_handle get_vertices(Edge_handle edge_handle) const { + auto edge((*this)[edge_handle]); + return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]); + } + + /** + * @brief Adds an edge between vertices a and b and all its cofaces. + */ + Edge_handle add_edge(Vertex_handle a, Vertex_handle b) { + assert(contains_vertex(a) && contains_vertex(b)); + assert(a != b); + + auto edge_handle((*this)[std::make_pair(a, b)]); + // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)]; + // Edge_handle edge_descr; + // bool edge_present = pair_descr_bool.second; + if (!edge_handle) { + edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first; + (*this)[*edge_handle].setId(get_id(a), get_id(b)); + degree_[a.vertex]++; + degree_[b.vertex]++; + if (visitor) + visitor->on_add_edge(a, b); + } + return *edge_handle; + } + + /** + * @brief Adds all edges and their cofaces of a simplex to the simplicial complex. + */ + void add_edges(const Simplex_handle & sigma) { + Simplex_handle_iterator i, j; + for (i = sigma.begin(); i != sigma.end(); ++i) + for (j = i, j++; j != sigma.end(); ++j) + add_edge(*i, *j); + } + + /** + * @brief Removes an edge from the simplicial complex and all its cofaces. + * @details returns the former Edge_handle representing the edge + */ + virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) { + bool found; + Edge_handle edge; + tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton); + if (found) { + if (visitor) + visitor->on_remove_edge(a, b); + boost::remove_edge(a.vertex, b.vertex, skeleton); + degree_[a.vertex]--; + degree_[b.vertex]--; + } + return edge; + } + + /** + * @brief Removes edge and its cofaces from the simplicial complex. + */ + void remove_edge(Edge_handle edge) { + assert(contains_vertex(first_vertex(edge))); + assert(contains_vertex(second_vertex(edge))); + remove_edge(first_vertex(edge), second_vertex(edge)); + } + + /** + * @brief The complex is reduced to its set of vertices. + * All the edges and blockers are removed. + */ + void keep_only_vertices() { + remove_blockers(); + + for (auto u : vertex_range()) { + while (this->degree(u) > 0) { + Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); + this->remove_edge(u, v); + } + } + } + + /** + * @return true iff the simplicial complex contains an edge between + * vertices a and b + */ + bool contains_edge(Vertex_handle a, Vertex_handle b) const { + // if (a.vertex<0 || b.vertex <0) return false; + return boost::edge(a.vertex, b.vertex, skeleton).second; + } + + /** + * @return true iff the simplicial complex contains all vertices + * and all edges of simplex sigma + */ + bool contains_edges(const Simplex_handle & sigma) const { + for (auto i = sigma.begin(); i != sigma.end(); ++i) { + if (!contains_vertex(*i)) + return false; + for (auto j = i; ++j != sigma.end();) { + if (!contains_edge(*i, *j)) + return false; + } + } + return true; + } + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Blockers operations + */ + //@{ + + /** + * @brief Adds the simplex to the set of blockers and + * returns a Blocker_handle toward it if was not present before and 0 otherwise. + */ + Blocker_handle add_blocker(const Simplex_handle& blocker) { + assert(blocker.dimension() > 1); + if (contains_blocker(blocker)) { + // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return 0; + } else { + if (visitor) + visitor->on_add_blocker(blocker); + Blocker_handle blocker_pt = new Simplex_handle(blocker); + num_blockers_++; + auto vertex = blocker_pt->begin(); + while (vertex != blocker_pt->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker_pt)); + ++vertex; + } + return blocker_pt; + } + } + + protected: + /** + * @brief Adds the simplex to the set of blockers + */ + void add_blocker(Blocker_handle blocker) { + if (contains_blocker(*blocker)) { + // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return; + } else { + if (visitor) + visitor->on_add_blocker(*blocker); + num_blockers_++; + auto vertex = blocker->begin(); + while (vertex != blocker->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker)); + ++vertex; + } + } + } + + protected: + /** + * Removes sigma from the blocker map of vertex v + */ + void remove_blocker(const Blocker_handle sigma, Vertex_handle v) { + Complex_blocker_around_vertex_iterator blocker; + for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end(); + ++blocker) { + if (*blocker == sigma) + break; + } + if (*blocker != sigma) { + std::cerr + << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; + assert(false); + } else { + blocker_map_.erase(blocker.current_position()); + } + } + + public: + /** + * @brief Removes the simplex from the set of blockers. + * @remark sigma has to belongs to the set of blockers + */ + void remove_blocker(const Blocker_handle sigma) { + for (auto vertex : *sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + + /** + * @brief Remove all blockers, in other words, it expand the simplicial + * complex to the smallest flag complex that contains it. + */ + void remove_blockers() { + // Desallocate the blockers + while (!blocker_map_.empty()) { + delete_blocker(blocker_map_.begin()->second); + } + num_blockers_ = 0; + blocker_map_.clear(); + } + + protected: + /** + * Removes the simplex sigma from the set of blockers. + * sigma has to belongs to the set of blockers + * + * @remark contrarily to delete_blockers does not call the destructor + */ + void remove_blocker(const Simplex_handle& sigma) { + assert(contains_blocker(sigma)); + for (auto vertex : sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + + public: + /** + * Removes the simplex s from the set of blockers + * and desallocate s. + */ + void delete_blocker(Blocker_handle sigma) { + if (visitor) + visitor->on_delete_blocker(sigma); + remove_blocker(sigma); + delete sigma; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Blocker_handle s) const { + if (s->dimension() < 2) + return false; + + Vertex_handle a = s->first_vertex(); + + for (const auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Simplex_handle & s) const { + if (s.dimension() < 2) + return false; + + Vertex_handle a = s.first_vertex(); + + for (auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + + private: + /** + * @return true iff a blocker of the simplicial complex + * is a face of sigma. + */ + bool blocks(const Simplex_handle & sigma) const { + for (auto s : sigma) + for (auto blocker : const_blocker_range(s)) + if (sigma.contains(*blocker)) + return true; + return false; + } + + //@} + + protected: + /** + * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$. + * If keep_only_superior is true then only vertices that are greater than v are added. + */ + virtual void add_neighbours(Vertex_handle v, Simplex_handle & n, + bool keep_only_superior = false) const { + boost_adjacency_iterator ai, ai_end; + for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; + ++ai) { + if (keep_only_superior) { + if (*ai > v.vertex) { + n.add_vertex(Vertex_handle(*ai)); + } + } else { + n.add_vertex(Vertex_handle(*ai)); + } + } + } + + /** + * @details Add to simplex res all vertices which are + * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. + * todo revoir + * + */ + virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res, + bool keep_only_superior = false) const { + res.clear(); + auto alpha_vertex = alpha.begin(); + add_neighbours(*alpha_vertex, res, keep_only_superior); + for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end(); + ++alpha_vertex) + keep_neighbours(*alpha_vertex, res, keep_only_superior); + } + + /** + * @details Remove from simplex n all vertices which are + * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. + */ + virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res, + bool keep_only_superior = false) const { + Simplex_handle nv; + add_neighbours(v, nv, keep_only_superior); + res.intersection(nv); + } + + /** + * @details Remove from simplex all vertices which are + * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are added. + */ + virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res, + bool keep_only_superior = false) const { + Simplex_handle nv; + add_neighbours(v, nv, keep_only_superior); + res.difference(nv); + } + + public: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Vertex_handle v) const { + return Link_complex(*this, Simplex_handle(v)); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Edge_handle edge) const { + return Link_complex(*this, edge); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(const Simplex_handle& simplex) const { + return Link_complex(*this, simplex); + } + + /** + * @brief Compute the local vertices of 's' in the current complex + * If one of them is not present in the complex then the return value is uninitialized. + * + * + */ + // xxx rename get_address et place un using dans sub_complex + + boost::optional<Simplex_handle> get_simplex_address( + const Root_simplex_handle& s) const { + boost::optional<Simplex_handle> res; + + Simplex_handle s_address; + // Root_simplex_const_iterator i; + for (auto i = s.begin(); i != s.end(); ++i) { + boost::optional<Vertex_handle> address = get_address(*i); + if (!address) + return res; + else + s_address.add_vertex(*address); + } + res = s_address; + return res; + } + + /** + * @brief returns a simplex with vertices which are the id of vertices of the + * argument. + */ + Root_simplex_handle get_id(const Simplex_handle& local_simplex) const { + Root_simplex_handle global_simplex; + for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) { + global_simplex.add_vertex(get_id(*x)); + } + return global_simplex; + } + + /** + * @brief returns true iff the simplex s belongs to the simplicial + * complex. + */ + virtual bool contains(const Simplex_handle & s) const { + if (s.dimension() == -1) { + return false; + } else if (s.dimension() == 0) { + return contains_vertex(s.first_vertex()); + } else { + return (contains_edges(s) && !blocks(s)); + } + } + + /* + * @brief returnrs true iff the complex is empty. + */ + bool empty() const { + return num_vertices() == 0; + } + + /* + * @brief returns the number of vertices in the complex. + */ + int num_vertices() const { + // remark boost::num_vertices(skeleton) counts deactivated vertices + return num_vertices_; + } + + /* + * @brief returns the number of edges in the complex. + * @details currently in O(n) + */ + // todo cache the value + + int num_edges() const { + return boost::num_edges(skeleton); + } + + int num_triangles() const { + auto triangles = triangle_range(); + return std::distance(triangles.begin(), triangles.end()); + } + + /* + * @brief returns the number of blockers in the complex. + */ + int num_blockers() const { + return num_blockers_; + } + + /* + * @brief returns true iff the graph of the 1-skeleton of the complex is complete. + */ + bool complete() const { + return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges(); + } + + /** + * @brief returns the number of connected components in the graph of the 1-skeleton. + */ + int num_connected_components() const { + int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); + std::vector<int> component(skeleton.vertex_set().size()); + return boost::connected_components(this->skeleton, &component[0]) + - num_vert_collapsed; + } + + /** + * @brief %Test if the complex is a cone. + * @details Runs in O(n) where n is the number of vertices. + */ + bool is_cone() const { + if (num_vertices() == 0) + return false; + if (num_vertices() == 1) + return true; + for (auto vi : vertex_range()) { + // xxx todo faire une methode bool is_in_blocker(Vertex_handle) + if (blocker_map_.find(vi) == blocker_map_.end()) { + // no blocker passes through the vertex, we just need to + // check if the current vertex is linked to all others vertices of the complex + if (degree_[vi.vertex] == num_vertices() - 1) + return true; + } + } + return false; + } + + //@} + /** @Simplification operations + */ + //@{ + + /** + * Returns true iff the blocker 'sigma' is popable. + * To define popable, let us call 'L' the complex that + * consists in the current complex without the blocker 'sigma'. + * A blocker 'sigma' is then "popable" if the link of 'sigma' + * in L is reducible. + * + */ + bool is_popable_blocker(Blocker_handle sigma) const; + + /** + * Removes all the popable blockers of the complex and delete them. + * @returns the number of popable blockers deleted + */ + void remove_popable_blockers(); + + /** + * Removes all the popable blockers of the complex passing through v and delete them. + */ + void remove_popable_blockers(Vertex_handle v); + + /** + * @brief Removes all the popable blockers of the complex passing through v and delete them. + * Also remove popable blockers in the neighborhood if they became popable. + * + */ + void remove_all_popable_blockers(Vertex_handle v); + + /** + * Remove the star of the vertex 'v' + */ + void remove_star(Vertex_handle v); + + private: + /** + * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. + * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added + * whenever the dimension of tau is at least 2. + */ + void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed); + + public: + /** + * Remove the star of the edge connecting vertices a and b. + * @returns the number of blocker that have been removed + */ + void remove_star(Vertex_handle a, Vertex_handle b); + + /** + * Remove the star of the edge 'e'. + */ + void remove_star(Edge_handle e); + + /** + * Remove the star of the simplex 'sigma' which needs to belong to the complex + */ + void remove_star(const Simplex_handle& sigma); + + /** + * @brief add a maximal simplex plus all its cofaces. + * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge). + */ + void add_simplex(const Simplex_handle& sigma); + + private: + /** + * remove all blockers that contains sigma + */ + void remove_blocker_containing_simplex(const Simplex_handle& sigma); + + /** + * remove all blockers that contains sigma + */ + void remove_blocker_include_in_simplex(const Simplex_handle& sigma); + + public: + enum simplifiable_status { + NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ + }; + + simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) { + // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of + // Skeleton_blocker_geometric_complex + // then call it there to build the link and return the value of link.is_contractible() + return MAYBE_HOMOTOPY_EQ; + } + + enum contractible_status { + NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE + }; + + /** + * @brief %Test if the complex is reducible using a strategy defined in the class + * (by default it tests if the complex is a cone) + * @details Note that NO could be returned if some invariant ensures that the complex + * is not a point (for instance if the euler characteristic is different from 1). + * This function will surely have to return MAYBE in some case because the + * associated problem is undecidable but it in practice, it can often + * be solved with the help of geometry. + */ + virtual contractible_status is_contractible() const { + if (this->is_cone()) { + return CONTRACTIBLE; + } else { + return MAYBE_CONTRACTIBLE; + } + } + //@} + + /** @Edge contraction operations + */ + //@{ + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const { + for (auto blocker : this->const_blocker_range(a)) + if (blocker->contains(b)) { + // false if ignore_popable_blockers is false + // otherwise the blocker has to be popable + return ignore_popable_blockers && is_popable_blocker(blocker); + } + return true; + } + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const { + const Graph_edge& edge = (*this)[e]; + assert(this->get_address(edge.first())); + assert(this->get_address(edge.second())); + Vertex_handle a(*this->get_address(edge.first())); + Vertex_handle b(*this->get_address(edge.second())); + return link_condition(a, b, ignore_popable_blockers); + } + + protected: + /** + * Compute simplices beta such that a.beta is an order 0 blocker + * that may be used to construct a new blocker after contracting ab. + * It requires that the link condition is satisfied. + */ + void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const; + + private: + /** + * @brief "Replace" the edge 'bx' by the edge 'ax'. + * Assume that the edge 'bx' was present whereas 'ax' was not. + * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. + * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added + * and just before edge 'bx' had been removed. That way, it can + * eventually access to information of 'ax'. + */ + void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x); + + private: + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_vertex(Vertex_handle v); + + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b); + + public: + /** + * Contracts the edge. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Edge_handle edge) { + contract_edge(this->first_vertex(edge), this->second_vertex(edge)); + } + + /** + * Contracts the edge connecting vertices a and b. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Vertex_handle a, Vertex_handle b); + + private: + void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, + std::set<Simplex_handle>& blockers_to_add); + /** + * delete all blockers that passes through a or b + */ + void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b); + void update_edges_after_contraction(Vertex_handle a, Vertex_handle b); + void notify_changed_edges(Vertex_handle a); + //@} + + public: + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertex iterators + */ + //@{ + typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator; + + // + // Range over the vertices of the simplicial complex. + // Methods .begin() and .end() return a Complex_vertex_iterator. + // + typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range; + + /** + * @brief Returns a Complex_vertex_range over all vertices of the complex + */ + Complex_vertex_range vertex_range() const { + auto begin = Complex_vertex_iterator(this); + auto end = Complex_vertex_iterator(this, 0); + return Complex_vertex_range(begin, end); + } + + typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator; + + + typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v + */ + Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const { + auto begin = Complex_neighbors_vertices_iterator(this, v); + auto end = Complex_neighbors_vertices_iterator(this, v, 0); + return Complex_neighbors_vertices_range(begin, end); + } + + //@} + + /** @name Edge iterators + */ + //@{ + + typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator; + + + typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex + */ + Complex_edge_range edge_range() const { + auto begin = Complex_edge_iterator(this); + auto end = Complex_edge_iterator(this, 0); + return Complex_edge_range(begin, end); + } + + + typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator; + + + typedef boost::iterator_range <Complex_edge_around_vertex_iterator> Complex_edge_around_vertex_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes + * through 'v' + */ + Complex_edge_around_vertex_range edge_range(Vertex_handle v) const { + auto begin = Complex_edge_around_vertex_iterator(this, v); + auto end = Complex_edge_around_vertex_iterator(this, v, 0); + return Complex_edge_around_vertex_range(begin, end); + } + + //@} + + /** @name Triangles iterators + */ + //@{ + private: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; + typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; + + public: + typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> + Superior_triangle_around_vertex_iterator; + typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > + Complex_triangle_around_vertex_range; + + /** + * @brief Range over triangles around a vertex of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const { + auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v); + auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0); + return Complex_triangle_around_vertex_range(begin, end); + } + + typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; + typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator; + + /** + * @brief Range over triangles of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_range triangle_range() const { + auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0); + if (empty()) { + return Complex_triangle_range(end, end); + } else { + auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); + return Complex_triangle_range(begin, end); + } + } + + //@} + + /** @name Simplices iterators + */ + //@{ + typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator; + + /** + * @brief Range over the simplices of the simplicial complex around a vertex. + * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. + */ + typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; + + /** + * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex + */ + Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const { + assert(contains_vertex(v)); + return Complex_simplex_around_vertex_range( + Complex_simplex_around_vertex_iterator(this, v), + Complex_simplex_around_vertex_iterator(this, v, true)); + } + + // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; + typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; + + typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; + + /** + * @brief Returns a Complex_simplex_range over all the simplices of the complex + */ + Complex_simplex_range simplex_range() const { + Complex_simplex_iterator end(this, true); + if (empty()) { + return Complex_simplex_range(end, end); + } else { + Complex_simplex_iterator begin(this); + return Complex_simplex_range(begin, end); + } + } + + //@} + + /** @name Blockers iterators + */ + //@{ + private: + /** + * @brief Iterator over the blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_around_vertex_iterator; + + /** + * @brief Iterator over (constant) blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_around_vertex_iterator; + + typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; + typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range; + + public: + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) { + auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Complex_blocker_around_vertex_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const { + auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Const_complex_blocker_around_vertex_range(begin, end); + } + + private: + /** + * @brief Iterator over the blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::iterator, + Blocker_handle> + Complex_blocker_iterator; + + /** + * @brief Iterator over the (constant) blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_iterator; + + typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; + typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; + + public: + /** + * @brief Returns a range of the blockers of the complex + */ + Complex_blocker_range blocker_range() { + auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Complex_blocker_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex + */ + Const_complex_blocker_range const_blocker_range() const { + auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Const_complex_blocker_range(begin, end); + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Print and IO methods + */ + //@{ + public: + std::string to_string() const { + std::ostringstream stream; + stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl; + stream << num_edges() << " edges:\n" << edges_to_string() << std::endl; + stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl; + return stream.str(); + } + + std::string vertices_to_string() const { + std::ostringstream stream; + for (auto vertex : vertex_range()) { + stream << "{" << (*this)[vertex].get_id() << "} "; + } + stream << std::endl; + return stream.str(); + } + + std::string edges_to_string() const { + std::ostringstream stream; + for (auto edge : edge_range()) + stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} "; + stream << std::endl; + return stream.str(); + } + + std::string blockers_to_string() const { + std::ostringstream stream; + + for (auto b : const_blocker_range()) + stream << *b << std::endl; + return stream.str(); + } + //@} }; - - /** * build a simplicial complex from a collection * of top faces. * return the total number of simplices */ template<typename Complex, typename SimplexHandleIterator> -Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, bool is_flag_complex = false) { - typedef typename Complex::Simplex_handle Simplex_handle; - std::vector<Simplex_handle> simplices; - for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) { - auto subfaces_topface = subfaces(*top_face); - simplices.insert(simplices.end(),subfaces_topface.begin(),subfaces_topface.end()); - } - return Complex(simplices.begin(),simplices.end(),is_flag_complex); +Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, + bool is_flag_complex = false) { + typedef typename Complex::Simplex_handle Simplex_handle; + std::vector<Simplex_handle> simplices; + for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) { + auto subfaces_topface = subfaces(*top_face); + simplices.insert(simplices.end(), subfaces_topface.begin(), subfaces_topface.end()); + } + return Complex(simplices.begin(), simplices.end(), is_flag_complex); } } // namespace skbl diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h index 69de1085..437482cb 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h @@ -37,195 +37,186 @@ namespace skbl { */ template<typename SkeletonBlockerGeometricDS> class Skeleton_blocker_geometric_complex : - public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> { - public: - typedef typename SkeletonBlockerGeometricDS::GT GT; - - typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; - - typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; - typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; - typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; - typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle; - - typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; - - typedef typename SkeletonBlockerGeometricDS::Point Point; - - - - Skeleton_blocker_geometric_complex(){ - } - - - - /** - * constructor given a list of points - */ - template<typename PointIterator> - explicit Skeleton_blocker_geometric_complex(int num_vertices,PointIterator begin,PointIterator end){ - for(auto point = begin; point != end; ++point) - add_vertex(*point); - } - - - /** - * @brief Constructor with a list of simplices. - * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). - */ - template<typename SimpleHandleOutputIterator,typename PointIterator> - Skeleton_blocker_geometric_complex( - SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, - PointIterator points_begin,PointIterator points_end, - bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){ - unsigned current = 0; - for(auto point = points_begin; point != points_end; ++point) - (*this)[Vertex_handle(current++)].point() = Point(point->begin(),point->end()); - } - - template<typename SimpleHandleOutputIterator,typename PointIterator> - friend Skeleton_blocker_geometric_complex make_complex_from_top_faces( - SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, - PointIterator points_begin,PointIterator points_end, - bool is_flag_complex = false){ - Skeleton_blocker_geometric_complex complex; - unsigned current = 0; - complex=make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin,simplex_end,is_flag_complex); - for(auto point = points_begin; point != points_end; ++point) - // complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end()); - complex.point(Vertex_handle(current++)) = Point(*point); - return complex; - } - - - /** - * @brief Constructor with a list of simplices. - * Points of every vertex are the point constructed with default constructor. - * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). - */ - template<typename SimpleHandleOutputIterator> - Skeleton_blocker_geometric_complex( - SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, - bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){ - } - - - /** - * @brief Add a vertex to the complex with a default constructed associated point. - */ - Vertex_handle add_vertex() { - return SimplifiableSkeletonblocker::add_vertex(); - } - - /** - * @brief Add a vertex to the complex with its associated point. - */ - Vertex_handle add_vertex(const Point& point) { - Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); - (*this)[ad].point() = point; - return ad; - } - - /** - * @brief Returns the Point associated to the vertex v. - */ - const Point& point(Vertex_handle v) const { - assert(this->contains_vertex(v)); - return (*this)[v].point(); - } - - /** - * @brief Returns the Point associated to the vertex v. - */ - Point& point(Vertex_handle v) { - assert(this->contains_vertex(v)); - return (*this)[v].point(); - } - - const Point& point(Root_vertex_handle global_v) const { - Vertex_handle local_v((*this)[global_v]); - assert(this->contains_vertex(local_v)); - return (*this)[local_v].point(); - } - - Point& point(Root_vertex_handle global_v) { - Vertex_handle local_v((*this)[global_v]); - assert(this->contains_vertex(local_v)); - return (*this)[local_v].point(); - } - - typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(Vertex_handle v) const { - Geometric_link link(*this, Simplex_handle(v)); - // we now add the point info - add_points_to_link(link); - return link; - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(const Simplex_handle& simplex) const { - Geometric_link link(*this, simplex); - // we now add the point info - add_points_to_link(link); - return link; - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Geometric_link link(Edge_handle edge) const { - Geometric_link link(*this, edge); - // we now add the point info - add_points_to_link(link); - return link; - } - - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link; - - /** - * Constructs the abstract link of v (without points coordinates). - */ - Abstract_link abstract_link(Vertex_handle v) const { - return Abstract_link(*this, Simplex_handle(v)); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Abstract_link abstract_link(const Simplex_handle& simplex) const { - return Abstract_link(*this, simplex); - } - - /** - * Constructs the link of 'simplex' with points coordinates. - */ - Abstract_link abstract_link(Edge_handle edge) const { - return Abstract_link(*this, edge); - } - - - - private: - - void add_points_to_link(Geometric_link& link) const { - for (Vertex_handle v : link.vertex_range()) { - Root_vertex_handle v_root(link.get_id(v)); - link.point(v) = (*this).point(v_root); - } - } - - - +public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> { + public: + typedef typename SkeletonBlockerGeometricDS::GT GT; + + typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; + + typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; + typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; + typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; + typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle; + + typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; + + typedef typename SkeletonBlockerGeometricDS::Point Point; + + Skeleton_blocker_geometric_complex() { } + + /** + * constructor given a list of points + */ + template<typename PointIterator> + explicit Skeleton_blocker_geometric_complex(int num_vertices, PointIterator begin, PointIterator end) { + for (auto point = begin; point != end; ++point) + add_vertex(*point); + } + + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be + * computed or not). + */ + template<typename SimpleHandleOutputIterator, typename PointIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + PointIterator points_begin, PointIterator points_end, + bool is_flag_complex = false) + : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) { + unsigned current = 0; + for (auto point = points_begin; point != points_end; ++point) + (*this)[Vertex_handle(current++)].point() = Point(point->begin(), point->end()); + } + + template<typename SimpleHandleOutputIterator, typename PointIterator> + friend Skeleton_blocker_geometric_complex make_complex_from_top_faces( + SimpleHandleOutputIterator simplex_begin, + SimpleHandleOutputIterator simplex_end, + PointIterator points_begin, + PointIterator points_end, + bool is_flag_complex = false) { + Skeleton_blocker_geometric_complex complex; + unsigned current = 0; + complex = + make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin, simplex_end, is_flag_complex); + for (auto point = points_begin; point != points_end; ++point) + // complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end()); + complex.point(Vertex_handle(current++)) = Point(*point); + return complex; + } + + /** + * @brief Constructor with a list of simplices. + * Points of every vertex are the point constructed with default constructor. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + bool is_flag_complex = false) + : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) { } + + /** + * @brief Add a vertex to the complex with a default constructed associated point. + */ + Vertex_handle add_vertex() { + return SimplifiableSkeletonblocker::add_vertex(); + } + + /** + * @brief Add a vertex to the complex with its associated point. + */ + Vertex_handle add_vertex(const Point& point) { + Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); + (*this)[ad].point() = point; + return ad; + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + const Point& point(Vertex_handle v) const { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + Point& point(Vertex_handle v) { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + const Point& point(Root_vertex_handle global_v) const { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + Point& point(Root_vertex_handle global_v) { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Vertex_handle v) const { + Geometric_link link(*this, Simplex_handle(v)); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(const Simplex_handle& simplex) const { + Geometric_link link(*this, simplex); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Edge_handle edge) const { + Geometric_link link(*this, edge); + // we now add the point info + add_points_to_link(link); + return link; + } + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link; + + /** + * Constructs the abstract link of v (without points coordinates). + */ + Abstract_link abstract_link(Vertex_handle v) const { + return Abstract_link(*this, Simplex_handle(v)); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Abstract_link abstract_link(const Simplex_handle& simplex) const { + return Abstract_link(*this, simplex); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Abstract_link abstract_link(Edge_handle edge) const { + return Abstract_link(*this, edge); + } + + private: + void add_points_to_link(Geometric_link& link) const { + for (Vertex_handle v : link.vertex_range()) { + Root_vertex_handle v_root(link.get_id(v)); + link.point(v) = (*this).point(v_root); + } + } }; -} // namespace skbl +} // namespace skbl -} // namespace Gudhi +} // namespace Gudhi #endif // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h index 9eab7f1e..86a12d90 100644 --- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -32,316 +32,52 @@ namespace Gudhi { namespace skbl { -///** -// * \brief Class that allows simplification operation on a simplicial complex represented -// * by a skeleton/blockers pair. -// * \ingroup skbl -// * @extends Skeleton_blocker_complex -// */ -//template<typename SkeletonBlockerDS> -//class Skeleton_blocker_complex : public Skeleton_blocker_complex<SkeletonBlockerDS> { -// template<class ComplexType> friend class Skeleton_blocker_sub_complex; -// -//public: -// typedef Skeleton_blocker_complex<SkeletonBlockerDS> SkeletonBlockerComplex; -// -// typedef typename SkeletonBlockerComplex::Graph_edge Graph_edge; -// -// typedef typename SkeletonBlockerComplex::boost_adjacency_iterator boost_adjacency_iterator; -// typedef typename SkeletonBlockerComplex::Edge_handle Edge_handle; -// typedef typename SkeletonBlockerComplex::boost_vertex_handle boost_vertex_handle; -// typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; -// typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle; -// typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle; -// typedef typename SkeletonBlockerComplex::Root_simplex_handle Root_simplex_handle; -// typedef typename SkeletonBlockerComplex::Blocker_handle Blocker_handle; -// -// -// typedef typename SkeletonBlockerComplex::Root_simplex_iterator Root_simplex_iterator; -// typedef typename SkeletonBlockerComplex::Simplex_handle_iterator Simplex_handle_iterator; -// typedef typename SkeletonBlockerComplex::BlockerMap BlockerMap; -// typedef typename SkeletonBlockerComplex::BlockerPair BlockerPair; -// typedef typename SkeletonBlockerComplex::BlockerMapIterator BlockerMapIterator; -// typedef typename SkeletonBlockerComplex::BlockerMapConstIterator BlockerMapConstIterator; -// -// typedef typename SkeletonBlockerComplex::Visitor Visitor; -// -// -// /** @name Constructors / Destructors / Initialization -// */ -// //@{ -// -// explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) : -// Skeleton_blocker_complex<SkeletonBlockerDS>(num_vertices_, visitor_) { } -// -// /** -// * @brief Constructor with a list of simplices -// * @details The list of simplices must be the list -// * of simplices of a simplicial complex, sorted with increasing dimension. -// */ -// //soon deprecated -// explicit Skeleton_blocker_complex(std::list<Simplex_handle>& simplices, Visitor* visitor_ = NULL) : -// Skeleton_blocker_complex<SkeletonBlockerDS>(simplices, visitor_) { } -// -// /** -// * @brief Constructor with a list of simplices -// * @details The list of simplices must be the list of simplices of a simplicial complex. -// */ -// template<typename SimpleHandleOutputIterator> -// explicit Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin,SimpleHandleOutputIterator simplex_end,bool is_flag_complex = false,Visitor* visitor_ = NULL) : -// Skeleton_blocker_complex<SkeletonBlockerDS>(simplex_begin,simplex_end,is_flag_complex, visitor_) { } -// -// -// virtual ~Skeleton_blocker_complex() { } -// -// //@} -// -// /** -// * Returns true iff the blocker 'sigma' is popable. -// * To define popable, let us call 'L' the complex that -// * consists in the current complex without the blocker 'sigma'. -// * A blocker 'sigma' is then "popable" if the link of 'sigma' -// * in L is reducible. -// * -// */ -// bool is_popable_blocker(Blocker_handle sigma) const; -// -// /** -// * Removes all the popable blockers of the complex and delete them. -// * @returns the number of popable blockers deleted -// */ -// void remove_popable_blockers(); -// -// /** -// * Removes all the popable blockers of the complex passing through v and delete them. -// */ -// void remove_popable_blockers(Vertex_handle v); -// -// /** -// * @brief Removes all the popable blockers of the complex passing through v and delete them. -// * Also remove popable blockers in the neighborhood if they became popable. -// * -// */ -// void remove_all_popable_blockers(Vertex_handle v); -// -// /** -// * Remove the star of the vertex 'v' -// */ -// void remove_star(Vertex_handle v) ; -// -//private: -// /** -// * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. -// * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added -// * whenever the dimension of tau is at least 2. -// */ -// void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed); -// -//public: -// /** -// * Remove the star of the edge connecting vertices a and b. -// * @returns the number of blocker that have been removed -// */ -// void remove_star(Vertex_handle a, Vertex_handle b) ; -// -// /** -// * Remove the star of the edge 'e'. -// */ -// void remove_star(Edge_handle e) ; -// -// /** -// * Remove the star of the simplex 'sigma' which needs to belong to the complex -// */ -// void remove_star(const Simplex_handle& sigma); -// -// /** -// * @brief add a maximal simplex plus all its cofaces. -// * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge). -// */ -// void add_simplex(const Simplex_handle& sigma); -// -//private: -// /** -// * remove all blockers that contains sigma -// */ -// void remove_blocker_containing_simplex(const Simplex_handle& sigma) ; -// -// /** -// * remove all blockers that contains sigma -// */ -// void remove_blocker_include_in_simplex(const Simplex_handle& sigma); -// -//public: -// enum simplifiable_status { -// NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ -// }; -// -// simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) { -// // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex -// // then call it there to build the link and return the value of link.is_contractible() -// return MAYBE_HOMOTOPY_EQ; -// } -// -// enum contractible_status { -// NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE -// }; -// -// /** -// * @brief %Test if the complex is reducible using a strategy defined in the class -// * (by default it tests if the complex is a cone) -// * @details Note that NO could be returned if some invariant ensures that the complex -// * is not a point (for instance if the euler characteristic is different from 1). -// * This function will surely have to return MAYBE in some case because the -// * associated problem is undecidable but it in practice, it can often -// * be solved with the help of geometry. -// */ -// virtual contractible_status is_contractible() const { -// if (this->is_cone()) { -// return CONTRACTIBLE; -// } else { -// return MAYBE_CONTRACTIBLE; -// } -// } -// -// /** @Edge contraction operations -// */ -// //@{ -// -// /** -// * @return If ignore_popable_blockers is true -// * then the result is true iff the link condition at edge ab is satisfied -// * or equivalently iff no blocker contains ab. -// * If ignore_popable_blockers is false then the -// * result is true iff all blocker containing ab are popable. -// */ -// bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{ -// for (auto blocker : this->const_blocker_range(a)) -// if (blocker->contains(b)) { -// // false if ignore_popable_blockers is false -// // otherwise the blocker has to be popable -// return ignore_popable_blockers && is_popable_blocker(blocker); -// } -// return true; -// -// } -// -// /** -// * @return If ignore_popable_blockers is true -// * then the result is true iff the link condition at edge ab is satisfied -// * or equivalently iff no blocker contains ab. -// * If ignore_popable_blockers is false then the -// * result is true iff all blocker containing ab are popable. -// */ -// bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const { -// const Graph_edge& edge = (*this)[e]; -// assert(this->get_address(edge.first())); -// assert(this->get_address(edge.second())); -// Vertex_handle a(*this->get_address(edge.first())); -// Vertex_handle b(*this->get_address(edge.second())); -// return link_condition(a, b, ignore_popable_blockers); -// } -// -//protected: -// /** -// * Compute simplices beta such that a.beta is an order 0 blocker -// * that may be used to construct a new blocker after contracting ab. -// * It requires that the link condition is satisfied. -// */ -// void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const; -// -//private: -// /** -// * @brief "Replace" the edge 'bx' by the edge 'ax'. -// * Assume that the edge 'bx' was present whereas 'ax' was not. -// * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. -// * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added -// * and just before edge 'bx' had been removed. That way, it can -// * eventually access to information of 'ax'. -// */ -// void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x); -// -//private: -// /** -// * @brief removes all blockers passing through the edge 'ab' -// */ -// void delete_blockers_around_vertex(Vertex_handle v); -// -// /** -// * @brief removes all blockers passing through the edge 'ab' -// */ -// void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b); -// -//public: -// /** -// * Contracts the edge. -// * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, -// * it removes first all blockers passing through 'ab' -// */ -// void contract_edge(Edge_handle edge) { -// contract_edge(this->first_vertex(edge), this->second_vertex(edge)); -// } -// -// /** -// * Contracts the edge connecting vertices a and b. -// * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, -// * it removes first all blockers passing through 'ab' -// */ -// void contract_edge(Vertex_handle a, Vertex_handle b); -// -//private: -// void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add); -// -// /** -// * delete all blockers that passes through a or b -// */ -// void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b); -// void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ; -// -// void notify_changed_edges(Vertex_handle a) ; -// //@} -//}; - - - +/** + * Returns true iff the blocker 'sigma' is popable. + * To define popable, let us call 'L' the complex that + * consists in the current complex without the blocker 'sigma'. + * A blocker 'sigma' is then "popable" if the link of 'sigma' + * in L is reducible. + * + */ template<typename SkeletonBlockerDS> bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_handle sigma) const { - assert(this->contains_blocker(*sigma)); - Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma; - build_link_of_blocker(*this, *sigma, link_blocker_sigma); + assert(this->contains_blocker(*sigma)); + Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma; + build_link_of_blocker(*this, *sigma, link_blocker_sigma); - bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE; - return res; + bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE; + return res; } - /** * Removes all the popable blockers of the complex and delete them. * @returns the number of popable blockers deleted */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() { - std::list<Vertex_handle> vertex_to_check; - for (auto v : this->vertex_range()) - vertex_to_check.push_front(v); - - while (!vertex_to_check.empty()) { - Vertex_handle v = vertex_to_check.front(); - vertex_to_check.pop_front(); - - bool blocker_popable_found = true; - while (blocker_popable_found) { - blocker_popable_found = false; - for (auto block : this->blocker_range(v)) { - if (this->is_popable_blocker(block)) { - for (Vertex_handle nv : *block) - if (nv != v) vertex_to_check.push_back(nv); - this->delete_blocker(block); - blocker_popable_found = true; - break; - } - } - } - } + std::list<Vertex_handle> vertex_to_check; + for (auto v : this->vertex_range()) + vertex_to_check.push_front(v); + + while (!vertex_to_check.empty()) { + Vertex_handle v = vertex_to_check.front(); + vertex_to_check.pop_front(); + + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (this->is_popable_blocker(block)) { + for (Vertex_handle nv : *block) + if (nv != v) vertex_to_check.push_back(nv); + this->delete_blocker(block); + blocker_popable_found = true; + break; + } + } + } + } } /** @@ -349,16 +85,16 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() { */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex_handle v) { - bool blocker_popable_found = true; - while (blocker_popable_found) { - blocker_popable_found = false; - for (auto block : this->blocker_range(v)) { - if (is_popable_blocker(block)) { - this->delete_blocker(block); - blocker_popable_found = true; - } - } - } + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (is_popable_blocker(block)) { + this->delete_blocker(block); + blocker_popable_found = true; + } + } + } } /** @@ -368,73 +104,83 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_all_popable_blockers(Vertex_handle v) { - std::list<Vertex_handle> vertex_to_check; - vertex_to_check.push_front(v); - - while (!vertex_to_check.empty()) { - Vertex_handle v = vertex_to_check.front(); - vertex_to_check.pop_front(); - - bool blocker_popable_found = true; - while (blocker_popable_found) { - blocker_popable_found = false; - for (auto block : this->blocker_range(v)) { - if (this->is_popable_blocker(block)) { - for (Vertex_handle nv : *block) - if (nv != v) vertex_to_check.push_back(nv); - this->delete_blocker(block); - blocker_popable_found = true; - break; - } - } - } - } + std::list<Vertex_handle> vertex_to_check; + vertex_to_check.push_front(v); + + while (!vertex_to_check.empty()) { + Vertex_handle v = vertex_to_check.front(); + vertex_to_check.pop_front(); + + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (this->is_popable_blocker(block)) { + for (Vertex_handle nv : *block) + if (nv != v) vertex_to_check.push_back(nv); + this->delete_blocker(block); + blocker_popable_found = true; + break; + } + } + } + } } - - +/** + * Remove the star of the vertice 'v' + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) { - // we remove the blockers that are not consistent anymore - update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v)); - while (this->degree(v) > 0) { - Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first)); - this->remove_edge(v, w); - } - this->remove_vertex(v); + // we remove the blockers that are not consistent anymore + update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v)); + while (this->degree(v) > 0) { + Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first)); + this->remove_edge(v, w); + } + this->remove_vertex(v); } +/** + * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. + * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added + * whenever the dimension of tau is at least 2. + */ template<typename SkeletonBlockerDS> -void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed) { - std::list <Blocker_handle> blockers_to_update; - if (simplex_to_be_removed.empty()) return; - - auto v0 = simplex_to_be_removed.first_vertex(); - for (auto blocker : this->blocker_range(v0)) { - if (blocker->contains(simplex_to_be_removed)) - blockers_to_update.push_back(blocker); - } - - for (auto blocker_to_update : blockers_to_update) { - Simplex_handle sub_blocker_to_be_added; - bool sub_blocker_need_to_be_added = - (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2; - if (sub_blocker_need_to_be_added) { - sub_blocker_to_be_added = *blocker_to_update; - sub_blocker_to_be_added.difference(simplex_to_be_removed); - } - this->delete_blocker(blocker_to_update); - if (sub_blocker_need_to_be_added) - this->add_blocker(sub_blocker_to_be_added); - } +void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge( + const Simplex_handle& simplex_to_be_removed) { + std::list <Blocker_handle> blockers_to_update; + if (simplex_to_be_removed.empty()) return; + + auto v0 = simplex_to_be_removed.first_vertex(); + for (auto blocker : this->blocker_range(v0)) { + if (blocker->contains(simplex_to_be_removed)) + blockers_to_update.push_back(blocker); + } + + for (auto blocker_to_update : blockers_to_update) { + Simplex_handle sub_blocker_to_be_added; + bool sub_blocker_need_to_be_added = + (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2; + if (sub_blocker_need_to_be_added) { + sub_blocker_to_be_added = *blocker_to_update; + sub_blocker_to_be_added.difference(simplex_to_be_removed); + } + this->delete_blocker(blocker_to_update); + if (sub_blocker_need_to_be_added) + this->add_blocker(sub_blocker_to_be_added); + } } - +/** + * Remove the star of the edge connecting vertices a and b. + * @returns the number of blocker that have been removed + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) { - update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b)); - // we remove the edge - this->remove_edge(a, b); + update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b)); + // we remove the edge + this->remove_edge(a, b); } /** @@ -442,22 +188,23 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, V */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Edge_handle e) { - return remove_star(this->first_vertex(e), this->second_vertex(e)); + return remove_star(this->first_vertex(e), this->second_vertex(e)); } - - +/** + * Remove the star of the simplex 'sigma' which needs to belong to the complex + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_handle& sigma) { - assert(this->contains(sigma)); - if (sigma.dimension() == 0) { - remove_star(sigma.first_vertex()); - } else if (sigma.dimension() == 1) { - remove_star(sigma.first_vertex(), sigma.last_vertex()); - } else { - remove_blocker_containing_simplex(sigma); - this->add_blocker(sigma); - } + assert(this->contains(sigma)); + if (sigma.dimension() == 0) { + remove_star(sigma.first_vertex()); + } else if (sigma.dimension() == 1) { + remove_star(sigma.first_vertex(), sigma.last_vertex()); + } else { + remove_blocker_containing_simplex(sigma); + this->add_blocker(sigma); + } } /** @@ -466,35 +213,34 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_hand */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_handle& sigma) { - assert(!this->contains(sigma)); - assert(sigma.dimension() > 1); - - int num_vertex_to_add = 0; - for(auto v : sigma) - if(!contains_vertex(v)) ++num_vertex_to_add; - while(num_vertex_to_add--) add_vertex(); - - for(auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it) - for(auto v_it = u_it; ++v_it != sigma.end(); /**/){ - std::cout <<"add edge"<<*u_it<<" "<<*v_it<<std::endl; - add_edge(*u_it,*v_it); - } - remove_blocker_include_in_simplex(sigma); + assert(!this->contains(sigma)); + assert(sigma.dimension() > 1); + + int num_vertex_to_add = 0; + for (auto v : sigma) + if (!contains_vertex(v)) ++num_vertex_to_add; + while (num_vertex_to_add--) add_vertex(); + + for (auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it) + for (auto v_it = u_it; ++v_it != sigma.end(); /**/) { + std::cout << "add edge" << *u_it << " " << *v_it << std::endl; + add_edge(*u_it, *v_it); + } + remove_blocker_include_in_simplex(sigma); } - /** * remove all blockers that contains sigma */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex_handle& sigma) { - std::vector <Blocker_handle> blockers_to_remove; - for (auto blocker : this->blocker_range(sigma.first_vertex())) { - if (blocker->contains(sigma)) - blockers_to_remove.push_back(blocker); - } - for (auto blocker_to_update : blockers_to_remove) - this->delete_blocker(blocker_to_update); + std::vector <Blocker_handle> blockers_to_remove; + for (auto blocker : this->blocker_range(sigma.first_vertex())) { + if (blocker->contains(sigma)) + blockers_to_remove.push_back(blocker); + } + for (auto blocker_to_update : blockers_to_remove) + this->delete_blocker(blocker_to_update); } /** @@ -502,173 +248,195 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simp */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex_handle& sigma) { - std::vector <Blocker_handle> blockers_to_remove; - for (auto blocker : this->blocker_range(sigma.first_vertex())) { - if (sigma.contains(*blocker)) - blockers_to_remove.push_back(blocker); - } - for (auto blocker_to_update : blockers_to_remove) - this->delete_blocker(blocker_to_update); + std::vector <Blocker_handle> blockers_to_remove; + for (auto blocker : this->blocker_range(sigma.first_vertex())) { + if (sigma.contains(*blocker)) + blockers_to_remove.push_back(blocker); + } + for (auto blocker_to_update : blockers_to_remove) + this->delete_blocker(blocker_to_update); } - +/** + * Compute simplices beta such that a.beta is an order 0 blocker + * that may be used to construct a new blocker after contracting ab. + * It requires that the link condition is satisfied. + */ template<typename SkeletonBlockerDS> -void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const { - for (auto const & blocker : this->const_blocker_range(a)) { - Simplex_handle beta = (*blocker); - beta.remove_vertex(a); - buffer.push_back(beta); - } - - Simplex_handle n; - this->add_neighbours(b, n); - this->remove_neighbours(a, n); - n.remove_vertex(a); - - - for (Vertex_handle y : n) { - Simplex_handle beta; - beta.add_vertex(y); - buffer.push_back(beta); - } +void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b, + std::vector<Simplex_handle> & buffer) const { + for (auto const & blocker : this->const_blocker_range(a)) { + Simplex_handle beta = (*blocker); + beta.remove_vertex(a); + buffer.push_back(beta); + } + + Simplex_handle n; + this->add_neighbours(b, n); + this->remove_neighbours(a, n); + n.remove_vertex(a); + + + for (Vertex_handle y : n) { + Simplex_handle beta; + beta.add_vertex(y); + buffer.push_back(beta); + } } +/** + * @brief "Replace" the edge 'bx' by the edge 'ax'. + * Assume that the edge 'bx' was present whereas 'ax' was not. + * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. + * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added + * and just before edge 'bx' had been removed. That way, it can + * eventually access to information of 'ax'. + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) { - this->add_edge(a, x); - if (this->visitor) this->visitor->on_swaped_edge(a, b, x); - this->remove_edge(b, x); + this->add_edge(a, x); + if (this->visitor) this->visitor->on_swaped_edge(a, b, x); + this->remove_edge(b, x); } template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertex(Vertex_handle v) { - std::list <Blocker_handle> blockers_to_delete; - for (auto blocker : this->blocker_range(v)) { - blockers_to_delete.push_back(blocker); - } - while (!blockers_to_delete.empty()) { - this->remove_blocker(blockers_to_delete.back()); - blockers_to_delete.pop_back(); - } + std::list <Blocker_handle> blockers_to_delete; + for (auto blocker : this->blocker_range(v)) { + blockers_to_delete.push_back(blocker); + } + while (!blockers_to_delete.empty()) { + this->remove_blocker(blockers_to_delete.back()); + blockers_to_delete.pop_back(); + } } +/** + * @brief removes all blockers passing through the edge 'ab' + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_edge(Vertex_handle a, Vertex_handle b) { - std::list<Blocker_handle> blocker_to_delete; - for (auto blocker : this->blocker_range(a)) - if (blocker->contains(b)) blocker_to_delete.push_back(blocker); - while (!blocker_to_delete.empty()) { - this->delete_blocker(blocker_to_delete.back()); - blocker_to_delete.pop_back(); - } + std::list<Blocker_handle> blocker_to_delete; + for (auto blocker : this->blocker_range(a)) + if (blocker->contains(b)) blocker_to_delete.push_back(blocker); + while (!blocker_to_delete.empty()) { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } } +/** + * Contracts the edge connecting vertices a and b. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vertex_handle b) { - assert(this->contains_vertex(a)); - assert(this->contains_vertex(b)); - assert(this->contains_edge(a, b)); + assert(this->contains_vertex(a)); + assert(this->contains_vertex(b)); + assert(this->contains_edge(a, b)); - // if some blockers passes through 'ab', we remove them. - if (!link_condition(a, b)) - delete_blockers_around_edge(a, b); + // if some blockers passes through 'ab', we remove them. + if (!link_condition(a, b)) + delete_blockers_around_edge(a, b); - std::set<Simplex_handle> blockers_to_add; + std::set<Simplex_handle> blockers_to_add; - get_blockers_to_be_added_after_contraction(a, b, blockers_to_add); + get_blockers_to_be_added_after_contraction(a, b, blockers_to_add); - delete_blockers_around_vertices(a, b); + delete_blockers_around_vertices(a, b); - update_edges_after_contraction(a, b); + update_edges_after_contraction(a, b); - this->remove_vertex(b); + this->remove_vertex(b); - notify_changed_edges(a); + notify_changed_edges(a); - for (auto block : blockers_to_add) - this->add_blocker(block); + for (auto block : blockers_to_add) + this->add_blocker(block); - assert(this->contains_vertex(a)); - assert(!this->contains_vertex(b)); + assert(this->contains_vertex(a)); + assert(!this->contains_vertex(b)); } template<typename SkeletonBlockerDS> void -Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add) { - blockers_to_add.clear(); - - typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType; - - LinkComplexType link_a(*this, a); - LinkComplexType link_b(*this, b); - - std::vector<Simplex_handle> vector_alpha, vector_beta; - - tip_blockers(a, b, vector_alpha); - tip_blockers(b, a, vector_beta); - - for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) { - for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) { - Simplex_handle sigma = *alpha; - sigma.union_vertices(*beta); - Root_simplex_handle sigma_id = this->get_id(sigma); - if (this->contains(sigma) && - proper_faces_in_union<Skeleton_blocker_complex<SkeletonBlockerDS>>(sigma_id, link_a, link_b)) { - // Blocker_handle blocker = new Simplex_handle(sigma); - sigma.add_vertex(a); - blockers_to_add.insert(sigma); - } - } - } +Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, + std::set<Simplex_handle>& blockers_to_add) { + blockers_to_add.clear(); + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType; + + LinkComplexType link_a(*this, a); + LinkComplexType link_b(*this, b); + + std::vector<Simplex_handle> vector_alpha, vector_beta; + + tip_blockers(a, b, vector_alpha); + tip_blockers(b, a, vector_beta); + + for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) { + for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) { + Simplex_handle sigma = *alpha; + sigma.union_vertices(*beta); + Root_simplex_handle sigma_id = this->get_id(sigma); + if (this->contains(sigma) && + proper_faces_in_union<Skeleton_blocker_complex < SkeletonBlockerDS >> (sigma_id, link_a, link_b)) { + // Blocker_handle blocker = new Simplex_handle(sigma); + sigma.add_vertex(a); + blockers_to_add.insert(sigma); + } + } + } } +/** + * delete all blockers that passes through a or b + */ template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b) { - std::vector<Blocker_handle> blocker_to_delete; - for (auto bl : this->blocker_range(a)) - blocker_to_delete.push_back(bl); - for (auto bl : this->blocker_range(b)) - blocker_to_delete.push_back(bl); - while (!blocker_to_delete.empty()) { - this->delete_blocker(blocker_to_delete.back()); - blocker_to_delete.pop_back(); - } + std::vector<Blocker_handle> blocker_to_delete; + for (auto bl : this->blocker_range(a)) + blocker_to_delete.push_back(bl); + for (auto bl : this->blocker_range(b)) + blocker_to_delete.push_back(bl); + while (!blocker_to_delete.empty()) { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } } - - template<typename SkeletonBlockerDS> void Skeleton_blocker_complex<SkeletonBlockerDS>::update_edges_after_contraction(Vertex_handle a, Vertex_handle b) { - // We update the set of edges - this->remove_edge(a, b); - - // For all edges {b,x} incident to b, - // we remove {b,x} and add {a,x} if not already there. - while (this->degree(b) > 0) { - Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first)); - if (!this->contains_edge(a, x)) - // we 'replace' the edge 'bx' by the edge 'ax' - this->swap_edge(a, b, x); - else - this->remove_edge(b, x); - } + // We update the set of edges + this->remove_edge(a, b); + + // For all edges {b,x} incident to b, + // we remove {b,x} and add {a,x} if not already there. + while (this->degree(b) > 0) { + Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first)); + if (!this->contains_edge(a, x)) + // we 'replace' the edge 'bx' by the edge 'ax' + this->swap_edge(a, b, x); + else + this->remove_edge(b, x); + } } template<typename SkeletonBlockerDS> void -Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a){ - // We notify the visitor that all edges incident to 'a' had changed - boost_adjacency_iterator v, v_end; - - for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v) - if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v)); +Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a) { + // We notify the visitor that all edges incident to 'a' had changed + boost_adjacency_iterator v, v_end; + for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v) + if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v)); } diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 1a96128f..c6968e44 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -273,16 +273,12 @@ private: if (!map.empty()) { std::cout << map.begin()->first; - if (map.begin()->second.children() == root()) - std::cout << "Sweet potato"; if (has_children(map.begin())) print_sc(map.begin()->second.children()); typename Dictionary::iterator it; for (it = map.begin()+1; it != map.end(); ++it) { std::cout << "," << it->first; - if (map.begin()->second.children() == root()) - std::cout << "Sweet potato"; if (has_children(it)) print_sc(it->second.children()); } diff --git a/src/Witness_complex/include/gudhi/Witness_complex1.h b/src/Witness_complex/include/gudhi/Witness_complex1.h index 6950c3bf..49ba7529 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex1.h +++ b/src/Witness_complex/include/gudhi/Witness_complex1.h @@ -334,10 +334,10 @@ private: if (j != i) { std::cout << "+++ We are at vertex=" << knn[witness_id][j] << std::endl; - if (curr_sibl->find(knn[witness_id][j]) == null_simplex()) + if (curr_sibl->members().find(knn[witness_id][j]) == null_simplex()) return false; std::cout << "++++ the simplex is there\n"; - curr_sh = curr_sibl->find(knn[witness_id][j]); + curr_sh = curr_sibl->members().find(knn[witness_id][j]); std::cout << "++++ curr_sh affectation is OK\n"; if (has_children(curr_sh)) curr_sibl = curr_sh->second.children(); |