diff options
Diffstat (limited to 'include/gudhi/Witness_complex.h')
| -rw-r--r-- | include/gudhi/Witness_complex.h | 333 |
1 files changed, 137 insertions, 196 deletions
diff --git a/include/gudhi/Witness_complex.h b/include/gudhi/Witness_complex.h index 489cdf11..63f03687 100644 --- a/include/gudhi/Witness_complex.h +++ b/include/gudhi/Witness_complex.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2015 INRIA (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 @@ -23,65 +23,44 @@ #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 <gudhi/Active_witness/Active_witness.h> +#include <gudhi/Witness_complex/all_faces_in.h> -#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> +/** + * \private + * \class Witness_complex + * \brief Constructs (weak) witness complex for a given table of nearest landmarks with respect to witnesses. + * \ingroup witness_complex + * + * \tparam Nearest_landmark_table_ needs to be a range of a range of pairs of nearest landmarks and distances. + * The class Nearest_landmark_table_::value_type must be a copiable range. + * The range of pairs must admit a member type 'iterator'. The dereference type + * of the pair range iterator needs to be 'std::pair<std::size_t, double>'. +*/ +template< class Nearest_landmark_table_ > 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 + typedef typename Nearest_landmark_table_::value_type Nearest_landmark_range; + typedef std::size_t Witness_id; + typedef std::size_t Landmark_id; + typedef std::pair<Landmark_id, double> Id_distance_pair; + typedef Active_witness<Id_distance_pair, Nearest_landmark_range> ActiveWitness; + typedef std::list< ActiveWitness > ActiveWitnessList; + typedef std::vector< Landmark_id > typeVectorVertex; + typedef std::vector<Nearest_landmark_range> Nearest_landmark_table_internal; + typedef Landmark_id Vertex_handle; + + protected: + Nearest_landmark_table_internal nearest_landmark_table_; public: ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -90,174 +69,136 @@ class Witness_complex { //@{ - // Witness_range<Closest_landmark_range<Vertex_handle>> + Witness_complex() { + } - /* - * \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] + /** + * \brief Initializes member variables before constructing simplicial complex. + * \details Records nearest landmark table. + * @param[in] nearest_landmark_table needs to be a range of a range of pairs of nearest landmarks and distances. + * The class Nearest_landmark_table_::value_type must be a copiable range. + * The range of pairs must admit a member type 'iterator'. The dereference type + * of the pair range iterator needs to be 'std::pair<std::size_t, double>'. */ - 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 + + Witness_complex(Nearest_landmark_table_ const & nearest_landmark_table) + : nearest_landmark_table_(std::begin(nearest_landmark_table), std::end(nearest_landmark_table)) { + } + + /** \brief Outputs the (weak) witness complex of relaxation 'max_alpha_square' + * in a simplicial complex data structure. + * \details The function returns true if the construction is successful and false otherwise. + * @param[out] complex Simplicial complex data structure compatible which is a model of + * SimplicialComplexForWitness concept. + * @param[in] max_alpha_square Maximal squared relaxation parameter. + * @param[in] limit_dimension Represents the maximal dimension of the simplicial complex + * (default value = no limit). + */ + template < typename SimplicialComplexForWitness > + bool create_complex(SimplicialComplexForWitness& complex, + double max_alpha_square, + std::size_t limit_dimension = std::numeric_limits<std::size_t>::max()) const { + if (complex.num_vertices() > 0) { + std::cerr << "Witness complex cannot create complex - complex is not empty.\n"; + return false; } - 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 - } + if (max_alpha_square < 0) { + std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n"; + return false; + } + if (limit_dimension < 0) { + std::cerr << "Witness complex cannot create complex - limit dimension must be non-negative.\n"; + return false; + } + ActiveWitnessList active_witnesses; + Landmark_id k = 0; /* current dimension in iterative construction */ + for (auto w : nearest_landmark_table_) + active_witnesses.push_back(ActiveWitness(w)); + while (!active_witnesses.empty() && k <= limit_dimension) { + typename ActiveWitnessList::iterator aw_it = active_witnesses.begin(); + std::vector<Landmark_id> simplex; + simplex.reserve(k+1); + while (aw_it != active_witnesses.end()) { + bool ok = add_all_faces_of_dimension(k, + max_alpha_square, + std::numeric_limits<double>::infinity(), + aw_it->begin(), + simplex, + complex, + aw_it->end()); + assert(simplex.empty()); + if (!ok) + active_witnesses.erase(aw_it++); // First increase the iterator and then erase the previous element + else + aw_it++; } k++; } + complex.set_dimension(k-1); + return true; } //@} 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 + /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness. + * Iterator is needed to know until how far we can take landmarks to form simplexes. + * simplex is the prefix of the simplexes to insert. + * The output value indicates if the witness rests active or not. */ - 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]); + template < typename SimplicialComplexForWitness > + bool add_all_faces_of_dimension(int dim, + double alpha2, + double norelax_dist2, + typename ActiveWitness::iterator curr_l, + std::vector<Landmark_id>& simplex, + SimplicialComplexForWitness& sc, + typename ActiveWitness::iterator end) const { + if (curr_l == end) + return false; + bool will_be_active = false; + typename ActiveWitness::iterator l_it = curr_l; + if (dim > 0) { + for (; l_it != end && l_it->second - alpha2 <= norelax_dist2; ++l_it) { + simplex.push_back(l_it->first); + if (sc.find(simplex) != sc.null_simplex()) { + typename ActiveWitness::iterator next_it = l_it; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha2, + norelax_dist2, + ++next_it, + simplex, + sc, + end) || will_be_active; } - } // 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; + assert(!simplex.empty()); + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (l_it->second <= norelax_dist2) + norelax_dist2 = l_it->second; } - } - 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; + } else if (dim == 0) { + for (; l_it != end && l_it->second - alpha2 <= norelax_dist2; ++l_it) { + simplex.push_back(l_it->first); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (l_it->second > norelax_dist2) + filtration_value = l_it->second - norelax_dist2; + if (all_faces_in(simplex, &filtration_value, sc)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); } - assert(is_witnessed); - return false; + assert(!simplex.empty()); + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (l_it->second < norelax_dist2) + norelax_dist2 = l_it->second; } } - return true; + return will_be_active; } }; - /** - * \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 |