diff options
Diffstat (limited to 'src/Witness_complex/include')
| -rw-r--r-- | src/Witness_complex/include/gudhi/A0_complex.h | 337 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Dim_list_iterator.h | 155 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Dim_lists.h | 195 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Good_links.h | 449 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Relaxed_witness_complex.h | 379 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Strange_witness_complex.h | 232 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h | 337 | ||||
| -rw-r--r-- | src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h | 379 |
8 files changed, 2463 insertions, 0 deletions
diff --git a/src/Witness_complex/include/gudhi/A0_complex.h b/src/Witness_complex/include/gudhi/A0_complex.h new file mode 100644 index 00000000..2018e1c8 --- /dev/null +++ b/src/Witness_complex/include/gudhi/A0_complex.h @@ -0,0 +1,337 @@ +/* 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 A0_COMPLEX_H_ +#define A0_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#include <CGAL/Cartesian_d.h> +#include <CGAL/Search_traits.h> +#include <CGAL/Search_traits_adapter.h> +#include <CGAL/property_map.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Orthogonal_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class A0_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 Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + 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 + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * 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 KNearestNeighbours > + A0_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha2, + unsigned limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + 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 (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, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + for (int i=0; i != nbW; ++i) { + // int i_end = limD+1; + // if (knn[i].size() < limD+1) + // i_end = knn[i].size(); + // double dist_wL = *(distances[i].begin()); + // while (distances[i][i_end] > dist_wL + alpha2) + // i_end--; + // add_all_witnessed_faces(distances[i].begin(), + // knn[i].begin(), + // knn[i].begin() + i_end + 1); + unsigned j_end = 0; + while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { + std::vector<int> simplex; + for (unsigned j = 0; j <= j_end; ++j) + simplex.push_back(knn[i][j]); + assert(distances[i][j_end] - distances[i][0] >= 0); + sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); + j_end++; + } + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \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 + */ + void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>::const_iterator end) + { + std::vector<int> simplex; + std::vector<int>::const_iterator l_end = curr_l; + for (; l_end != end; ++l_end) { + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + simplex = {}; + for (; l_it != l_end; ++l_it, ++d_it) + simplex.push_back(*l_it); + sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); + } + } + + /** \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 + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Dim_list_iterator.h b/src/Witness_complex/include/gudhi/Dim_list_iterator.h new file mode 100644 index 00000000..3f23e8c9 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Dim_list_iterator.h @@ -0,0 +1,155 @@ +/* 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 DIM_LISTS_ITERATOR_H_ +#define DIM_LISTS_ITERATOR_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplex_handle > +class Dim_lists_iterator { + +private: + + typedef typename std::list<Simplex_handle>::iterator List_iterator; + typedef typename std::vector<std::list<Simplex_handle>>::reverse_iterator Vector_iterator; + typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator; + + + List_iterator sh_; + Vector_iterator curr_list_; + typename std::vector<std::list<Simplex_handle>>& table_; + +public: + + Dim_lists_iterator(List_iterator sh, + Vector_iterator curr_list, + typename std::vector<std::list<Simplex_handle>>& table) + : sh_(sh), curr_list_(curr_list), table_(table) + { + } + + Simplex_handle operator*() + { + return *sh_; + } + + Iterator operator++() + { + increment(); + return (*this); + } + + Iterator operator++(int) + { + Iterator prev_it(sh_, curr_list_, table_); + increment(); + return prev_it; + } + + Iterator dim_begin() + { + return Iterator(curr_list_->begin(), curr_list_, table_); + } + + Iterator dim_end() + { + return Iterator(curr_list_->end(), curr_list_, table_); + } + + Iterator dimp1_begin() + { + return Iterator((curr_list_-1)->begin(), curr_list_-1, table_); + } + + Iterator dimp1_end() + { + return Iterator((curr_list_-1)->end(), curr_list_-1, table_); + } + + bool operator==(const Iterator& it2) const + { + return (sh_ == it2.sh_); + } + + bool operator!=(const Iterator& it2) const + { + return (sh_ != it2.sh_); + } + + void remove_incr() + { + + } + +private: + + void increment() + { + if (++sh_ == curr_list_->end()) + if (++curr_list_ != table_.rend()) + sh_ = curr_list_->begin(); + // The iterator of the end of the table is the end of the last list + } + + +}; //class Dim_lists_iterator + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Dim_lists.h b/src/Witness_complex/include/gudhi/Dim_lists.h new file mode 100644 index 00000000..1d1b03c3 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Dim_lists.h @@ -0,0 +1,195 @@ +/* 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 DIM_LISTS_H_ +#define DIM_LISTS_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include <gudhi/Dim_list_iterator.h> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Dim_lists { + +private: + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator; + + std::vector<std::list<Simplex_handle>> table_; + Simplicial_complex& sc_; + +public: + + Dim_lists(Simplicial_complex & sc, int dim, double alpha_max = 100) + : sc_(sc) + { + table_ = std::vector<std::list<Simplex_handle>>(dim+1); + for (auto sh: sc.filtration_simplex_range()) { + if (sc_.filtration(sh) < alpha_max) + table_[sc.dimension(sh)].push_front(sh); + } + auto t_it = table_.rbegin(); + while (t_it->empty()) { + t_it++; + table_.pop_back(); + } + } + + Iterator begin() + { + return Iterator(table_.rbegin()->begin(), table_.rbegin(), table_); + } + + Iterator end() + { + return Iterator(table_[0].end(), table_.rend(), table_); + } + + unsigned size() + { + unsigned curr_size = 0; + for (auto l: table_) + curr_size += l.size(); + return curr_size; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc_.simplex_vertex_range(face); + auto cfvr = sc_.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + void output_simplices() { + std::cout << "Size of vector: " << size() << std::endl; + for (auto line_it = table_.rbegin(); line_it != table_.rend(); ++line_it) + for (auto sh: *line_it) { + std::cout << sc_.dimension(sh) << " "; + for (auto v : sc_.simplex_vertex_range(sh)) + std::cout << v << " "; + std::cout << sc_.filtration(sh) << "\n"; + } + } + + void collapse() + { + auto coface_list_it = table_.rbegin(); + auto face_list_it = table_.rbegin()+1; + for ( ; + face_list_it != table_.rend(); + ++face_list_it) { + auto face_it = face_list_it->begin(); + while (face_it != face_list_it->end() && sc_.filtration(*face_it) != 0) { + int coface_count = 0; + auto reduced_coface = coface_list_it->begin(); + for (auto coface_it = coface_list_it->begin(); coface_it != coface_list_it->end() && sc_.filtration(*coface_it) != 0; ++coface_it) + if (is_face(*face_it, *coface_it)) { + coface_count++; + if (coface_count == 1) + reduced_coface = coface_it; + else + break; + } + if (coface_count == 1) { + /* + std::cout << "Erase ( "; + for (auto v: sc_.simplex_vertex_range(*reduced_coface)) + std::cout << v << " "; + */ + coface_list_it->erase(reduced_coface); + /* + std::cout << ") and then ( "; + for (auto v: sc_.simplex_vertex_range(*face_it)) + std::cout << v << " "; + std::cout << ")\n"; + */ + face_list_it->erase(face_it); + face_it = face_list_it->begin(); + } + else + face_it++; + } + if ((coface_list_it++)->empty()) + table_.pop_back(); + } + } + + bool is_pseudomanifold() + { + + return true; + } + +}; //class Dim_lists + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Good_links.h b/src/Witness_complex/include/gudhi/Good_links.h new file mode 100644 index 00000000..c5e993e5 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Good_links.h @@ -0,0 +1,449 @@ +/* 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 GOOD_LINKS_H_ +#define GOOD_LINKS_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include <gudhi/Dim_list_iterator.h> +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + +template < typename Simplicial_complex > +class Good_links { + + typedef typename Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef std::vector<Vertex_handle> Vertex_vector; + + // graph is an adjacency list + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph; + // map that gives to a certain simplex its node in graph and its dimension + //typedef std::pair<boost::vecS,Color> Reference; + typedef boost::graph_traits<Adj_graph>::vertex_descriptor Vertex_t; + typedef boost::graph_traits<Adj_graph>::edge_descriptor Edge_t; + typedef boost::graph_traits<Adj_graph>::adjacency_iterator Adj_it; + typedef std::pair<Adj_it, Adj_it> Out_edge_it; + + typedef boost::container::flat_map<Simplex_handle, Vertex_t> Graph_map; + typedef boost::container::flat_map<Vertex_t, Simplex_handle> Inv_graph_map; + +public: + Good_links(Simplicial_complex& sc): sc_(sc) + { + int dim = 0; + for (auto sh: sc_.complex_simplex_range()) + if (sc_.dimension(sh) > dim) + dim = sc_.dimension(sh); + dimension = dim; + count_good = Vertex_vector(dim); + count_bad = Vertex_vector(dim); + } + +private: + + Simplicial_complex& sc_; + unsigned dimension; + std::vector<int> count_good; + std::vector<int> count_bad; + + void add_vertices_to_link_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_vertices_to_link_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + curr_simplex, + curr_d+1, link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + + void add_edges_to_link_graph(Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map) + { + Simplex_handle sh; + // Add edges + //std::cout << "Entered add edges:\n"; + typename Graph_map::iterator map_it; + for (auto d_map_pair : d_map) { + //std::cout << "*"; + sh = d_map_pair.first; + Vertex_t d_vert = d_map_pair.second; + for (auto facet_sh : sc_.boundary_simplex_range(sh)) + //for (auto f_map_it : f_map) + { + //std::cout << "'"; + map_it = f_map.find(facet_sh); + //We must have all the facets in the graph at this point + assert(map_it != f_map.end()); + Vertex_t f_vert = map_it->second; + //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n"; + boost::add_edge(d_vert,f_vert,adj_graph); + } + } + } + + + + /* \brief Verifies if the simplices formed by vertices given by link_vertices + * form a pseudomanifold. + * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional + * faces and edges represent adjacency between them. + */ + bool link_is_pseudomanifold(Vertex_vector& star_vertices, + int dimension) + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Vertex_vector init_vector = {}; + add_vertices_to_link_graph(star_vertices, + star_vertices.begin()+1, + adj_graph, + d_map, + f_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2){ + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + int star_dim(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + int curr_d, + Vertex_vector& curr_simplex, + typename std::vector< int >::iterator curr_dc) + { + //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n"; + Simplex_handle sh; + int final_d = curr_d; + typename Vertex_vector::iterator it; + typename Vertex_vector::iterator dc_it; + //std::cout << "Current vertex is " << + for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it) + { + curr_simplex.push_back(*it); + Vertex_vector curr_simplex_copy(curr_simplex); + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " final_dim " << final_d; + */ + sh = sc_.find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first + if (sh != sc_.null_simplex()) + { + //std::cout << " -> " << *it << "\n"; + int d = star_dim(star_vertices, + it+1, + curr_d+1, + curr_simplex, + dc_it); + if (d >= final_d) + { + final_d = d; + //std::cout << d << " "; + //print_vector(curr_simplex); + //std::cout << std::endl; + } + if (d >= *dc_it) + *dc_it = d; + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); + } + return final_d; + } + +public: + + /** \brief Returns true if the link is good + */ + bool has_good_link(Vertex_handle v) + { + typedef Vertex_vector typeVectorVertex; + typeVectorVertex star_vertices; + // Fill star_vertices + star_vertices.push_back(v); + for (auto u: sc_.complex_vertex_range()) + { + typeVectorVertex edge = {u,v}; + if (u != v && sc_.find(edge) != sc_.null_simplex()) + star_vertices.push_back(u); + } + // Find the dimension + typeVectorVertex init_simplex = {star_vertices[0]}; + bool is_pure = true; + std::vector<int> dim_coface(star_vertices.size(), 1); + int d = star_dim(star_vertices, + star_vertices.begin()+1, + 0, + init_simplex, + dim_coface.begin()+1) - 1; //link_dim = star_dim - 1 + + assert(init_simplex.size() == 1); + // if (!is_pure) + // std::cout << "Found an impure star around " << v << "\n"; + for (int dc: dim_coface) + is_pure = (dc == dim_coface[0]); + /* + if (d == count_good.size()) + { + std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is "; + print_vector(star_vertices); std::cout << std::endl; + } + */ + //if (d == -1) count_bad[0]++; + bool b= (is_pure && link_is_pseudomanifold(star_vertices,d)); + if (d != -1) {if (b) count_good[d]++; else count_bad[d]++;} + if (!is_pure) count_bad[0]++; + return (d != -1 && b && is_pure); + } + +private: + void add_max_simplices_to_graph(Vertex_vector& star_vertices, + typename Vertex_vector::iterator curr_v, + Adj_graph& adj_graph, + Graph_map& d_map, + Graph_map& f_map, + Inv_graph_map& inv_d_map, + Vertex_vector& curr_simplex, + int curr_d, + int link_dimension) + { + Simplex_handle sh; + Vertex_t vert; + typename Vertex_vector::iterator it; + //std::pair<typename Graph_map::iterator,bool> resPair; + //typename Graph_map::iterator resPair; + //Add vertices + //std::cout << "Entered add vertices\n"; + for (it = curr_v; it != star_vertices.end(); ++it) { + curr_simplex.push_back(*it); //push next vertex in question + //curr_simplex.push_back(star_vertices[0]); //push the center of the star + /* + std::cout << "Searching for "; + print_vector(curr_simplex); + std::cout << " curr_dim " << curr_d << " d " << dimension << ""; + */ + Vertex_vector curr_simplex_copy(curr_simplex); + sh = sc_.find(curr_simplex_copy); //a simplex of the star + //curr_simplex.pop_back(); //pop the center of the star + curr_simplex_copy = Vertex_vector(curr_simplex); + if (sh != sc_.null_simplex()) { + //std::cout << " added\n"; + if (curr_d == link_dimension) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); //ASSERT! + vert = boost::add_vertex(adj_graph); + d_map.emplace(sh,vert); + inv_d_map.emplace(vert,sh); + } + else { + if (curr_d == link_dimension-1) { + sh = sc_.find(curr_simplex_copy); //a simplex of the link + assert(sh != sc_.null_simplex()); + vert = boost::add_vertex(adj_graph); + f_map.emplace(sh,vert); + } + //delete (&curr_simplex_copy); //Just so it doesn't stack + add_max_simplices_to_graph(star_vertices, + it+1, + adj_graph, + d_map, + f_map, + inv_d_map, + curr_simplex, + curr_d+1, + link_dimension); + } + } + /* + else + std::cout << "\n"; + */ + curr_simplex.pop_back(); //pop the vertex in question + } + } + +public: + bool complex_is_pseudomanifold() + { + Adj_graph adj_graph; + Graph_map d_map, f_map; + // d_map = map for d-dimensional simplices + // f_map = map for its facets + Inv_graph_map inv_d_map; + Vertex_vector init_vector = {}; + std::vector<int> star_vertices; + for (int v: sc_.complex_vertex_range()) + star_vertices.push_back(v); + add_max_simplices_to_graph(star_vertices, + star_vertices.begin(), + adj_graph, + d_map, + f_map, + inv_d_map, + init_vector, + 0, dimension); + //std::cout << "DMAP_SIZE: " << d_map.size() << "\n"; + //std::cout << "FMAP_SIZE: " << f_map.size() << "\n"; + add_edges_to_link_graph(adj_graph, + d_map, + f_map); + for (auto f_map_it : f_map) { + //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n"; + if (boost::out_degree(f_map_it.second, adj_graph) != 2) { + // if (boost::out_degree(f_map_it.second, adj_graph) >= 3) { + // // std::cout << "This simplex has 3+ cofaces: "; + // // for(auto v : sc_.simplex_vertex_range(f_map_it.first)) + // // std::cout << v << " "; + // // std::cout << std::endl; + // Adj_it ai, ai_end; + // for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) { + // auto it = inv_d_map.find(*ai); + // assert (it != inv_d_map.end()); + // typename Simplicial_complex::Simplex_handle sh = it->second; + // for(auto v : sc_.simplex_vertex_range(sh)) + // std::cout << v << " "; + // std::cout << std::endl; + // } + // } + count_bad[dimension]++; + return false; + } + } + // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices + // What is left is to check the connexity + //std::vector<int> components(boost::num_vertices(adj_graph)); + return true; //Forget the connexity + //return (boost::connected_components(adj_graph, &components[0]) == 1); + } + + int number_good_links(int dim) + { + return count_good[dim]; + } + + int number_bad_links(int dim) + { + return count_bad[dim]; + } + +}; + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h new file mode 100644 index 00000000..2971149a --- /dev/null +++ b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h @@ -0,0 +1,379 @@ +/* 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 RELAXED_WITNESS_COMPLEX_H_ +#define RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#include <CGAL/Cartesian_d.h> +#include <CGAL/Search_traits.h> +#include <CGAL/Search_traits_adapter.h> +#include <CGAL/property_map.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Orthogonal_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Relaxed_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 Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + 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 + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * 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 KNearestNeighbours > + Relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha, + int limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + 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 (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, Filtration_value(0.0)); + /* 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); + while (!active_w.empty() && k <= limD && k < nbL ) + { + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) + { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, + alpha, + std::numeric_limits<double>::infinity(), + distances[*aw_it].begin(), + knn[*aw_it].begin(), + simplex, + knn[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; + } + k++; + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \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 + */ + bool add_all_faces_of_dimension(int dim, + double alpha, + double norelax_dist, + std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>& simplex, + std::vector<int>::const_iterator end) + { + if (curr_l == end) + return false; + bool will_be_active = false; + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + if (dim > 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + if (sc.find(simplex) != sc.null_simplex()) + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it <= norelax_dist) + norelax_dist = *d_it; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + } + else if (dim == 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (*d_it > norelax_dist) + filtration_value = *d_it - norelax_dist; + if (all_faces_in(simplex, &filtration_value)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); + } + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it < norelax_dist) + norelax_dist = *d_it; + } + return will_be_active; + } + + /** \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 + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Strange_witness_complex.h b/src/Witness_complex/include/gudhi/Strange_witness_complex.h new file mode 100644 index 00000000..c1c2912b --- /dev/null +++ b/src/Witness_complex/include/gudhi/Strange_witness_complex.h @@ -0,0 +1,232 @@ +/* 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/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.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 { + +/** + \class Witness_complex + \brief Constructs the witness complex for the given set of witnesses and landmarks. + \ingroup witness_complex + */ +template< class Simplicial_complex> +class Strange_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 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 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 + Simplicial_complex& 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 > + Strange_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 + */ + 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 + } + for (int i = 0; i != nbW; ++i) { + std::vector<int> simplex; + simplex.reserve(dim+1); + for (int j = 0; j <= dim; j++) + simplex.push_back(knn[i][j]); + sc.insert_simplex_and_subfaces(simplex); + } + } + + //@} + + 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; + } + } + return true; + } +}; + +} // namespace witness_complex + +} // namespace Gudhi + +#endif // WITNESS_COMPLEX_H_ diff --git a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h new file mode 100644 index 00000000..ee863a42 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h @@ -0,0 +1,337 @@ +/* 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 STRONG_RELAXED_WITNESS_COMPLEX_H_ +#define STRONG_RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#include <CGAL/Cartesian_d.h> +#include <CGAL/Search_traits.h> +#include <CGAL/Search_traits_adapter.h> +#include <CGAL/property_map.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Orthogonal_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Strong_relaxed_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 Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + 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 + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * 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 KNearestNeighbours > + Strong_relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha2, + unsigned limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + 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 (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, Filtration_value(0.0)); + /* TODO Error if not inserted : normally no need here though*/ + } + for (int i=0; i != nbW; ++i) { + // int i_end = limD+1; + // if (knn[i].size() < limD+1) + // i_end = knn[i].size(); + // double dist_wL = *(distances[i].begin()); + // while (distances[i][i_end] > dist_wL + alpha2) + // i_end--; + // add_all_witnessed_faces(distances[i].begin(), + // knn[i].begin(), + // knn[i].begin() + i_end + 1); + unsigned j_end = 0; + while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) { + std::vector<int> simplex; + for (unsigned j = 0; j <= j_end; ++j) + simplex.push_back(knn[i][j]); + assert(distances[i][j_end] - distances[i][0] >= 0); + sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]); + j_end++; + } + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \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 + */ + void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>::const_iterator end) + { + std::vector<int> simplex; + std::vector<int>::const_iterator l_end = curr_l; + for (; l_end != end; ++l_end) { + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + simplex = {}; + for (; l_it != l_end; ++l_it, ++d_it) + simplex.push_back(*l_it); + sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d); + } + } + + /** \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 + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Strong_relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif diff --git a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h new file mode 100644 index 00000000..a1aedd8e --- /dev/null +++ b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h @@ -0,0 +1,379 @@ +/* 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 WEAK_RELAXED_WITNESS_COMPLEX_H_ +#define WEAK_RELAXED_WITNESS_COMPLEX_H_ + +#include <boost/container/flat_map.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <algorithm> +#include <utility> +#include "gudhi/reader_utils.h" +#include "gudhi/distance_functions.h" +#include "gudhi/Simplex_tree.h" +#include <vector> +#include <list> +#include <set> +#include <queue> +#include <limits> +#include <math.h> +#include <ctime> +#include <iostream> + +// Needed for nearest neighbours +#include <CGAL/Cartesian_d.h> +#include <CGAL/Search_traits.h> +#include <CGAL/Search_traits_adapter.h> +#include <CGAL/property_map.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Orthogonal_k_neighbor_search.h> + +#include <boost/tuple/tuple.hpp> +#include <boost/iterator/zip_iterator.hpp> +#include <boost/iterator/counting_iterator.hpp> +#include <boost/range/iterator_range.hpp> + +// 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 { + + /** \addtogroup simplex_tree + * Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$: + * \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$ + * and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that + * w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. + */ +template< class Simplicial_complex > +class Weak_relaxed_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 Simplicial_complex::Simplex_handle Simplex_handle; + typedef typename Simplicial_complex::Vertex_handle Vertex_handle; + typedef typename Simplicial_complex::Filtration_value FT; + + 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 + Simplicial_complex& sc; // Simplicial complex + + public: + /** @name Constructor + */ + + //@{ + + /** + * \brief Iterative construction of the relaxed witness complex. + * \details The witness complex is written in sc_ basing on a matrix knn + * of k nearest neighbours of the form {witnesses}x{landmarks} and + * and a matrix distances of distances to these landmarks from witnesses. + * The parameter alpha defines relaxation and + * limD defines the + * + * The line lengths in one given matrix can differ, + * however both matrices have the same corresponding line lengths. + * + * 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 KNearestNeighbours > + Weak_relaxed_witness_complex(std::vector< std::vector<double> > const & distances, + KNearestNeighbours const & knn, + Simplicial_complex & sc_, + int nbL_, + double alpha, + int limD) : nbL(nbL_), sc(sc_) { + int nbW = knn.size(); + 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 (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, Filtration_value(0.0)); + /* 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); + while (!active_w.empty() && k <= limD && k < nbL ) + { + typename ActiveWitnessList::iterator aw_it = active_w.begin(); + while (aw_it != active_w.end()) + { + std::vector<int> simplex; + bool ok = add_all_faces_of_dimension(k, + alpha, + std::numeric_limits<double>::infinity(), + distances[*aw_it].begin(), + knn[*aw_it].begin(), + simplex, + knn[*aw_it].end()); + if (!ok) + active_w.erase(aw_it++); //First increase the iterator and then erase the previous element + else + aw_it++; + } + k++; + } + sc.set_dimension(limD); + } + + //@} + +private: + /* \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 + */ + bool add_all_faces_of_dimension(int dim, + double alpha, + double norelax_dist, + std::vector<double>::const_iterator curr_d, + std::vector<int>::const_iterator curr_l, + std::vector<int>& simplex, + std::vector<int>::const_iterator end) + { + if (curr_l == end) + return false; + bool will_be_active = false; + std::vector<int>::const_iterator l_it = curr_l; + std::vector<double>::const_iterator d_it = curr_d; + if (dim > 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + if (sc.find(simplex) != sc.null_simplex()) + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it <= norelax_dist) + norelax_dist = *d_it; + will_be_active = add_all_faces_of_dimension(dim-1, + alpha, + norelax_dist, + d_it+1, + l_it+1, + simplex, + end) || will_be_active; + } + else if (dim == 0) + for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) { + simplex.push_back(*l_it); + double filtration_value = 0; + // if norelax_dist is infinite, relaxation is 0. + if (*d_it > norelax_dist) + filtration_value = *d_it - norelax_dist; + if (all_faces_in(simplex, &filtration_value)) { + will_be_active = true; + sc.insert_simplex(simplex, filtration_value); + } + simplex.pop_back(); + // If norelax_dist is infinity, change to first omitted distance + if (*d_it < norelax_dist) + norelax_dist = *d_it; + } + return will_be_active; + } + + /** \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 + */ + bool all_faces_in(std::vector<int>& simplex, double* filtration_value) + { + std::vector< int > facet; + for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it) + { + facet.clear(); + for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it) + if (it != not_it) + facet.push_back(*it); + Simplex_handle facet_sh = sc.find(facet); + if (facet_sh == sc.null_simplex()) + return false; + else if (sc.filtration(facet_sh) > *filtration_value) + *filtration_value = sc.filtration(facet_sh); + } + return true; + } + + bool is_face(Simplex_handle face, Simplex_handle coface) + { + // vertex range is sorted in decreasing order + auto fvr = sc.simplex_vertex_range(face); + auto cfvr = sc.simplex_vertex_range(coface); + auto fv_it = fvr.begin(); + auto cfv_it = cfvr.begin(); + while (fv_it != fvr.end() && cfv_it != cfvr.end()) { + if (*fv_it < *cfv_it) + ++cfv_it; + else if (*fv_it == *cfv_it) { + ++cfv_it; + ++fv_it; + } + else + return false; + + } + return (fv_it == fvr.end()); + } + + // void erase_simplex(Simplex_handle sh) + // { + // auto siblings = sc.self_siblings(sh); + // auto oncles = siblings->oncles(); + // int prev_vertex = siblings->parent(); + // siblings->members().erase(sh->first); + // if (siblings->members().empty()) { + // typename typedef Simplicial_complex::Siblings Siblings; + // oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex)); + // assert(!sc.has_children(oncles->members().find(prev_vertex))); + // //delete siblings; + // } + + // } + + void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh) + { + erase_simplex(coface_sh); + erase_simplex(face_sh); + } + +public: + // void collapse(std::vector<Simplex_handle>& simplices) + // { + // // Get a vector of simplex handles ordered by filtration value + // std::cout << sc << std::endl; + // //std::vector<Simplex_handle> simplices; + // for (Simplex_handle sh: sc.filtration_simplex_range()) + // simplices.push_back(sh); + // // std::sort(simplices.begin(), + // // simplices.end(), + // // [&](Simplex_handle sh1, Simplex_handle sh2) + // // { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2); + // // return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); }); + // // Double iteration + // auto face_it = simplices.rbegin(); + // while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) { + // int coface_count = 0; + // auto reduced_coface = simplices.rbegin(); + // for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it) + // if (face_it != coface_it && is_face(*face_it, *coface_it)) { + // coface_count++; + // if (coface_count == 1) + // reduced_coface = coface_it; + // else + // break; + // } + // if (coface_count == 1) { + // std::cout << "Erase ( "; + // for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base()))) + // std::cout << v << " "; + + // simplices.erase(--(reduced_coface.base())); + // //elementary_collapse(*face_it, *reduced_coface); + // std::cout << ") and then ( "; + // for (auto v: sc.simplex_vertex_range(*(--face_it.base()))) + // std::cout << v << " "; + // std::cout << ")\n"; + // simplices.erase(--((face_it++).base())); + // //face_it = simplices.rbegin(); + // //std::cout << "Size of vector: " << simplices.size() << "\n"; + // } + // else + // face_it++; + // } + // sc.initialize_filtration(); + // //std::cout << sc << std::endl; + // } + + template <class Dim_lists> + void collapse(Dim_lists& dim_lists) + { + dim_lists.collapse(); + } + +private: + + /** Collapse recursively boundary faces of the given simplex + * if its filtration is bigger than alpha_lim. + */ + void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim) + { + for (Simplex_handle face_it : sc.boundary_simplex_range()) { + + } + + } + +}; //class Weak_relaxed_witness_complex + +} // namespace witness_complex + +} // namespace Gudhi + +#endif |