From e1eab27ee280826a5c9c44d530679626da2efb1f Mon Sep 17 00:00:00 2001 From: skachano Date: Thu, 29 Sep 2016 18:30:50 +0000 Subject: Weak Witness + Active Witness. It compiles git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1593 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7307e01148fc7cadfc22bd136ebdf0d55c8da937 --- .../example/witness_complex_sphere.cpp | 19 ++- .../include/gudhi/Active_witness/Active_witness.h | 70 ++++++++ .../gudhi/Active_witness/Active_witness_iterator.h | 96 +++++++++++ .../include/gudhi/Witness_complex.h | 188 +++++++++++---------- 4 files changed, 276 insertions(+), 97 deletions(-) create mode 100644 src/Witness_complex/include/gudhi/Active_witness/Active_witness.h create mode 100644 src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h (limited to 'src') diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp index 3d4a54aa..5a041c3e 100644 --- a/src/Witness_complex/example/witness_complex_sphere.cpp +++ b/src/Witness_complex/example/witness_complex_sphere.cpp @@ -54,6 +54,8 @@ void write_data(Data_range & data, std::string filename) { } int main(int argc, char * const argv[]) { + typedef Gudhi::witness_complex::Witness_complex> Witness_complex; + if (argc != 2) { std::cerr << "Usage: " << argv[0] << " number_of_landmarks \n"; @@ -77,21 +79,20 @@ int main(int argc, char * const argv[]) { // Choose landmarks start = clock(); - std::vector > knn; - Gudhi::subsampling::pick_n_random_points(point_vector, 100, std::back_inserter(landmarks)); - Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn); + Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks)); // Compute witness complex - // Gudhi::witness_complex::witness_complex(knn, number_of_landmarks, point_vector[0].size(), simplex_tree); - Gudhi::witness_complex::Witness_complex, CGAL::Epick_d>(landmarks.begin(), - landmarks.end(), - point_vector.begin(), - point_vector.end(), - simplex_tree); + Witness_complex witness_complex(landmarks.begin(), + landmarks.end(), + point_vector.begin(), + point_vector.end()); + witness_complex.create_complex(simplex_tree, 0); end = clock(); double time = static_cast(end - start) / CLOCKS_PER_SEC; std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " << time << " s. \n"; + assert(1 == 0); + std::cout << simplex_tree << "\n"; std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; l_time.push_back(std::make_pair(nbP, time)); } diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h new file mode 100644 index 00000000..87981c25 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness.h @@ -0,0 +1,70 @@ +/* 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 . + */ + +#ifndef ACTIVE_WITNESS_H_ +#define ACTIVE_WITNESS_H_ + +#include "Active_witness_iterator.h" +#include +#include + +namespace Gudhi { + +namespace witness_complex { + +template< typename Id_distance_pair, + typename INS_range > +class Active_witness { +public: + typedef Active_witness ActiveWitness; + typedef typename INS_range::iterator INS_iterator; + typedef Active_witness_iterator< ActiveWitness, Id_distance_pair, INS_iterator > iterator; + + std::vector nearest_landmark_table_; + INS_range search_range_; + INS_iterator iterator_last_; + INS_iterator iterator_end_; + + Active_witness(INS_range search_range) + : search_range_(search_range), iterator_end_(search_range.end()), iterator_last_(search_range.begin()) + { + nearest_landmark_table_.push_back(*iterator_last_); + } + + iterator begin() + { + return iterator(this, nearest_landmark_table_.begin()); + } + + iterator end() + { + return iterator(this); + } + + std::vector end_element_table_; + typename std::vector::iterator end_pointer = end_element_table_.end(); +}; + +} +} + +#endif diff --git a/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.h new file mode 100644 index 00000000..cd4f4b92 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Active_witness/Active_witness_iterator.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) 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 . + */ + +#ifndef ACTIVE_WITNESS_ITERATOR_H_ +#define ACTIVE_WITNESS_ITERATOR_H_ + +//#include "Active_witness.h" +#include +#include + +namespace Gudhi { + +namespace witness_complex { + +template< typename Active_witness, + typename Id_distance_pair, + typename INS_iterator > +class Active_witness_iterator + : public boost::iterator_facade< Active_witness_iterator +, Id_distance_pair const +, boost::forward_traversal_tag +, Id_distance_pair const> { + friend class boost::iterator_core_access; + + //typedef Active_witness Active_witness; + typedef typename std::vector::iterator Pair_iterator; + typedef typename Gudhi::witness_complex::Active_witness_iterator Iterator; + + + Active_witness *aw_; + Pair_iterator lh_; // landmark handle + //INS_iterator iterator_last; + //INS_iterator iterator_end; + +public: + Active_witness_iterator(Active_witness* aw) + : aw_(aw), lh_(aw_->end_pointer) + { + } + + Active_witness_iterator(Active_witness* aw, Pair_iterator lh_) + { + } + +private : + + Id_distance_pair dereference() const + { + return *lh_; + } + + bool equal(const Iterator& other) const + { + return (lh_ == other.lh_); + } + + void increment() + { + // if the id of the current landmark is the same as the last one + if (lh_->first == aw_->iterator_last_->first) { + // if the next iterator is end, lh_it = nullptr + if (++(aw_->iterator_last_) == aw_->iterator_end_) { + lh_ = aw_->end_pointer; + return; + } + else + aw_->nearest_landmark_table_.push_back(*(aw_->iterator_last_)); + } + lh_++; + } + +}; + +} +} + +#endif diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 4c489e7a..c6dff3bf 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -30,6 +30,7 @@ #include +#include "Active_witness/Active_witness.h" #include #include @@ -55,9 +56,9 @@ namespace witness_complex { // \brief Constructs the witness complex for the given set of witnesses and landmarks. // \ingroup witness_complex // */ -template< class SimplicialComplex, - class Kernel_ > +template< class Kernel_ > class Witness_complex { +private: typedef Kernel_ K; typedef typename K::Point_d Point_d; typedef typename K::FT FT; @@ -65,40 +66,35 @@ class Witness_complex { typedef gss::Kd_tree_search Kd_tree; typedef typename Kd_tree::INS_range Nearest_landmark_range; typedef typename std::vector Nearest_landmark_table; - typedef typename Nearest_landmark_table::const_iterator Nearest_landmark_table_iterator; + typedef typename Nearest_landmark_range::iterator Nearest_landmark_row_iterator; //typedef std::vector> Nearest_landmarks; - 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; + // 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_) { } + // }; + 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 FT Filtration_value; - typedef int Witness_id; - typedef int Landmark_id; - typedef std::list< Vertex_handle > ActiveWitnessList; + + typedef std::size_t Witness_id; + typedef std::size_t Landmark_id; + typedef std::pair Id_distance_pair; + typedef Active_witness ActiveWitness; + typedef std::list< ActiveWitness > ActiveWitnessList; + typedef std::vector< Landmark_id > typeVectorVertex; + typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex; private: - int nbL_; // Number of landmarks - SimplicialComplex& sc_; // Simplicial complex Point_range witnesses_, landmarks_; Kd_tree landmark_tree_; - std::vector nearest_landmarks_; public: ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -128,28 +124,38 @@ class Witness_complex { Witness_complex(InputIteratorLandmarks landmarks_first, InputIteratorLandmarks landmarks_last, InputIteratorWitnesses witnesses_first, - InputIteratorWitnesses witnesses_last, - SimplicialComplex& sc) - : witnesses_(witnesses_first, witnesses_last), landmarks_(landmarks_first, landmarks_last), landmark_tree_(landmarks_), sc_(sc) + InputIteratorWitnesses witnesses_last) + : witnesses_(witnesses_first, witnesses_last), landmarks_(landmarks_first, landmarks_last), landmark_tree_(landmarks_) { - for (Point_d w : witnesses_) - nearest_landmarks_.push_back(landmark_tree_.query_incremental_nearest_neighbors(w)); } + /** \brief Returns the point corresponding to the given vertex. + */ Point_d get_point( std::size_t vertex ) const { return landmarks_[vertex]; } + /** \brief Outputs the (weak) witness complex with + * squared relaxation parameter 'max_alpha_square' + * to simplicial complex 'complex'. + */ template < typename SimplicialComplexForWitness > bool create_complex(SimplicialComplexForWitness& complex, - FT max_alpha_square) const + FT max_alpha_square) { - unsigned nbW = witnesses_.size(), - nbL = landmarks_.size(); + unsigned nbL = landmarks_.size(); + if (complex.num_vertices() > 0) { + std::cerr << "Witness complex cannot create complex - complex is not empty.\n"; + return false; + } + if (max_alpha_square < 0) { + std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n"; + return false; + } typeVectorVertex vv; - ActiveWitnessList active_w;// = new ActiveWitnessList(); - for (int i = 0; i != nbL; ++i) { + ActiveWitnessList active_witnesses;// = new ActiveWitnessList(); + for (auto 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++; @@ -158,26 +164,27 @@ class Witness_complex { /* TODO Error if not inserted : normally no need here though*/ } unsigned k = 1; /* current dimension in iterative construction */ - for (int i=0; i != nbW; ++i) - active_w.push_back(i); - while (!active_w.empty() && k < nbL ) { - typename ActiveWitnessList::iterator aw_it = active_w.begin(); - while (aw_it != active_w.end()) { + for (auto w: witnesses_) + active_witnesses.push_back(ActiveWitness(landmark_tree_.query_incremental_nearest_neighbors(w))); + while (!active_witnesses.empty() && k < nbL ) { + typename ActiveWitnessList::iterator aw_it = active_witnesses.begin(); + while (aw_it != active_witnesses.end()) { std::vector simplex; bool ok = add_all_faces_of_dimension(k, max_alpha_square, std::numeric_limits::infinity(), - nearest_landmarks_[*aw_it].begin(), + aw_it->begin(), simplex, complex, - nearest_landmarks_[*aw_it].end()); + aw_it->end()); if (!ok) - active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + active_witnesses.erase(aw_it++); //First increase the iterator and then erase the previous element else aw_it++; } k++; } + return true; } //@} @@ -193,34 +200,37 @@ class Witness_complex { bool add_all_faces_of_dimension(int dim, double alpha2, double norelax_dist2, - Nearest_landmark_table_iterator curr_l, + typename ActiveWitness::iterator curr_l, std::vector& simplex, SimplicialComplexForWitness& sc, - Nearest_landmark_table_iterator end) + typename ActiveWitness::iterator end) { if (curr_l == end) return false; bool will_be_active = false; - Nearest_landmark_table_iterator l_it = curr_l; + typename ActiveWitness::iterator l_it = curr_l; if (dim > 0) for (; l_it->second - alpha2 <= norelax_dist2 && l_it != end; ++l_it) { simplex.push_back(l_it->first); - if (sc.find(simplex) != sc.null_simplex()) + 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, - l_it+1, + ++next_it, simplex, sc, end) || will_be_active; + } simplex.pop_back(); // If norelax_dist is infinity, change to first omitted distance if (l_it->second <= norelax_dist2) norelax_dist2 = l_it->second; + typename ActiveWitness::iterator next_it = l_it; will_be_active = add_all_faces_of_dimension(dim-1, alpha2, norelax_dist2, - l_it+1, + ++next_it, simplex, sc, end) || will_be_active; @@ -253,6 +263,8 @@ class Witness_complex { double* filtration_value, SimplicialComplexForWitness& sc) { + typedef typename SimplicialComplexForWitness::Simplex_handle Simplex_handle; + std::vector< int > facet; for (std::vector::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) { @@ -297,45 +309,45 @@ class Witness_complex { // * \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; - } + // 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; + // } }; /** -- cgit v1.2.3