diff options
Diffstat (limited to 'include/gudhi/Witness_complex.h')
| -rw-r--r-- | include/gudhi/Witness_complex.h | 265 |
1 files changed, 265 insertions, 0 deletions
diff --git a/include/gudhi/Witness_complex.h b/include/gudhi/Witness_complex.h new file mode 100644 index 00000000..489cdf11 --- /dev/null +++ b/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_ |