diff options
Diffstat (limited to 'src')
20 files changed, 1159 insertions, 8 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 97d53306..a1b0b589 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -70,6 +70,7 @@ else() add_subdirectory(example/Persistent_cohomology) add_subdirectory(example/Skeleton_blocker) add_subdirectory(example/Contraction) + add_subdirectory(example/Witness_complex) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Doxyfile b/src/Doxyfile index 25f4f1ac..a629ea61 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -837,7 +837,7 @@ IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Contraction/ \ doc/Simplex_tree/ \ doc/Persistent_cohomology/ \ - + doc/Witness_complex/ # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/src/Persistent_cohomology/concept/FilteredComplex.h b/src/Persistent_cohomology/concept/FilteredComplex.h index e124d524..e47babe4 100644 --- a/src/Persistent_cohomology/concept/FilteredComplex.h +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -43,7 +43,7 @@ struct FilteredComplex * is model of IndexingTag. */ typedef unspecified Indexing_tag; -/** Returns a Simplex_hanlde that is different from all simplex handles +/** Returns a Simplex_handle that is different from all simplex handles * of the simplices. */ Simplex_handle null_simplex(); /** \brief Returns the number of simplices in the complex. @@ -61,12 +61,13 @@ struct FilteredComplex /** \brief Returns a key that is different from the keys associated * to the simplices. */ Simplex_key null_key (); -/** \brief Returns the key associated to a simplex. */ +/** \brief Returns the key associated to a simplex. + * + * This is never called on null_simplex(). */ Simplex_key key ( Simplex_handle sh ); -/** \brief Returns the simplex associated to a key. - * - * If key is different from null_key(), returns the simplex that - * has index idx in the filtration. */ +/** \brief Returns the simplex that has index idx in the filtration. + * + * This is never called on null_key(). */ Simplex_handle simplex ( Simplex_key idx ); /** \brief Assign a key to a simplex. */ void assign_key(Simplex_handle sh, Simplex_key key); diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index bea62de6..46670b85 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -578,7 +578,7 @@ class Simplex_tree { bool contiguous_vertices() const { if (root_.members_.empty()) return true; if (root_.members_.begin()->first != 0) return false; - if (std::prev(root_.members_.end())->first != root_.members_.size()-1) return false; + if (std::prev(root_.members_.end())->first != static_cast<Vertex_handle>(root_.members_.size() - 1)) return false; return true; } diff --git a/src/Witness_complex/concept/Simplicial_complex_for_witness.h b/src/Witness_complex/concept/Simplicial_complex_for_witness.h new file mode 100644 index 00000000..caaf0db6 --- /dev/null +++ b/src/Witness_complex/concept/Simplicial_complex_for_witness.h @@ -0,0 +1,87 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (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 CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ +#define CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ + +namespace Gudhi { + +namespace witness_complex { + +/** \brief The concept Simplicial_Complex describes the requirements + * for a type to implement a simplicial complex, + * used for example to build a 'Witness_complex'. + */ +struct SimplicialComplexForWitness { + /** Handle to specify a simplex. */ + typedef unspecified Simplex_handle; + /** Handle to specify a vertex. Must be a non-negative integer. */ + typedef unspecified Vertex_handle; + + /** Returns a Simplex_hanlde that is different from all simplex handles + * of the simplices. */ + Simplex_handle null_simplex(); + + /** \brief Iterator over the simplices of the complex, + * in an arbitrary order. + * + * 'value_type' must be 'Simplex_handle'.*/ + typedef unspecified Complex_simplex_range; + + /** + * \brief Returns a range over all the simplices of a + * complex. + */ + Complex_simplex_range complex_simplex_range(); + + /** \brief Iterator over vertices of a simplex. + * + * 'value type' must be 'Vertex_handle'.*/ + typedef unspecified Simplex_vertex_range; + + /** \brief Returns a range over vertices of a given + * simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle const & simplex); + + /** \brief Return type of an insertion of a simplex + */ + typedef unspecified Insertion_result_type; + + /** \brief Inserts a simplex with vertices from a given range + * 'vertex_range' in the simplicial complex. + * */ + template< typedef Input_vertex_range > + Insertion_result_type insert_simplex(Input_vertex_range const & vertex_range); + + /** \brief Finds a simplex with vertices given by a range + * + * If a simplex exists, its Simplex_handle is returned. + * Otherwise null_simplex() is returned. */ + template< typedef Input_vertex_range > + Simplex_handle find(Input_vertex_range const & vertex_range); +}; + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // CONCEPT_WITNESS_COMPLEX_SIMPLICIAL_COMPLEX_FOR_WITNESS_H_ diff --git a/src/Witness_complex/doc/Witness_complex_doc.h b/src/Witness_complex/doc/Witness_complex_doc.h new file mode 100644 index 00000000..60dfd27b --- /dev/null +++ b/src/Witness_complex/doc/Witness_complex_doc.h @@ -0,0 +1,42 @@ +#ifndef WITNESS_COMPLEX_DOC_H_ +#define WITNESS_COMPLEX_DOC_H_ + +/** + \defgroup witness_complex Witness complex + + \author Siargey Kachanovich + + \image html "Witness_complex_representation.png" "Witness complex representation" + + \section Definitions + + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$: + + \li \f$W\f$ set of **witnesses** and + \li \f$L \subseteq W\f$ set of **landmarks**. + + The simplices are based on landmarks + and a simplex belongs to the witness complex if and only if it is witnessed, that is: + + \f$ \sigma \subset L \f$ is witnessed if there exists a point \f$w \in W\f$ such that + w is closer to the vertices of \f$ \sigma \f$ than other points in \f$ L \f$ and all of its faces are witnessed as well. + + The data structure is described in \cite boissonnatmariasimplextreealgorithmica . + + \section Implementation + + The principal class of this module is Gudhi::Witness_complex. + + In both cases, the constructor for this class takes a {witness}x{closest_landmarks} table, where each row represents a witness and consists of landmarks sorted by distance to this witness. + This table can be constructed by two additional classes Landmark_choice_by_furthest_point and Landmark_choice_by_random_point also included in the module. + + *\image html "bench_Cy8.png" "Running time as function on number of landmarks" width=10cm + *\image html "bench_sphere.png" "Running time as function on number of witnesses for |L|=300" width=10cm + + + \copyright GNU General Public License v3. + + + */ + +#endif // WITNESS_COMPLEX_DOC_H_ diff --git a/src/Witness_complex/doc/Witness_complex_representation.png b/src/Witness_complex/doc/Witness_complex_representation.png Binary files differnew file mode 100644 index 00000000..1d31a490 --- /dev/null +++ b/src/Witness_complex/doc/Witness_complex_representation.png diff --git a/src/Witness_complex/doc/bench_Cy8.png b/src/Witness_complex/doc/bench_Cy8.png Binary files differnew file mode 100644 index 00000000..d9045294 --- /dev/null +++ b/src/Witness_complex/doc/bench_Cy8.png diff --git a/src/Witness_complex/doc/bench_sphere.png b/src/Witness_complex/doc/bench_sphere.png Binary files differnew file mode 100644 index 00000000..ba6bb381 --- /dev/null +++ b/src/Witness_complex/doc/bench_sphere.png diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt new file mode 100644 index 00000000..b304479e --- /dev/null +++ b/src/Witness_complex/example/CMakeLists.txt @@ -0,0 +1,44 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplex) + +# A simple example + add_executable( witness_complex_from_file witness_complex_from_file.cpp ) + add_test( witness_complex_from_bunny ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 100) + +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.6.0) + message(STATUS "CGAL version: ${CGAL_VERSION}.") + + include( ${CGAL_USE_FILE} ) + # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten. + # cf. http://doc.cgal.org/latest/Manual/installation.html#title40 + # A workaround is to add "-std=c++11" again. + # A fix would be to use https://cmake.org/cmake/help/v3.1/prop_gbl/CMAKE_CXX_KNOWN_FEATURES.html + # or even better https://cmake.org/cmake/help/v3.1/variable/CMAKE_CXX_STANDARD.html + # but it implies to use cmake version 3.1 at least. + if(NOT MSVC) + include(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11) + if(COMPILER_SUPPORTS_CXX11) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + endif() + endif() + # - End of workaround + + find_package(Eigen3 3.1.0) + if (EIGEN3_FOUND) + message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.") + include( ${EIGEN3_USE_FILE} ) + message(STATUS "Eigen3 use file: ${EIGEN3_USE_FILE}.") + include_directories (BEFORE "../../include") + + add_executable ( witness_complex_sphere witness_complex_sphere.cpp ) + target_link_libraries(witness_complex_sphere ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_test( witness_complex_sphere_10 ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_sphere 10) + else() + message(WARNING "Eigen3 not found. Version 3.1.0 is required for witness_complex_sphere example.") + endif() + else() + message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile witness_complex_sphere example. Version 4.6.0 is required.") + endif () +endif() diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h new file mode 100644 index 00000000..ac445261 --- /dev/null +++ b/src/Witness_complex/example/generators.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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ +#define EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> + +#include <fstream> +#include <string> +#include <vector> + +typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K; +typedef K::FT FT; +typedef K::Point_d Point_d; +typedef std::vector<Point_d> Point_Vector; +typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator; +typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator; + +/** + * \brief Rock age method of reading off file + * + */ +inline void +off_reader_cust(std::string file_name, std::vector<Point_d> & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + // Line OFF. No need in it + if (!getline(in_file, line)) { + std::cerr << "No line OFF\n"; + return; + } + // Line with 3 numbers. No need + if (!getline(in_file, line)) { + std::cerr << "No line with 3 numbers\n"; + return; + } + // Reading points + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + points.push_back(Point_d(point)); + } + in_file.close(); +} + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, Point_Vector & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +/** \brief Generate points on a grid in a cube of side 2 + * having {+-1}^D as vertices and insert them in W. + * The grid has "width" points on each side. + * If torus is true then it is supposed that the cube represents + * a flat torus, hence the opposite borders are associated. + * The points on border in this case are not placed twice. + */ +void generate_points_grid(Point_Vector& W, int width, int D, bool torus) { + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) { + std::vector<double> point; + int cell_i = i; + for (int l = 0; l < D; ++l) { + if (torus) + point.push_back(-1 + (2.0 / (width - 1))*(cell_i % width)); + else + point.push_back(-1 + (2.0 / width)*(cell_i % width)); + // attention: the bottom and the right are covered too! + cell_i /= width; + } + W.push_back(point); + } +} + +/** \brief Generate nbP points uniformly in a cube of side 2 + * having {+-1}^dim as its vertices and insert them in W. + */ +void generate_points_random_box(Point_Vector& W, int nbP, int dim) { + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) { + W.push_back(*rp++); + } +} + +/** \brief Generate nbP points uniformly on a (dim-1)-sphere + * and insert them in W. + */ +void generate_points_sphere(Point_Vector& W, int nbP, int dim) { + CGAL::Random_points_on_sphere_d<Point_d> rp(dim, 1); + for (int i = 0; i < nbP; i++) + W.push_back(*rp++); +} + +#endif // EXAMPLE_WITNESS_COMPLEX_GENERATORS_H_ diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp new file mode 100644 index 00000000..53207ad2 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_from_file.cpp @@ -0,0 +1,100 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/reader_utils.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef std::vector< std::vector <double> > Point_Vector; + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + if (point.size() != 1) + points.push_back(point); + } + in_file.close(); +} + +int main(int argc, char * const argv[]) { + if (argc != 3) { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file nbL \n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]); + clock_t start, end; + + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_Vector point_vector; + read_points_cust(file_name, point_vector); + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, nbL, knn); + end = clock(); + std::cout << "Landmark choice for " << nbL << " landmarks took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + + // Compute witness complex + start = clock(); + Gudhi::witness_complex::witness_complex(knn, nbL, point_vector[0].size(), simplex_tree); + end = clock(); + std::cout << "Witness complex took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; +} diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp new file mode 100644 index 00000000..b26c9f36 --- /dev/null +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -0,0 +1,90 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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_PARAMETER_MAX_ARITY 12 + + +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/reader_utils.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <utility> +#include <string> +#include <vector> + +#include "generators.h" + +/** Write a gnuplot readable file. + * Data range is a random access range of pairs (arg, value) + */ +template < typename Data_range > +void write_data(Data_range & data, std::string filename) { + std::ofstream ofs(filename, std::ofstream::out); + for (auto entry : data) + ofs << entry.first << ", " << entry.second << "\n"; + ofs.close(); +} + +int main(int argc, char * const argv[]) { + if (argc != 2) { + std::cerr << "Usage: " << argv[0] + << " number_of_landmarks \n"; + return 0; + } + + int number_of_landmarks = atoi(argv[1]); + clock_t start, end; + + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + + std::vector< std::pair<int, double> > l_time; + + // Read the point file + for (int nbP = 500; nbP < 10000; nbP += 500) { + Point_Vector point_vector; + generate_points_sphere(point_vector, nbP, 4); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + std::vector<std::vector< int > > knn; + Gudhi::witness_complex::landmark_choice_by_random_point(point_vector, number_of_landmarks, knn); + + // Compute witness complex + Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); + end = clock(); + double time = static_cast<double>(end - start) / CLOCKS_PER_SEC; + std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " + << time << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; + l_time.push_back(std::make_pair(nbP, time)); + } + write_data(l_time, "w_time.dat"); +} diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h new file mode 100644 index 00000000..df93155b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h @@ -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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ +#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ + +#include <boost/range/size.hpp> + +#include <limits> // for numeric_limits<> +#include <iterator> +#include <algorithm> // for sort +#include <vector> + +namespace Gudhi { + +namespace witness_complex { + + typedef std::vector<int> typeVectorVertex; + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by iteratively adding the furthest witness from the + * current landmark set as the new landmark. + * \details It chooses nbL landmarks from a random access range `points` and + * writes {witness}*{closest landmarks} matrix in `knn`. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges + * + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_furthest_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nb_points = boost::size(points); + assert(nb_points >= nbL); + // distance matrix witness x landmarks + std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>()); + // landmark list + typeVectorVertex chosen_landmarks; + + knn = KNearestNeighbours(nb_points, std::vector<int>()); + int current_number_of_landmarks = 0; // counter for landmarks + double curr_max_dist = 0; // used for defining the furhest point from L + const double infty = std::numeric_limits<double>::infinity(); // infinity (see next entry) + std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from points + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + // or better yet std::uniform_int_distribution + int rand_int = rand() % nb_points; + int curr_max_w = rand_int; // For testing purposes a pseudo-random number is used here + + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + // curr_max_w at this point is the next landmark + chosen_landmarks.push_back(curr_max_w); + unsigned i = 0; + for (auto& p : points) { + double curr_dist = euclidean_distance(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks])); + wit_land_dist[i].push_back(curr_dist); + knn[i].push_back(current_number_of_landmarks); + if (curr_dist < dist_to_L[i]) + dist_to_L[i] = curr_dist; + ++i; + } + curr_max_dist = 0; + for (i = 0; i < dist_to_L.size(); i++) + if (dist_to_L[i] > curr_max_dist) { + curr_max_dist = dist_to_L[i]; + curr_max_w = i; + } + } + for (int i = 0; i < nb_points; ++i) + std::sort(std::begin(knn[i]), + std::end(knn[i]), + [&wit_land_dist, i](int a, int b) { + return wit_land_dist[i][a] < wit_land_dist[i][b]; }); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h new file mode 100644 index 00000000..ebf6aad1 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 LANDMARK_CHOICE_BY_RANDOM_POINT_H_ +#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_ + +#include <boost/range/size.hpp> + +#include <queue> // for priority_queue<> +#include <utility> // for pair<> +#include <iterator> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace witness_complex { + + /** + * \ingroup witness_complex + * \brief Landmark choice strategy by taking random vertices for landmarks. + * \details It chooses nbL distinct landmarks from a random access range `points` + * and outputs a matrix {witness}*{closest landmarks} in knn. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges and + * Vertex_handle is the label type of a vertex in a simplicial complex. + * Closest_landmark_range needs to have push_back operation. + */ + + template <typename KNearestNeighbours, + typename Point_random_access_range> + void landmark_choice_by_random_point(Point_random_access_range const &points, + int nbL, + KNearestNeighbours &knn) { + int nbP = boost::size(points); + assert(nbP >= nbL); + std::set<int> landmarks; + int current_number_of_landmarks = 0; // counter for landmarks + + // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety + int chosen_landmark = rand() % nbP; + for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) { + while (landmarks.find(chosen_landmark) != landmarks.end()) + chosen_landmark = rand() % nbP; + landmarks.insert(chosen_landmark); + } + + int dim = boost::size(*std::begin(points)); + typedef std::pair<double, int> dist_i; + typedef bool (*comp)(dist_i, dist_i); + knn = KNearestNeighbours(nbP); + for (int points_i = 0; points_i < nbP; points_i++) { + std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) { + return j1.first > j2.first; + }); + std::set<int>::iterator landmarks_it; + int landmarks_i = 0; + for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end(); + ++landmarks_it, landmarks_i++) { + dist_i dist = std::make_pair(euclidean_distance(points[points_i], points[*landmarks_it]), landmarks_i); + l_heap.push(dist); + } + for (int i = 0; i < dim + 1; i++) { + dist_i dist = l_heap.top(); + knn[points_i].push_back(dist.second); + l_heap.pop(); + } + } + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // LANDMARK_CHOICE_BY_RANDOM_POINT_H_ diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h new file mode 100644 index 00000000..489cdf11 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -0,0 +1,265 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (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 WITNESS_COMPLEX_H_ +#define WITNESS_COMPLEX_H_ + +// Needed for the adjacency graph in bad link search +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +#include <boost/range/size.hpp> + +#include <gudhi/distance_functions.h> + +#include <algorithm> +#include <utility> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <ctime> +#include <iostream> + +namespace Gudhi { + +namespace witness_complex { + +// /* +// * \private +// \class Witness_complex +// \brief Constructs the witness complex for the given set of witnesses and landmarks. +// \ingroup witness_complex +// */ +template< class SimplicialComplex> +class Witness_complex { + private: + struct Active_witness { + int witness_id; + int landmark_id; + + Active_witness(int witness_id_, int landmark_id_) + : witness_id(witness_id_), + landmark_id(landmark_id_) { } + }; + + private: + typedef typename SimplicialComplex::Simplex_handle Simplex_handle; + typedef typename SimplicialComplex::Vertex_handle Vertex_handle; + + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; + + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; + + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; + + private: + int nbL_; // Number of landmarks + SimplicialComplex& sc_; // Simplicial complex + + public: + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /* @name Constructor + */ + + //@{ + + // Witness_range<Closest_landmark_range<Vertex_handle>> + + /* + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Constructor takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + template< typename KNearestNeighbors > + Witness_complex(KNearestNeighbors const & knn, + int nbL, + int dim, + SimplicialComplex & sc) : nbL_(nbL), sc_(sc) { + // Construction of the active witness list + int nbW = boost::size(knn); + typeVectorVertex vv; + int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations + */ + ActiveWitnessList active_w; // = new ActiveWitnessList(); + for (Vertex_handle i = 0; i != nbL_; ++i) { + // initial fill of 0-dimensional simplices + // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore + counter++; + vv = {i}; + sc_.insert_simplex(vv); + // TODO(SK) Error if not inserted : normally no need here though + } + int k = 1; /* current dimension in iterative construction */ + for (int i = 0; i != nbW; ++i) + active_w.push_back(i); + while (!active_w.empty() && k < dim) { + typename ActiveWitnessList::iterator it = active_w.begin(); + while (it != active_w.end()) { + typeVectorVertex simplex_vector; + /* THE INSERTION: Checking if all the subfaces are in the simplex tree*/ + bool ok = all_faces_in(knn, *it, k); + if (ok) { + for (int i = 0; i != k + 1; ++i) + simplex_vector.push_back(knn[*it][i]); + sc_.insert_simplex(simplex_vector); + // TODO(SK) Error if not inserted : normally no need here though + ++it; + } else { + active_w.erase(it++); // First increase the iterator and then erase the previous element + } + } + k++; + } + } + + //@} + + private: + /* \brief Check if the facets of the k-dimensional simplex witnessed + * by witness witness_id are already in the complex. + * inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id + */ + template <typename KNearestNeighbors> + bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) { + std::vector< Vertex_handle > facet; + // CHECK ALL THE FACETS + for (int i = 0; i != k + 1; ++i) { + facet = {}; + for (int j = 0; j != k + 1; ++j) { + if (j != i) { + facet.push_back(knn[witness_id][j]); + } + } // endfor + if (sc_.find(facet) == sc_.null_simplex()) + return false; + } // endfor + return true; + } + + template <typename T> + static void print_vector(const std::vector<T>& v) { + std::cout << "["; + if (!v.empty()) { + std::cout << *(v.begin()); + for (auto it = v.begin() + 1; it != v.end(); ++it) { + std::cout << ","; + std::cout << *it; + } + } + std::cout << "]"; + } + + public: + // /* + // * \brief Verification if every simplex in the complex is witnessed by witnesses in knn. + // * \param print_output =true will print the witnesses for each simplex + // * \remark Added for debugging purposes. + // */ + template< class KNearestNeighbors > + bool is_witness_complex(KNearestNeighbors const & knn, bool print_output) { + for (Simplex_handle sh : sc_.complex_simplex_range()) { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; // number of verticed in the simplex + for (Vertex_handle v : sc_.simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w : knn) { + bool has_vertices = true; + for (Vertex_handle v : simplex) + if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) { + has_vertices = false; + } + if (has_vertices) { + is_witnessed = true; + if (print_output) { + std::cout << "The simplex "; + print_vector(simplex); + std::cout << " is witnessed by the witness "; + print_vector(w); + std::cout << std::endl; + } + break; + } + } + if (!is_witnessed) { + if (print_output) { + std::cout << "The following simplex is not witnessed "; + print_vector(simplex); + std::cout << std::endl; + } + assert(is_witnessed); + return false; + } + } + return true; + } +}; + + /** + * \ingroup witness_complex + * \brief Iterative construction of the witness complex. + * \details The witness complex is written in simplicial complex sc_ + * basing on a matrix knn of + * nearest neighbours of the form {witnesses}x{landmarks}. + * + * The type KNearestNeighbors can be seen as + * Witness_range<Closest_landmark_range<Vertex_handle>>, where + * Witness_range and Closest_landmark_range are random access ranges. + * + * Procedure takes into account at most (dim+1) + * first landmarks from each landmark range to construct simplices. + * + * Landmarks are supposed to be in [0,nbL_-1] + */ + template <class KNearestNeighbors, class SimplicialComplexForWitness> + void witness_complex(KNearestNeighbors const & knn, + int nbL, + int dim, + SimplicialComplexForWitness & sc) { + Witness_complex<SimplicialComplexForWitness>(knn, nbL, dim, sc); + } + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // WITNESS_COMPLEX_H_ diff --git a/src/Witness_complex/test/CMakeLists.txt b/src/Witness_complex/test/CMakeLists.txt new file mode 100644 index 00000000..9422c152 --- /dev/null +++ b/src/Witness_complex/test/CMakeLists.txt @@ -0,0 +1,34 @@ +cmake_minimum_required(VERSION 2.6) +project(GUDHIWitnessComplexUT) + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg") +endif() + +add_executable ( simple_witness_complexUT simple_witness_complex.cpp ) +target_link_libraries(simple_witness_complexUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests definition and xml result file generation +add_test(NAME simple_witness_complexUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complexUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/simple_witness_complexUT.xml --log_level=test_suite --report_level=no) + +add_executable ( witness_complex_pointsUT witness_complex_points.cpp ) +target_link_libraries(witness_complex_pointsUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +# Unitary tests definition and xml result file generation +add_test(NAME witness_complex_pointsUT + COMMAND ${CMAKE_CURRENT_BINARY_DIR}/witness_complex_pointsUT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/witness_complex_pointsUT.xml --log_level=test_suite --report_level=no) + diff --git a/src/Witness_complex/test/simple_witness_complex.cpp b/src/Witness_complex/test/simple_witness_complex.cpp new file mode 100644 index 00000000..03df78ee --- /dev/null +++ b/src/Witness_complex/test/simple_witness_complex.cpp @@ -0,0 +1,59 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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_DYN_LINK +#define BOOST_TEST_MODULE "simple_witness_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> + +#include <iostream> +#include <ctime> +#include <vector> + +typedef Gudhi::Simplex_tree<> Simplex_tree; +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex; + +BOOST_AUTO_TEST_CASE(simple_witness_complex) { + Simplex_tree complex; + std::vector< typeVectorVertex > knn; + + knn.push_back({1, 0, 5, 2, 6, 3, 4}); + knn.push_back({2, 6, 4, 5, 0, 1, 3}); + knn.push_back({3, 4, 2, 1, 5, 6, 0}); + knn.push_back({4, 2, 1, 3, 5, 6, 0}); + knn.push_back({5, 1, 6, 0, 2, 3, 4}); + knn.push_back({6, 0, 5, 2, 1, 3, 4}); + knn.push_back({0, 5, 6, 1, 2, 3, 4}); + knn.push_back({2, 6, 4, 5, 3, 1, 0}); + knn.push_back({1, 2, 5, 4, 3, 6, 0}); + knn.push_back({3, 4, 0, 6, 5, 1, 2}); + knn.push_back({5, 0, 1, 3, 6, 2, 4}); + knn.push_back({5, 6, 1, 0, 2, 3, 4}); + knn.push_back({1, 6, 0, 5, 2, 3, 4}); + WitnessComplex witnessComplex(knn, 7, 7, complex); + + BOOST_CHECK(witnessComplex.is_witness_complex(knn, false)); +} diff --git a/src/Witness_complex/test/witness_complex_points.cpp b/src/Witness_complex/test/witness_complex_points.cpp new file mode 100644 index 00000000..bd3df604 --- /dev/null +++ b/src/Witness_complex/test/witness_complex_points.cpp @@ -0,0 +1,64 @@ +/* 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): Siargey Kachanovich + * + * Copyright (C) 2016 INRIA Sophia Antipolis-Méditerranée (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_DYN_LINK +#define BOOST_TEST_MODULE "witness_complex_points" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Landmark_choice_by_random_point.h> +#include <gudhi/Landmark_choice_by_furthest_point.h> + +#include <iostream> +#include <vector> + +typedef std::vector<double> Point; +typedef std::vector< Vertex_handle > typeVectorVertex; +typedef Gudhi::Simplex_tree<> Simplex_tree; +typedef Gudhi::witness_complex::Witness_complex<Simplex_tree> WitnessComplex; + +BOOST_AUTO_TEST_CASE(witness_complex_points) { + std::vector< typeVectorVertex > knn; + std::vector< Point > points; + // Add grid points as witnesses + for (double i = 0; i < 10; i += 1.0) + for (double j = 0; j < 10; j += 1.0) + for (double k = 0; k < 10; k += 1.0) + points.push_back(Point({i, j, k})); + + bool b_print_output = false; + // First test: random choice + Simplex_tree complex1; + Gudhi::witness_complex::landmark_choice_by_random_point(points, 100, knn); + assert(!knn.empty()); + WitnessComplex witnessComplex1(knn, 100, 3, complex1); + BOOST_CHECK(witnessComplex1.is_witness_complex(knn, b_print_output)); + + // Second test: furthest choice + knn.clear(); + Simplex_tree complex2; + Gudhi::witness_complex::landmark_choice_by_furthest_point(points, 100, knn); + WitnessComplex witnessComplex2(knn, 100, 3, complex2); + BOOST_CHECK(witnessComplex2.is_witness_complex(knn, b_print_output)); +} diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index 35313409..1db1ea8a 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -57,6 +57,22 @@ </td> </tr> </table> + \subsection WitnessComplexDataStructure Witness complex + \image html "Witness_complex_representation.png" "Witness complex representation" +<table border="0"> + <tr> + <td width="25%"> + <b>Introduced in:</b> GUDHI 1.3.0<br> + <b>Copyright:</b> GPL v3<br> + </td> + <td width="75%"> + <i>Siargey Kachanovich</i><br> + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$. + The data structure is described in \cite boissonnatmariasimplextreealgorithmica .<br> + <b>User manual:</b> \ref witness_complex - <b>Reference manual:</b> Gudhi::witness_complex::SimplicialComplexForWitness + </td> + </tr> +</table> \section Toolbox Toolbox \subsection PersistentCohomologyToolbox Persistent Cohomology |