diff options
Diffstat (limited to 'src/Witness_complex/include/gudhi/Witness_complex.h')
| -rw-r--r-- | src/Witness_complex/include/gudhi/Witness_complex.h | 414 |
1 files changed, 192 insertions, 222 deletions
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 819a0583..34cbc882 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -23,257 +23,227 @@ #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/container/flat_map.hpp> #include <boost/iterator/transform_iterator.hpp> + +#include <gudhi/distance_functions.h> + #include <algorithm> #include <utility> -#include <gudhi/distance_functions.h> #include <vector> #include <list> #include <set> #include <queue> #include <limits> -#include <math.h> #include <ctime> #include <iostream> -// 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> - namespace Gudhi { - namespace witness_complex { - - /** - \class Witness_complex - \brief Constructs the witness complex for the given set of witnesses and landmarks. - \ingroup witness_complex - */ - template< class Simplicial_complex> - class Witness_complex { - - private: - - struct Active_witness { - int witness_id; - int landmark_id; - - Active_witness(int witness_id_, int landmark_id_) +namespace witness_complex { + +/** + \class Witness_complex + \brief Constructs the witness complex for the given set of witnesses and landmarks. + \ingroup witness_complex + */ +template< class Simplicial_complex> +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_) - {} - }; + landmark_id(landmark_id_) { } + }; + private: + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - private: - - typedef typename Simplicial_complex::Simplex_handle Simplex_handle; - typedef typename Simplicial_complex::Vertex_handle Vertex_handle; - - typedef std::vector< double > Point_t; - typedef std::vector< Point_t > Point_Vector; + typedef std::vector< double > Point_t; + typedef std::vector< Point_t > Point_Vector; - //typedef typename Simplicial_complex::Filtration_value Filtration_value; - 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 - double density; // Desired density - Simplicial_complex& sc; // Simplicial complex - - public: + // typedef typename Simplicial_complex::Filtration_value Filtration_value; + typedef std::vector< Vertex_handle > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; + typedef std::pair< Simplex_handle, bool > typePairSimplexBool; - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /** @name Constructor - */ + typedef int Witness_id; + typedef int Landmark_id; + typedef std::list< Vertex_handle > ActiveWitnessList; - //@{ + private: + int nbL; // Number of landmarks + Simplicial_complex& sc; // Simplicial complex - // 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, - Simplicial_complex & sc_, - int nbL_, - int dim ): nbL(nbL_), sc(sc_) - { - //Construction of the active witness list - int nbW = knn.size(); - typeVectorVertex vv; - typeSimplex simplex; - typePairSimplexBool returnValue; - int counter = 0; - /* The list of still useful witnesses - * it will diminuish in the course of iterations - */ - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int 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}; - returnValue = sc.insert_simplex(vv); - /* TODO 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); - //std::cout << "Successfully added edges" << std::endl; - while (!active_w.empty() && k < dim ) - { - //std::cout << "Started the step k=" << k << std::endl; - 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]); - returnValue = sc.insert_simplex(simplex_vector,0.0); - it++; - } - else - active_w.erase(it++); //First increase the iterator and then erase the previous element - } - k++; - } - } + public: + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @name Constructor + */ - //@} - - 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 + //@{ + + // 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, + Simplicial_complex & sc_, + int nbL_, + int dim) : nbL(nbL_), sc(sc_) { + // Construction of the active witness list + int nbW = knn.size(); + typeVectorVertex vv; + int counter = 0; + /* The list of still useful witnesses + * it will diminuish in the course of iterations */ - template <typename KNearestNeighbors> - bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) - { - //std::cout << "All face in with the landmark " << inserted_vertex << std::endl; - std::vector< Vertex_handle > facet; - //Vertex_handle curr_vh = curr_sh->first; - // 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; - //std::cout << "++++ finished loop safely\n"; - } //endfor - return true; + ActiveWitnessList active_w; // = new ActiveWitnessList(); + for (int 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); + // std::cout << "Successfully added edges" << std::endl; + while (!active_w.empty() && k < dim) { + // std::cout << "Started the step k=" << k << std::endl; + 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, 0.0); + // 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++; } + } - template <typename T> - 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; - } + //@} + + 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::cout << "All face in with the landmark " << inserted_vertex << std::endl; + std::vector< Vertex_handle > facet; + // Vertex_handle curr_vh = curr_sh->first; + // 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; + // std::cout << "++++ finished loop safely\n"; + } // 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) - { - //bool final_result = true; - for (Simplex_handle sh: sc.complex_simplex_range()) - { - bool is_witnessed = false; - typeVectorVertex simplex; - int nbV = 0; //number of verticed in the simplex - for (int v: sc.simplex_vertex_range(sh)) - simplex.push_back(v); - nbV = simplex.size(); - for (typeVectorVertex w: knn) - { - bool has_vertices = true; - for (int v: simplex) - if (std::find(w.begin(), w.begin()+nbV, v) == w.begin()+nbV) - { - has_vertices = false; - //break; - } - 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; - } + 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) { + // bool final_result = true; + for (Simplex_handle sh : sc.complex_simplex_range()) { + bool is_witnessed = false; + typeVectorVertex simplex; + int nbV = 0; // number of verticed in the simplex + for (int v : sc.simplex_vertex_range(sh)) + simplex.push_back(v); + nbV = simplex.size(); + for (typeVectorVertex w : knn) { + bool has_vertices = true; + for (int v : simplex) + if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) { + has_vertices = false; + // break; + } + 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; } - return true; + assert(is_witnessed); + return false; + } } + return true; + } +}; - - }; //class Witness_complex +} // namespace witness_complex - } //namespace witness_complex - -} // namespace Guhdi +} // namespace Gudhi -#endif +#endif // WITNESS_COMPLEX_H_ |