diff options
author | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2016-02-11 14:18:49 +0000 |
---|---|---|
committer | skachano <skachano@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2016-02-11 14:18:49 +0000 |
commit | e7c6a8ae371de46469f0509697572c9f30071281 (patch) | |
tree | 11b2f24039a583dd6d9357a2f8b1df6e0e1bbb18 /src/Witness_complex/include | |
parent | 1a8ad726eafa817ea07a057d490db3bfc4d6ce05 (diff) |
Added relaxed witness complex + persistence example
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/relaxed-witness@1016 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 06f7fec6464e35a6ca0aa18f7a25cbb7a7f68ad5
Diffstat (limited to 'src/Witness_complex/include')
-rw-r--r-- | src/Witness_complex/include/gudhi/Relaxed_witness_complex.h | 270 |
1 files changed, 270 insertions, 0 deletions
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..6fc71a05 --- /dev/null +++ b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h @@ -0,0 +1,270 @@ +/* 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 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; + } + +}; //class Witness_complex + +} // namespace witness_complex + +} // namespace Guhdi + +#endif |